[PWGJE] Adding Outlier rejection task for MC

PWGJE/DataModel/JetReducedData.h

@@ -1,4 +1,4 @@
 // Copyright 2019-2020 CERN and copyright holders of ALICE O2.
 // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 // All rights not expressly granted are reserved.
 //
@@ -12,7 +12,7 @@
 ///
 /// \brief Table definitions for reduced data model for jets
 ///
 /// \author Nima Zardoshti <nima.zardoshti@cern.ch>
 
 #ifndef PWGJE_DATAMODEL_JETREDUCEDDATA_H_
 #define PWGJE_DATAMODEL_JETREDUCEDDATA_H_
@@ -32,7 +32,7 @@
 
 namespace jbc
 {
 DECLARE_SOA_INDEX_COLUMN(BC, bc);
 DECLARE_SOA_COLUMN(RunNumber, runNumber, int);
 DECLARE_SOA_COLUMN(GlobalBC, globalBC, uint64_t);
 DECLARE_SOA_COLUMN(Timestamp, timestamp, uint64_t);
@@ -67,7 +67,7 @@
 namespace jcollision
 {
 DECLARE_SOA_INDEX_COLUMN(Collision, collision);
 DECLARE_SOA_INDEX_COLUMN(JBC, bc);
 DECLARE_SOA_COLUMN(PosX, posX, float);
 DECLARE_SOA_COLUMN(PosY, posY, float);
 DECLARE_SOA_COLUMN(PosZ, posZ, float);
@@ -98,7 +98,8 @@
 DECLARE_SOA_COLUMN(ReadCountsWithTVXAndZVertexAndSel7KINT7, readCountsWithTVXAndZVertexAndSel7KINT7, std::vector<int>);
 DECLARE_SOA_COLUMN(ReadCountsWithCustom, readCountsWithCustom, std::vector<int>);
 DECLARE_SOA_COLUMN(IsAmbiguous, isAmbiguous, bool);
 DECLARE_SOA_COLUMN(IsEMCALReadout, isEmcalReadout, bool);
+DECLARE_SOA_COLUMN(IsOutlier, isOutlier, bool);
 } // namespace jcollision
 
 DECLARE_SOA_TABLE_STAGED(JCollisions, "JCOLLISION",
@@ -122,6 +123,9 @@
                          jcollision::Weight,
                          jcollision::SubGeneratorId);
 
+DECLARE_SOA_TABLE_STAGED(JCollisionOutliers, "JCOLLISIONOUTLR",
+                         jcollision::IsOutlier);
+
 DECLARE_SOA_TABLE_STAGED(JEMCCollisionLbs, "JEMCCOLLISIONLB",
                          jcollision::IsAmbiguous,
                          jcollision::IsEMCALReadout);
@@ -170,6 +174,8 @@
 DECLARE_SOA_COLUMN(Attempted, attempted, uint64_t);
 DECLARE_SOA_COLUMN(XsectGen, xsectGen, float);
 DECLARE_SOA_COLUMN(XsectErr, xsectErr, float);
+DECLARE_SOA_COLUMN(PtHard, ptHard, float);
+DECLARE_SOA_COLUMN(IsOutlier, isOutlier, bool);
 } // namespace jmccollision
 DECLARE_SOA_TABLE_STAGED(JMcCollisions, "JMCCOLLISION",
                          o2::soa::Index<>,
@@ -181,7 +187,8 @@
                          jmccollision::Accepted,
                          jmccollision::Attempted,
                          jmccollision::XsectGen,
-                         jmccollision::XsectErr);
+                         jmccollision::XsectErr,
+                         jmccollision::PtHard);
 
 using JMcCollision = JMcCollisions::iterator;
 using StoredJMcCollision = StoredJMcCollisions::iterator;
@@ -191,22 +198,25 @@
 
 namespace jmccollisionlb
 {
 DECLARE_SOA_INDEX_COLUMN(JMcCollision, mcCollision);
 }
 
 DECLARE_SOA_TABLE_STAGED(JMcCollisionLbs, "JMCCOLLISIONLB",
                          jmccollisionlb::JMcCollisionId);
 
