add custom track selection for cut variation for systematics

Author: Aimeric Landou

PWGJE/Tasks/trackEfficiency.cxx

@@ -13,6 +13,10 @@
 /// \author Aimeric Landou <aimeric.landou@cern.ch>
 /// \brief task that creates the histograms necessary for computation of efficiency and purity functions in offline postprocess macros; also can make mcparticle and track QC histograms
 
+
+#include "Common/Core/TrackSelection.h"
+#include "Common/Core/TrackSelectionDefaults.h"
+
 #include "PWGJE/Core/JetDerivedDataUtilities.h"
 #include "PWGJE/DataModel/Jet.h"
 #include "PWGJE/DataModel/JetReducedData.h"
@@ -84,6 +88,17 @@ struct TrackEfficiency {
   Configurable<bool> getPtHatFromHepMCXSection{"getPtHatFromHepMCXSection", true, "test configurable, configurable should be removed once well tested"};
   Configurable<bool> useTrueTrackWeight{"useTrueTrackWeight", true, "test configurable, should be set to 1 then config removed once well tested"};
 
+  // systematics variation
+  TrackSelection customTrackSelection;
+  Configurable<bool> useCustomTrackSelection{"useCustomTrackSelection", false, "whether to use the custom cuts (used for cut variation for tracking efficiency systematics)"};
+  Configurable<float> effSystMinNCrossedRowsTPC{"effSystMinNCrossedRowsTPC", 70, "min number of crossed rows TPC"};
+  Configurable<float> effSystMinNCrossedRowsOverFindableClustersTPC{"effSystMinNCrossedRowsOverFindableClustersTPC", 0.8, "min ratio of crossed rows over findable clusters TPC"};
+  Configurable<float> effSystMaxChi2PerClusterTPC{"effSystMaxChi2PerClusterTPC", 4.0, "max chi2 per cluster TPC"};
+  Configurable<float> effSystMaxChi2PerClusterITS{"effSystMaxChi2PerClusterITS", 36.0, "max chi2 per cluster ITS	"};
+  // Configurable<float> effSystMaxDcaXY{"effSystMaxDcaXY", 0.0105 * 0.035 / pT^1.1 ????, "max DCA to vertex xy"};
+  Configurable<float> effSystMaxDcaZ{"effSystMaxDcaZ", 2.0, "max DCA to vertex z"};
+  Configurable<float> effSystMinNrequiredHits{"effSystMinNrequiredHits", 1, "minimum number of hits among the 3 innermost layers of the ITS"};
+
   std::vector<int> eventSelectionBits;
   int trackSelection = -1;
 
@@ -92,6 +107,21 @@ struct TrackEfficiency {
     SplitOkCheckAnyAssocColl,      // 1
     SplitOkCheckFirstAssocCollOnly // 2
   };
+ 
+  template <typename TJetTrack>
+  bool isAcceptedTrack(TJetTrack const& jetTrack) {
+    if (!useCustomTrackSelection) {
+      if (jetderiveddatautilities::selectTrack(jetTrack, trackSelection) && jetderiveddatautilities::selectTrackDcaZ(jetTrack, trackDcaZmax)) { // if track selection is uniformTrack, dcaZ cuts need to be added as they aren't in the selection so that they can be studied here
+        return true;
+      }
+    } else {
+      const auto& aodTrack = jetTrack.template track_as<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA>>(); // might need the aodTracks to have the TracksExtra table as well; should check; check what is needed for the track selection
+      if (customTrackSelection.IsSelected(aodTrack)) {
+        return true;
+      } 
+    }
+    return false;
+  }
 
   bool isChargedParticle(int code)
   {
@@ -104,11 +134,11 @@ struct TrackEfficiency {
     return std::abs(charge) >= chargeUnit;
   }
 
-  template <typename TCollision, typename TTracks>
-  void fillTrackHistograms(TCollision const& collision, TTracks const& tracks, float weight = 1.0)
+  template <typename TCollision, typename TJetTracks>
+  void fillTrackHistograms(TCollision const& collision, TJetTracks const& jetTracks, float weight = 1.0)
   {
-    for (auto const& track : tracks) {
-      if (!(jetderiveddatautilities::selectTrack(track, trackSelection) && jetderiveddatautilities::selectTrackDcaZ(track, trackDcaZmax))) {
+    for (auto const& track : jetTracks) {
+      if (!isAcceptedTrack(track)) {
         continue;
       }
 
@@ -155,6 +185,44 @@ struct TrackEfficiency {
     eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
     trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast<std::string>(trackSelections));
 
+    if (useCustomTrackSelection) {
+      // Custom track cuts
+      LOGP(info, "Using custom track selection from values:");
+      LOGP(info, "\tminNCrossedRowsTPC= %f", effSystMinNCrossedRowsTPC.value);
+      LOGP(info, "\tminNCrossedRowsOverFindableClustersTPC= %f", effSystMinNCrossedRowsOverFindableClustersTPC.value);
+      LOGP(info, "\tmaxChi2PerClusterTPC= %f", effSystMaxChi2PerClusterTPC.value);
+      LOGP(info, "\tmaxChi2PerClusterITS= %f", effSystMaxChi2PerClusterITS.value);
+      // LOGP(info, "\tmaxDcaXY= %f", effSystMaxDcaXY.value);
+      LOGP(info, "\tmaxDcaZ= %f", effSystMaxDcaZ.value);
+      LOGP(info, "\tRequireHitsInITSLayers= %i", effSystMinNrequiredHits.value);
+
+      LOGP(info, "\trequireITS= true");
+      LOGP(info, "\trequireTPC= true");
+
+      // customTrackSelection = getGlobalTrackSelectionRun3ITSMatch(itsPattern.value); does this need to be setup?
+      LOGP(info, "Customizing track selection:");
+      int dcaSetup = 0; //default dca setup
+      customTrackSelection = getGlobalTrackSelectionRun3ITSMatch(TrackSelection::GlobalTrackRun3ITSMatching::Run3ITSibAny, dcaSetup); // takes global tracks configuration, then edit some cuts
+      customTrackSelection.SetEtaRange(-999, 999);
+      customTrackSelection.SetPtRange(0, 1e10f);
+      
+      customTrackSelection.SetMinNCrossedRowsTPC(effSystMinNCrossedRowsTPC.value);
+      customTrackSelection.SetMinNCrossedRowsOverFindableClustersTPC(effSystMinNCrossedRowsOverFindableClustersTPC.value);
+      customTrackSelection.SetMaxChi2PerClusterTPC(effSystMaxChi2PerClusterTPC.value);
+      customTrackSelection.SetMaxChi2PerClusterITS(effSystMaxChi2PerClusterITS.value);
+      // customTrackSelection.SetMaxDcaXY(effSystMaxDcaXY.value);
+      customTrackSelection.SetMaxDcaZ(effSystMaxDcaZ.value);
+      customTrackSelection.SetRequireHitsInITSLayers(effSystMinNrequiredHits.value, {0, 1, 2}); // one hit in any SPD layer (#hits, {layer0, layer1,...})
+
+      // customTrackSelection.SetRequireITSRefit(true); already set by default
+      // customTrackSelection.SetRequireTPCRefit(true); already set by default
+      // customTrackSelection.SetRequireGoldenChi2(requireGoldenChi2.value); already set by default
+
+      customTrackSelection.print();
+    } else {
+      LOGP(info, "Using standard track selection: %s", trackSelections.value);
+    }
+
     if (doprocessEFficiencyPurity || doprocessEFficiencyPurityWeighted) {
 
       registry.add("hMcCollCutsCounts", "McColl cuts count checks", {HistType::kTH1F, {{10, 0., 10.}}});
@@ -279,6 +347,22 @@ struct TrackEfficiency {
       registry.add("h2_centrality_mccollisions_weighted", "centrality vs mccollisions; centrality; collisions", {HistType::kTH2F, {centAxis, {4, 0.0, 4.0}}});
       registry.add("h2_mccollision_pthardfromweight_pthardfromhepmcxsection_weighted", "ptHard from weight vs ptHard from HepMCXSections; ptHard_weight; ptHard_hepmcxsections", {HistType::kTH2F, {{200, 0.0, 200.0}, {200, 0.0, 200.0}}});
     }
+
+    if (doprocessTrackSelectionHistograms) {
+      registry.add("h_trackselplot_tpccrossedrows", "track selection variable: number of tpc crossed rows", {HistType::kTH1F, {{165, -0.5, 164.5}}});
+      registry.add("h_trackselplot_tpccrossedrowsoverfindable", "track selection variable: ratio of of tpc crossed rows over number of findable clusters", {HistType::kTH1F, {{120, 0.0, 1.2}}});
+      registry.add("h_trackselplot_chi2ncls_tpc", "track selection variable: Chi2 / cluster for the TPC track segment", {HistType::kTH1F, {{100, 0.0, 10.0}}});
+      registry.add("h_trackselplot_chi2ncls_its", "track selection variable: Chi2 / cluster for the ITS track segment", {HistType::kTH1F, {{200, 0.0, 40.0}}});
+      registry.add("h_trackselplot_dcaxy", "track selection variable: dca XY", {HistType::kTH1F, {{1000, -1.0, 1.0}}});
+      registry.add("h_trackselplot_dcaz", "track selection variable: dca Z", {HistType::kTH1F, {{1000, -1.0, 1.0}}});
+
+      registry.add("h2_trackselplot_pt_tpccrossedrows", "track selection variable: pt vs number of tpc crossed rows", {HistType::kTH2F, {{200, 0., 200.},{165, -0.5, 164.5}}});
+      registry.add("h2_trackselplot_pt_tpccrossedrowsoverfindable", "track selection variable: pt vs ratio of of tpc crossed rows over number of findable clusters", {HistType::kTH2F, {{200, 0., 200.}, {120, 0.0, 1.2}}});
+      registry.add("h2_trackselplot_pt_chi2ncls_tpc", "track selection variable: pt vs Chi2 / cluster for the TPC track segment", {HistType::kTH2F, {{200, 0., 200.}, {100, 0.0, 10.0}}});
+      registry.add("h2_trackselplot_pt_chi2ncls_its", "track selection variable: pt vs Chi2 / cluster for the ITS track segment", {HistType::kTH2F, {{200, 0., 200.}, {200, 0.0, 40.0}}});
+      registry.add("h2_trackselplot_pt_dcaxy", "track selection variable: pt vs dca XY", {HistType::kTH2F, {{200, 0., 200.}, {1000, -1.0, 1.0}}});
+      registry.add("h2_trackselplot_pt_dcaz", "track selection variable: pt vs dca Z", {HistType::kTH2F, {{200, 0., 200.}, {1000, -1.0, 1.0}}});
+    }
   }
 
   Preslice<aod::JetTracksMCD> tracksPerJCollision = o2::aod::jtrack::collisionId;
@@ -291,7 +375,8 @@ struct TrackEfficiency {
   void processEFficiencyPurity(soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>::iterator const& mcCollision,
                                soa::Join<aod::McCollisions, aod::HepMCXSections> const&,
                                soa::SmallGroups<aod::JetCollisionsMCD> const& collisions, // smallgroups gives only the collisions associated to the current mccollision, thanks to the mccollisionlabel pre-integrated in jetcollisionsmcd
-                               soa::Join<aod::JetTracksMCD, aod::JTrackExtras> const& jetTracks,
+                               soa::Join<aod::JetTracksMCD, aod::JTrackExtras, aod::JTrackPIs> const& jetTracks,
+                               soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA> const&,
                                JetParticlesWithOriginal const& jMcParticles)
   {
     // missing:
@@ -407,7 +492,7 @@ struct TrackEfficiency {
       for (auto const& track : collTracks) {
         registry.fill(HIST("hTrackCutsCounts"), 0.5);
 
-        if (!(jetderiveddatautilities::selectTrack(track, trackSelection) && jetderiveddatautilities::selectTrackDcaZ(track, trackDcaZmax))) { // if track selection is uniformTrack, dcaZ cuts need to be added as they aren't in the selection so that they can be studied here
+        if (!isAcceptedTrack(track)) {
           continue;
         }
         registry.fill(HIST("hTrackCutsCounts"), 1.5);
@@ -472,7 +557,8 @@ struct TrackEfficiency {
   void processEFficiencyPurityWeighted(soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>::iterator const& mcCollision,
                                        soa::Join<aod::McCollisions, aod::HepMCXSections> const&,
                                        soa::SmallGroups<aod::JetCollisionsMCD> const& collisions, // smallgroups gives only the collisions associated to the current mccollision, thanks to the mccollisionlabel pre-integrated in jetcollisionsmcd
-                                       soa::Join<aod::JetTracksMCD, aod::JTrackExtras> const& jetTracks,
+                                       soa::Join<aod::JetTracksMCD, aod::JTrackExtras, aod::JTrackPIs> const& jetTracks,
+                                       soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA> const&,
                                        JetParticlesWithOriginal const& jMcParticles)
   {
     // missing:
@@ -576,7 +662,7 @@ struct TrackEfficiency {
       for (auto const& track : collTracks) {
         registry.fill(HIST("hTrackCutsCounts"), 0.5, mcCollision.weight());
 
-        if (!(jetderiveddatautilities::selectTrack(track, trackSelection) && jetderiveddatautilities::selectTrackDcaZ(track, trackDcaZmax))) { // if track selection is uniformTrack, dcaZ cuts need to be added as they aren't in the selection so that they can be studied here
+        if (!isAcceptedTrack(track)) {
           continue;
         }
         registry.fill(HIST("hTrackCutsCounts"), 1.5, mcCollision.weight());
@@ -595,8 +681,8 @@ struct TrackEfficiency {
         registry.fill(HIST("hTrackCutsCounts"), 3.5, mcCollision.weight());
 
         auto mcParticle = track.mcParticle_as<JetParticlesWithOriginal>();
-        auto trueTrackMcCollision = mcParticle.mcCollision_as<aod::JetMcCollisions>();
-        float trueTrackCollEventWeight = useTrueTrackWeight ? trueTrackMcCollision.weight() : mcCollEventWeight; // test1
+        auto trueTrackMcCollision = mcParticle.mcCollision_as<soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>>();
+        float trueTrackCollEventWeight = useTrueTrackWeight ? trueTrackMcCollision.weight() : mcCollEventWeight;
 
         auto jMcParticleFromTrack = track.mcParticle_as<JetParticlesWithOriginal>();
         if (!jMcParticleFromTrack.isPhysicalPrimary()) {
@@ -648,7 +734,8 @@ struct TrackEfficiency {
   PROCESS_SWITCH(TrackEfficiency, processEFficiencyPurityWeighted, "Histograms for efficiency and purity quantities for weighted simulations", false);
 
   void processTracksFromData(soa::Filtered<aod::JetCollisions>::iterator const& collision,
-                             soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras>> const& tracks)
+                             soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras, aod::JTrackPIs>> const& jetTracks,
+                             soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA> const&)
   {
     if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents)) {
       return;
@@ -657,14 +744,15 @@ struct TrackEfficiency {
       return;
     }
 
-    fillTrackHistograms(collision, tracks);
+    fillTrackHistograms(collision, jetTracks);
   }
   PROCESS_SWITCH(TrackEfficiency, processTracksFromData, "QA for charged tracks in data", false);
 
-  void processTracksFromMc(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>::iterator const& collision,
+  void processTracksFromMc(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>::iterator const& collision, // a filter should probably be added here to stay consistent with processTracksFromData
                            soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs> const&,
                            soa::Join<aod::McCollisions, aod::HepMCXSections> const&,
-                           soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras>> const& tracks)
+                           soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras, aod::JTrackPIs>> const& jetTracks,
+                           soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA> const&)
   {
     if (!collision.has_mcCollision()) { // the collision is fake and has no associated mc coll; skip as .mccollision() cannot be called
       return;
@@ -681,19 +769,20 @@ struct TrackEfficiency {
       return;
     }
 
-    fillTrackHistograms(collision, tracks);
+    fillTrackHistograms(collision, jetTracks);
   }
   PROCESS_SWITCH(TrackEfficiency, processTracksFromMc, "QA for charged tracks in MC without weights", false);
 
-  void processTracksFromMcWeighted(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>::iterator const& collision,
+  void processTracksFromMcWeighted(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>::iterator const& collision, // a filter should probably be added here to stay consistent with processTracksFromData
                                    soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs> const&,
                                    soa::Join<aod::McCollisions, aod::HepMCXSections> const&,
-                                   soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras>> const& tracks)
+                                   soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras, aod::JTrackPIs>> const& jetTracks,
+                                   soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA> const&)
   {
     if (!collision.has_mcCollision()) { // the collision is fake and has no associated mc coll; skip as .mccollision() cannot be called
       return;
     }
-    float eventWeight = collision.mcCollision().weight();
+    float eventWeight = collision.mcCollision_as<soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>>().weight();
     if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents)) {
       return;
     }
@@ -706,7 +795,7 @@ struct TrackEfficiency {
       return;
     }
 
-    fillTrackHistograms(collision, tracks, eventWeight);
+    fillTrackHistograms(collision, jetTracks, eventWeight);
   }
   PROCESS_SWITCH(TrackEfficiency, processTracksFromMcWeighted, "QA for charged tracks in weighted MC", false);
 
@@ -836,7 +925,7 @@ struct TrackEfficiency {
   }
   PROCESS_SWITCH(TrackEfficiency, processCollisionsFromData, "QA for reconstructed collisions in data", false);
 
-  void processCollisionsFromMc(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>::iterator const& collision,
+  void processCollisionsFromMc(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>::iterator const& collision, // a filter should probably be added here to stay consistent with processTracksFromData
                                soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs> const&,
                                soa::Join<aod::McCollisions, aod::HepMCXSections> const&)
   {
@@ -866,15 +955,15 @@ struct TrackEfficiency {
   }
   PROCESS_SWITCH(TrackEfficiency, processCollisionsFromMc, "QA for reconstructed collisions in MC without weights", false);
 
-  void processCollisionsFromMcWeighted(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>::iterator const& collision,
+  void processCollisionsFromMcWeighted(soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>::iterator const& collision, // a filter should probably be added here to stay consistent with processTracksFromData
                                        soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs> const&,
                                        soa::Join<aod::McCollisions, aod::HepMCXSections> const&)
   {
     if (!collision.has_mcCollision()) { // the collision is fake and has no associated mc coll; skip as .mccollision() cannot be called
       registry.fill(HIST("h_fakecollisions"), 0.5);
       return;
     }
-    float eventWeight = collision.mcCollision().weight();
+    float eventWeight = collision.mcCollision_as<soa::Join<aod::JetMcCollisions, aod::JMcCollisionPIs>>().weight();
     registry.fill(HIST("h_collisions"), 0.5);
     registry.fill(HIST("h_collisions_weighted"), 0.5, eventWeight);
     if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents)) {
@@ -1015,6 +1104,23 @@ struct TrackEfficiency {
     registry.fill(HIST("h_mccollisions_weighted"), 2.5, eventWeight);
   }
   PROCESS_SWITCH(TrackEfficiency, processMcCollisionsWeighted, "QA for McCollisions in weighted MC", false);
+
+  void processTrackSelectionHistograms(soa::Join<aod::Tracks, aod::TracksExtra, o2::aod::TracksDCA>::iterator const& track) { // should probably select collisions 
+    registry.fill(HIST("h_trackselplot_tpccrossedrows"), track.tpcNClsCrossedRows());
+    registry.fill(HIST("h_trackselplot_tpccrossedrowsoverfindable"), track.tpcCrossedRowsOverFindableCls());
+    registry.fill(HIST("h_trackselplot_chi2ncls_tpc"), track.tpcChi2NCl	());
+    registry.fill(HIST("h_trackselplot_chi2ncls_its"), track.itsChi2NCl	());
+    registry.fill(HIST("h_trackselplot_dcaxy"), track.dcaXY());
+    registry.fill(HIST("h_trackselplot_dcaz"), track.dcaZ());
+
+    registry.fill(HIST("h2_trackselplot_pt_tpccrossedrows"), track.pt(), track.tpcNClsCrossedRows());
+    registry.fill(HIST("h2_trackselplot_pt_tpccrossedrowsoverfindable"), track.pt(), track.tpcCrossedRowsOverFindableCls());
+    registry.fill(HIST("h2_trackselplot_pt_chi2ncls_tpc"), track.pt(), track.tpcChi2NCl());
+    registry.fill(HIST("h2_trackselplot_pt_chi2ncls_its"), track.pt(), track.itsChi2NCl());
+    registry.fill(HIST("h2_trackselplot_pt_dcaxy"), track.pt(), track.dcaXY());
+    registry.fill(HIST("h2_trackselplot_pt_dcaz"), track.pt(), track.dcaZ());
+  }
+  PROCESS_SWITCH(TrackEfficiency, processTrackSelectionHistograms, "plots distributions of variables that are cut on during track selection", false);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)