[PWGJE] include origin of track in outlier ID (#13929)

Author: jaimenorman

PWGJE/TableProducer/mcOutlierRejector.cxx

@@ -33,6 +33,7 @@ using namespace o2::framework::expressions;
 struct McOutlierRejectorTask {
   Produces<aod::JCollisionOutliers> collisionOutliers;
   Produces<aod::JMcCollisionOutliers> mcCollisionOutliers;
+  Preslice<soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>> perColParticle = aod::jmccollision::mcCollisionId;
 
   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"};
@@ -55,22 +56,36 @@ struct McOutlierRejectorTask {
   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)
+  void collisionSelection(int32_t collisionIndex, int32_t mcCollisionId, T const& selectionObjects, float ptHard, std::vector<bool>& flagArray, std::optional<std::reference_wrapper<const soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>>> mcCollisionsOpt = std::nullopt)
   {
-
     if (selectionObjects.size() != 0) {
+      bool isTrueOutlier = true;
       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>) {
+      if constexpr (std::is_same_v<std::decay_t<T>, aod::JetTracksMCD> || std::is_same_v<std::decay_t<T>, aod::JetParticles>) {
         for (auto selectionObject : selectionObjects) {
           if (selectionObject.pt() > maxSelectionObjectPt) {
             maxSelectionObjectPt = selectionObject.pt();
+            // may be slow - could save only MC particle then check difference only for tracks IDd as outliers?
+            if constexpr (std::is_same_v<std::decay_t<T>, aod::JetTracksMCD>) {
+              auto& mcCollisions = mcCollisionsOpt.value().get();
+              auto mcParticle = selectionObject.template mcParticle_as<soa::Join<aod::JetParticles, aod::JMcParticlePIs>>();
+              auto mcCollision = mcCollisions.sliceBy(perColParticle, mcParticle.mcCollisionId());
+              int subGenID = mcCollision.begin().subGeneratorId();
+              int diffCollisionID = mcParticle.mcCollisionId() - mcCollisionId;
+              if (subGenID == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap &&
+                  diffCollisionID != 0) {
+                isTrueOutlier = true;
+              } else {
+                isTrueOutlier = false;
+              }
+            }
           }
         }
       } else {
         maxSelectionObjectPt = selectionObjects.iteratorAt(0).pt();
       }
 
-      if (maxSelectionObjectPt > ptHatMax * ptHard) {
+      if (maxSelectionObjectPt > ptHatMax * ptHard && isTrueOutlier) {
         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
       }
     }
@@ -80,13 +95,20 @@ struct McOutlierRejectorTask {
   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);
+    collisionSelection(collision.globalIndex(), 1., selectionObjects, mcCollision.ptHard(), collisionFlag);
+  }
+
+  template <typename T>
+  void processSelectionObjectsTracks(aod::JetCollisionMCD const& collision, T const& selectionObjects, soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs> const& mcCollisions, soa::Join<aod::JetParticles, aod::JMcParticlePIs> const&)
+  {
+    auto mcCollision = collision.mcCollision_as<soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>>();
+    collisionSelection(collision.globalIndex(), collision.mcCollisionId(), selectionObjects, mcCollision.ptHard(), collisionFlag, mcCollisions);
   }
 
   template <typename T>
   void processSelectionMcObjects(aod::JetMcCollision const& mcCollision, T const& selectionMcObjects)
   {
-    collisionSelection(mcCollision.globalIndex(), selectionMcObjects, mcCollision.ptHard(), mcCollisionFlag);
+    collisionSelection(mcCollision.globalIndex(), 1., selectionMcObjects, mcCollision.ptHard(), mcCollisionFlag);
   }
 
   PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::ChargedMCDetectorLevelJets>, processSelectingChargedMCDetectorLevelJets, "process mc detector level charged jets", true);
@@ -101,7 +123,7 @@ struct McOutlierRejectorTask {
   PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::BplusChargedMCDetectorLevelJets>, processSelectingBplusChargedMCDetectorLevelJets, "process mc detector level Bplus charged jets", false);
   PROCESS_SWITCH_FULL(McOutlierRejectorTask, processSelectionObjects<aod::XicToXiPiPiChargedMCDetectorLevelJets>, processSelectingXicToXiPiPiChargedMCDetectorLevelJets, "process mc detector level XicToXiPiPi 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, processSelectionObjectsTracks<aod::JetTracksMCD>, 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);

PWGJE/Tasks/jetFinderQA.cxx

@@ -329,14 +329,14 @@ struct JetFinderQATask {
     }
 
     if (doprocessTracks || doprocessTracksWeighted) {
-      registry.add("h_collisions", "event status;event status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}});
+      registry.add("h_collisions", "event status;event status;entries", {HistType::kTH1F, {{5, 0.0, 5.0}}});
       registry.add("h2_centrality_collisions", "centrality vs collisions; centrality; collisions", {HistType::kTH2F, {{1200, -10.0, 110.0}, {4, 0.0, 4.0}}});
       registry.add("h3_centrality_track_pt_track_phi", "centrality vs track pT vs track #varphi; centrality; #it{p}_{T,track} (GeV/#it{c}); #varphi_{track}", {HistType::kTH3F, {{1200, -10.0, 110.0}, {200, 0., 200.}, {160, -1.0, 7.}}});
       registry.add("h3_centrality_track_pt_track_eta", "centrality vs track pT vs track #eta; centrality; #it{p}_{T,track} (GeV/#it{c}); #eta_{track}", {HistType::kTH3F, {{1200, -10.0, 110.0}, {200, 0., 200.}, trackEtaAxis}});
       registry.add("h3_centrality_track_pt_track_dcaxy", "centrality vs track pT vs track DCA_{xy}; centrality; #it{p}_{T,track} (GeV/#it{c}); track DCA_{xy}", {HistType::kTH3F, {{120, -10.0, 110.0}, {20, 0., 100.}, {200, -0.15, 0.15}}});
       registry.add("h3_track_pt_track_eta_track_phi", "track pT vs track #eta vs track #varphi; #it{p}_{T,track} (GeV/#it{c}); #eta_{track}; #varphi_{track}", {HistType::kTH3F, {{200, 0., 200.}, trackEtaAxis, {160, -1.0, 7.}}});
       if (doprocessTracksWeighted) {
-        registry.add("h_collisions_weighted", "event status;event status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}});
+        registry.add("h_collisions_weighted", "event status;event status;entries", {HistType::kTH1F, {{5, 0.0, 5.0}}});
       }
     }
     if (doprocessTracksSub) {
@@ -904,11 +904,14 @@ struct JetFinderQATask {
   }
   PROCESS_SWITCH(JetFinderQATask, processJetsMCD, "jet finder QA mcd", false);
 
-  void processJetsMCDWeighted(soa::Filtered<aod::JetCollisions>::iterator const& collision, soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetEventWeights> const& jets, aod::JetTracks const&)
+  void processJetsMCDWeighted(soa::Filtered<soa::Join<aod::JetCollisions, aod::JCollisionOutliers>>::iterator const& collision, soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetEventWeights> const& jets, aod::JetTracks const&)
   {
     if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {
       return;
     }
+    if (collision.isOutlier()) {
+      return;
+    }
     for (auto const& jet : jets) {
       if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
         continue;
@@ -946,7 +949,7 @@ struct JetFinderQATask {
   }
   PROCESS_SWITCH(JetFinderQATask, processJetsMCP, "jet finder QA mcp", false);
 
-  void processJetsMCPWeighted(soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetEventWeights>::iterator const& jet, aod::JetParticles const&, aod::JetMcCollisions const&, soa::Filtered<aod::JetCollisionsMCD> const& collisions)
+  void processJetsMCPWeighted(soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetEventWeights>::iterator const& jet, aod::JetParticles const&, aod::JetMcCollisions const&, soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::JCollisionOutliers>> const& collisions)
   {
     if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
       return;
@@ -962,7 +965,7 @@ struct JetFinderQATask {
     }
     if (checkMcCollisionIsMatched) {
       auto collisionspermcpjet = collisions.sliceBy(CollisionsPerMCPCollision, jet.mcCollisionId());
-      if (collisionspermcpjet.size() >= 1 && jetderiveddatautilities::selectCollision(collisionspermcpjet.begin(), eventSelectionBits)) {
+      if (collisionspermcpjet.size() >= 1 && jetderiveddatautilities::selectCollision(collisionspermcpjet.begin(), eventSelectionBits) && !collisionspermcpjet.begin().isOutlier()) {
         fillMCPHistograms(jet, jet.eventWeight());
       }
     } else {
@@ -997,14 +1000,17 @@ struct JetFinderQATask {
   }
   PROCESS_SWITCH(JetFinderQATask, processJetsMCPMCDMatched, "jet finder QA matched mcp and mcd", false);
 
-  void processJetsMCPMCDMatchedWeighted(soa::Filtered<aod::JetCollisions>::iterator const& collision,
+  void processJetsMCPMCDMatchedWeighted(soa::Filtered<soa::Join<aod::JetCollisions, aod::JCollisionOutliers>>::iterator const& collision,
                                         soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets, aod::ChargedMCDetectorLevelJetEventWeights> const& mcdjets,
                                         soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetsMatchedToChargedMCDetectorLevelJets, aod::ChargedMCParticleLevelJetEventWeights> const&,
                                         aod::JetTracks const&, aod::JetParticles const&)
   {
     if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {
       return;
     }
+    if (collision.isOutlier()) {
+      return;
+    }
     for (const auto& mcdjet : mcdjets) {
       if (!jetfindingutilities::isInEtaAcceptance(mcdjet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
         continue;
@@ -1188,7 +1194,7 @@ struct JetFinderQATask {
   }
   PROCESS_SWITCH(JetFinderQATask, processTracks, "QA for charged tracks", false);
 
-  void processTracksWeighted(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>::iterator const& collision,
+  void processTracksWeighted(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs, aod::JCollisionOutliers>::iterator const& collision,
                              aod::JetMcCollisions const&,
                              soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras>> const& tracks)
   {
@@ -1208,6 +1214,11 @@ struct JetFinderQATask {
     }
     registry.fill(HIST("h_collisions"), 2.5);
     registry.fill(HIST("h_collisions_weighted"), 2.5, eventWeight);
+    if (collision.isOutlier()) {
+      return;
+    }
+    registry.fill(HIST("h_collisions"), 3.5);
+    registry.fill(HIST("h_collisions_weighted"), 3.5, eventWeight);
     fillTrackHistograms(collision, tracks, eventWeight);
   }
   PROCESS_SWITCH(JetFinderQATask, processTracksWeighted, "QA for charged tracks weighted", false);