+DECLARE_SOA_TABLE_STAGED(JMcCollisionOutliers, "JMCCOLLISIONOUTLR",
+                         jmccollision::IsOutlier);
+
 namespace jtrack
 {
 DECLARE_SOA_INDEX_COLUMN(JCollision, collision);
 DECLARE_SOA_INDEX_COLUMN(Track, track);
 DECLARE_SOA_COLUMN(Pt, pt, float);
 DECLARE_SOA_COLUMN(Eta, eta, float);
 DECLARE_SOA_COLUMN(Phi, phi, float);
 DECLARE_SOA_COLUMN(DCAX, dcaX, float);
 DECLARE_SOA_COLUMN(DCAY, dcaY, float);
 DECLARE_SOA_COLUMN(DCAZ, dcaZ, float);
 DECLARE_SOA_COLUMN(DCAXY, dcaXY, float);
 DECLARE_SOA_COLUMN(DCAXYZ, dcaXYZ, float);
 DECLARE_SOA_COLUMN(SigmaDCAZ, sigmadcaZ, float);

PWGJE/TableProducer/CMakeLists.txt

@@ -53,6 +53,11 @@ o2physics_add_dpl_workflow(jet-eventweight-mcp
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(mc-outlier-rejector
+                    SOURCES mcOutlierRejector.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(jet-track-derived
                     SOURCES jetTrackDerived.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore

PWGJE/TableProducer/derivedDataProducer.cxx

@@ -256,14 +256,14 @@ struct JetDerivedDataProducerTask {
 
   void processMcCollisions(aod::McCollision const& mcCollision)
   {
-    products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), 1, 1, 1.0, 1.0);
+    products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), 1, 1, 1.0, 1.0, 999.0);
     products.jMcCollisionsParentIndexTable(mcCollision.globalIndex());
   }
   PROCESS_SWITCH(JetDerivedDataProducerTask, processMcCollisions, "produces derived MC collision table", false);
 
   void processMcCollisionsWithXsection(soa::Join<aod::McCollisions, aod::HepMCXSections>::iterator const& mcCollision)
   {
-    products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr());
+    products.jMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.getSubGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr(), mcCollision.ptHard());
     products.jMcCollisionsParentIndexTable(mcCollision.globalIndex());
   }
   PROCESS_SWITCH(JetDerivedDataProducerTask, processMcCollisionsWithXsection, "produces derived MC collision table with cross section information", false);

PWGJE/TableProducer/derivedDataWriter.cxx

@@ -432,7 +432,7 @@ struct JetDerivedDataWriter {
     mcCollisionMapping.resize(mcCollisions.size(), -1);
     for (auto const& mcCollision : mcCollisions) {
       if (mcCollision.isMcCollisionSelected()) {
-        products.storedJMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.subGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr());
+        products.storedJMcCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.weight(), mcCollision.subGeneratorId(), mcCollision.accepted(), mcCollision.attempted(), mcCollision.xsectGen(), mcCollision.xsectErr(), mcCollision.ptHard());
         products.storedJMcCollisionsParentIndexTable(mcCollision.mcCollisionId());
         mcCollisionMapping[mcCollision.globalIndex()] = products.storedJMcCollisionsTable.lastIndex();
       }

PWGJE/TableProducer/mcOutlierRejector.cxx

@@ -0,0 +1,132 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+// Task to produce a table joinable to the jcollision table which contains an outlier rejector
+//
+/// \author Nima Zardoshti <nima.zardoshti@cern.ch>
+
+#include "PWGJE/DataModel/Jet.h"
+#include "PWGJE/DataModel/JetReducedData.h"
+
+#include "Framework/ASoA.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/O2DatabasePDGPlugin.h"
+#include <Framework/AnalysisHelpers.h>
+#include <Framework/Configurable.h>
+#include <Framework/InitContext.h>
+#include <Framework/runDataProcessing.h>
+
+#include <vector>
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+
+struct McOutlierRejectorTask {
+  Produces<aod::JCollisionOutliers> collisionOutliers;
+  Produces<aod::JMcCollisionOutliers> mcCollisionOutliers;
+
+  Configurable<bool> checkmcCollisionForCollision{"checkmcCollisionForCollision", true, "additionally reject collision based on mcCollision"};
+  Configurable<float> ptHatMax{"ptHatMax", 4.0, "maximum factor of pt hat the leading jet in the event is allowed"};
+
+  std::vector<bool> collisionFlag;
+  std::vector<bool> mcCollisionFlag;
+
+  void processSetupCollisionSelection(aod::JCollisions const& collisions)
+  {
+    collisionFlag.clear();
+    collisionFlag.resize(collisions.size(), false);
+  }
+  PROCESS_SWITCH(McOutlierRejectorTask, processSetupCollisionSelection, "Setup collision processing", true);
+
+  void processSetupMcCollisionSelection(aod::JMcCollisions const& mcCollisions)
+  {
+    mcCollisionFlag.clear();
+    mcCollisionFlag.resize(mcCollisions.size(), false);
+  }
+  PROCESS_SWITCH(McOutlierRejectorTask, processSetupMcCollisionSelection, "Setup MC Collision processing", true);
+
+  template <typename T>
+  void collisionSelection(int32_t collisionIndex, T const& selectionObjects, float ptHard, std::vector<bool>& flagArray)
+  {
+
+    if (selectionObjects.size() != 0) {
+      float maxSelectionObjectPt = 0.0;
+      if constexpr (std::is_same_v<std::decay_t<T>, aod::JetTracks> || std::is_same_v<std::decay_t<T>, aod::JetParticles>) {
+        for (auto selectionObject : selectionObjects) {
+          if (selectionObject.pt() > maxSelectionObjectPt) {
+            maxSelectionObjectPt = selectionObject.pt();
+          }
+        }
+      } else {
+        maxSelectionObjectPt = selectionObjects.iteratorAt(0).pt();
+      }
+
+      if (maxSelectionObjectPt > ptHatMax * ptHard) {
+        flagArray[collisionIndex] = true; // Currently if running multiple different jet finders, then a single type of jet can veto an event for others. Decide if this is the best way
+      }
+    }
+  }
+
+  template <typename T>
+  void processSelectionObjects(aod::JetCollisionMCD const& collision, T const& selectionObjects, aod::JetMcCollisions const&)
+  {
+    auto mcCollision = collision.mcCollision_as<aod::JetMcCollisions>();
+    collisionSelection(collision.globalIndex(), selectionObjects, mcCollision.ptHard(), collisionFlag);
+  }
+
+  template <typename T>
+  void processSelectionMcObjects(aod::JetMcCollision const& mcCollision, T const& selectionMcObjects)
+  {
+    collisionSelection(mcCollision.globalIndex(), selectionMcObjects, mcCollision.ptHard(), mcCollisionFlag);
+  }
+
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::ChargedMCDetectorLevelJets>, processSelectingChargedMCDetectorLevelJets, "process mc detector level charged jets", true);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::NeutralMCDetectorLevelJets>, processSelectingNeutralMCDetectorLevelJets, "process mc detector level neutral jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::FullMCDetectorLevelJets>, processSelectingFullMCDetectorLevelJets, "process mc detector level full jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::D0ChargedMCDetectorLevelJets>, processSelectingD0ChargedMCDetectorLevelJets, "process mc detector level D0 charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::DplusChargedMCDetectorLevelJets>, processSelectingDplusChargedMCDetectorLevelJets, "process mc detector level Dplus charged jets", false);
+  // PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::DstarChargedMCDetectorLevelJets>, processSelectingDstarChargedMCDetectorLevelJets, "process mc detector level Dstar charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::LcChargedMCDetectorLevelJets>, processSelectingLcChargedMCDetectorLevelJets, "process mc detector level Lc charged jets", false);
+  // PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::B0ChargedMCDetectorLevelJets>, processSelectingB0ChargedMCDetectorLevelJets, "process mc detector level B0 charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::BplusChargedMCDetectorLevelJets>, processSelectingBplusChargedMCDetectorLevelJets, "process mc detector level Bplus charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::DielectronChargedMCDetectorLevelJets>, processSelectingDielectronChargedMCDetectorLevelJets, "process mc detector level Dielectron charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::JetTracks>, processSelectingTracks, "process tracks", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::ChargedMCParticleLevelJets>, processSelectingChargedMCParticleLevelJets, "process mc particle level charged jets", true);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::NeutralMCParticleLevelJets>, processSelectingNeutralMCParticleLevelJets, "process mc particle level neutral jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::FullMCParticleLevelJets>, processSelectingFullMCParticleLevelJets, "process mc particle level full jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::D0ChargedMCParticleLevelJets>, processSelectingD0ChargedMCParticleLevelJets, "process mc particle level D0 charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::DplusChargedMCParticleLevelJets>, processSelectingDplusChargedMCParticleLevelJets, "process mc particle level Dplus charged jets", false);
+  // PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::DstarChargedMCParticleLevelJets>, processSelectingDstarChargedMCParticleLevelJets, "process mc particle level Dstar charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::LcChargedMCParticleLevelJets>, processSelectingLcChargedMCParticleLevelJets, "process mc particle level Lc charged jets", false);
+  // PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::B0ChargedMCParticleLevelJets>, processSelectingB0ChargedMCParticleLevelJets, "process mc particle level B0 charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::BplusChargedMCParticleLevelJets>, processSelectingBplusChargedMCParticleLevelJets, "process mc particle level Bplus charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::DielectronChargedMCParticleLevelJets>, processSelectingDielectronChargedMCParticleLevelJets, "process mc particle level Dielectron charged jets", false);
+  PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionMcObjects<aod::JetParticles>, processSelectingParticles, "process mc particles", false);
+
+  void processStoreCollisionDecision(aod::JetCollisionMCD const& collision)
+  {
+    bool rejectCollision = collisionFlag[collision.globalIndex()];
+    if (!rejectCollision && checkmcCollisionForCollision && collision.has_mcCollision()) {
+      rejectCollision = mcCollisionFlag[collision.mcCollisionId()];
+    }
+    collisionOutliers(rejectCollision);
+  }
+  PROCESS_SWITCH(McOutlierRejectorTask, processStoreCollisionDecision, "write out decision of rejecting collision", true);
+
+  void processStoreMcCollisionDecision(aod::JetMcCollision const& mcCollision)
+  {
+    mcCollisionOutliers(mcCollisionFlag[mcCollision.globalIndex()]);
+  }
+  PROCESS_SWITCH(McOutlierRejectorTask, processStoreMcCollisionDecision, "write out decision of rejecting mcCollision", true);
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<McOutlierRejectorTask>(cfgc, TaskName{"mc-outlier-rejector"})}; }