[PWGJE] trackEfficiency: add weights to eff. and purity plots

Author: aimeric-landou

PWGJE/Tasks/trackEfficiency.cxx

@@ -76,6 +76,8 @@ struct TrackEfficiencyJets {
   Configurable<std::vector<double>> centralityBinning{"centralityBinning", {0., 10., 50., 70.}, "binning of centrality histograms"};
   Configurable<int> intRateNBins{"intRateNBins", 50, "number of bins for interaction rate axis"};
   Configurable<float> intRateMax{"intRateMax", 50000.0, "maximum value of interaction rate axis"};
+  Configurable<int> phiEffNBins{"phiEffNBins", 200, "number of bins for phi axis in efficiency plots"};
+  Configurable<int> etaEffNBins{"etaEffNBins", 200, "number of bins for eta axis in efficiency plots"};
 
   std::vector<int> eventSelectionBits;
   int trackSelection = -1;
@@ -126,7 +128,7 @@ struct TrackEfficiencyJets {
     eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
     trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast<std::string>(trackSelections));
 
-    if (doprocessEFficiencyPurity) {
+    if (doprocessEFficiencyPurity || doprocessEFficiencyPurityWeighted) {
 
       registry.add("hMcCollCutsCounts", "McColl cuts count checks", {HistType::kTH1F, {{10, 0., 10.}}});
       registry.get<TH1>(HIST("hMcCollCutsCounts"))->GetXaxis()->SetBinLabel(1, "allMcColl");
@@ -151,8 +153,8 @@ struct TrackEfficiencyJets {
 
       AxisSpec ptAxis_eff = {nBinsLowPt, 0., 10., "#it{p}_{T} (GeV/#it{c})"};
       AxisSpec ptAxisHigh_eff = {18, 10., 100., "#it{p}_{T} (GeV/#it{c})"};
-      AxisSpec etaAxis_eff = {100, -1.0, 1.0, "#eta"};
-      AxisSpec phiAxis_eff = {200, -1.0, 7., "#phi"};
+      AxisSpec etaAxis_eff = {etaEffNBins, -1.0, 1.0, "#eta"};
+      AxisSpec phiAxis_eff = {phiEffNBins, -1.0, 7., "#phi"};
 
       // ptAxisLow
       registry.add("h3_particle_pt_particle_eta_particle_phi_mcpartofinterest", "#it{p}_{T, mcpart} (GeV/#it{c}); #eta_{mcpart}; #phi_{mcpart}", {HistType::kTH3F, {ptAxis_eff, etaAxis_eff, phiAxis_eff}});
@@ -398,6 +400,168 @@ struct TrackEfficiencyJets {
   }
   PROCESS_SWITCH(TrackEfficiencyJets, processEFficiencyPurity, "Histograms for efficiency and purity quantities", true);
 
+  void processEFficiencyPurityWeighted(aod::JetMcCollision const& mcCollision,
+                                       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,
+                                       JetParticlesWithOriginal const& jMcParticles)
+  {
+    // missing:
+    //   * constexpr auto hasCentrality = CollisionMCRecTableCentFT0C::template contains<aod::CentFT0Cs>();
+    //           if constexpr (hasCentrality) {
+    //   * dividing in centrality bins
+    // I should maybe introduce the sel8 cuts on the collisoins (reco, but what about mccoll? maybe not htat way included in efficiency)
+
+    registry.fill(HIST("hMcCollCutsCounts"), 0.5); // all mcCollisions
+
+    if (!(abs(mcCollision.posZ()) < vertexZCut)) {
+      return;
+    }
+    registry.fill(HIST("hMcCollCutsCounts"), 1.5); // mcCollision.posZ() condition
+
+    if (collisions.size() < 1) {
+      return;
+    }
+    registry.fill(HIST("hMcCollCutsCounts"), 2.5); // mcCollisions with at least one reconstructed collision
+
+    if (acceptSplitCollisions == 0 && collisions.size() > 1) {
+      return;
+    }
+    registry.fill(HIST("hMcCollCutsCounts"), 3.5); // split mcCollisions condition
+
+    bool hasSel8Coll = false;
+    bool centralityCheck = false;
+    if (acceptSplitCollisions == 2) {                                                                          // check only that the first reconstructed collision passes the check
+      if (jetderiveddatautilities::selectCollision(collisions.begin(), eventSelectionBits, skipMBGapEvents)) { // Skipping MC events that have not a single selected reconstructed collision ; effect unclear if mcColl is split
+        hasSel8Coll = true;
+      }
+      if (!checkCentrality || ((centralityMin < collisions.begin().centrality()) && (collisions.begin().centrality() < centralityMax))) { // effect unclear if mcColl is split
+        centralityCheck = true;
+      }
+    } else { // check that at least one of the reconstructed collisions passes the checks
+      for (auto& collision : collisions) {
+        if (jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents)) { // Skipping MC events that have not a single selected reconstructed collision ; effect unclear if mcColl is split
+          hasSel8Coll = true;
+        }
+        if (!checkCentrality || ((centralityMin < collision.centrality()) && (collision.centrality() < centralityMax))) { // effect unclear if mcColl is split
+          centralityCheck = true;
+        }
+      }
+    }
+    if (!hasSel8Coll) {
+      return;
+    }
+    registry.fill(HIST("hMcCollCutsCounts"), 4.5); // at least one of the reconstructed collisions associated with this mcCollision is selected
+
+    if (!centralityCheck) {
+      return;
+    }
+    registry.fill(HIST("hMcCollCutsCounts"), 5.5); // at least one of the reconstructed collisions associated with this mcCollision is selected with regard to centrality
+
+    float eventWeight = mcCollision.weight();
+
+    for (auto& jMcParticle : jMcParticles) {
+      registry.fill(HIST("hMcPartCutsCounts"), 0.5); // allPartsInSelMcColl
+
+      if (!isChargedParticle(jMcParticle.pdgCode())) {
+        continue;
+      }
+      registry.fill(HIST("hMcPartCutsCounts"), 1.5); // isCharged
+
+      registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_mcpart_nonprimary"), jMcParticle.pt(), jMcParticle.eta(), jMcParticle.phi(), eventWeight);
+
+      if (checkPrimaryPart && !jMcParticle.isPhysicalPrimary()) { // global tracks should be mostly primaries
+        continue;
+      }
+      registry.fill(HIST("hMcPartCutsCounts"), 2.5); // isPrimary
+
+      registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_mcpartofinterest"), jMcParticle.pt(), jMcParticle.eta(), jMcParticle.phi(), eventWeight);
+
+      registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_mcpartofinterest"), jMcParticle.pt(), jMcParticle.eta(), jMcParticle.phi(), eventWeight);
+
+      if ((abs(jMcParticle.eta()) < trackEtaAcceptanceCountQA)) { // removed from actual cuts for now because all the histograms have an eta axis
+        registry.fill(HIST("hMcPartCutsCounts"), 3.5);            // etaAccept // not actually applied here but it will give an idea of what will be done in the post processing
+      }
+    }
+
+    std::vector<int> seenMcParticlesVector; // is reset every mc collision
+
+    int splitCollCounter = 0;
+    for (auto& collision : collisions) {
+      splitCollCounter++;
+      if (acceptSplitCollisions == 2 && splitCollCounter > 1) {
+        return;
+      }
+
+      if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents) || !(abs(collision.posZ()) < vertexZCut)) {
+        continue;
+      }
+
+      auto collTracks = jetTracks.sliceBy(tracksPerJCollision, collision.globalIndex());
+      for (auto& 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
+          continue;
+        }
+        registry.fill(HIST("hTrackCutsCounts"), 1.5);
+
+        if (!track.has_mcParticle()) {
+          registry.fill(HIST("h3_track_pt_track_eta_track_phi_nonassociatedtrack"), track.pt(), track.eta(), track.phi(), eventWeight);
+
+          registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_nonassociatedtrack"), track.pt(), track.eta(), track.phi(), eventWeight);
+          continue;
+        }
+        registry.fill(HIST("hTrackCutsCounts"), 2.5);
+
+        auto jMcParticleFromTrack = track.mcParticle_as<JetParticlesWithOriginal>();
+        if (!jMcParticleFromTrack.isPhysicalPrimary()) {
+          registry.fill(HIST("h3_track_pt_track_eta_track_phi_associatedtrack_nonprimary"), track.pt(), track.eta(), track.phi(), eventWeight);
+          registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_associatedtrack_nonprimary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi(), eventWeight);
+
+          registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_associatedtrack_nonprimary"), track.pt(), track.eta(), track.phi(), eventWeight);
+          registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_associatedtrack_nonprimary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi(), eventWeight);
+
+          if (std::find(seenMcParticlesVector.begin(), seenMcParticlesVector.end(), jMcParticleFromTrack.globalIndex()) != seenMcParticlesVector.end()) {
+            registry.fill(HIST("h3_track_pt_track_eta_track_phi_associatedtrack_split_nonprimary"), track.pt(), track.eta(), track.phi(), eventWeight);
+            registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_associatedtrack_split_nonprimary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi(), eventWeight);
+
+            registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_associatedtrack_split_nonprimary"), track.pt(), track.eta(), track.phi(), eventWeight);
+            registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_associatedtrack_split_nonprimary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi(), eventWeight);
+          } else {
+            seenMcParticlesVector.push_back(jMcParticleFromTrack.globalIndex());
+          }
+
+          continue;
+        }
+
+        registry.fill(HIST("hTrackCutsCounts"), 3.5);
+
+        registry.fill(HIST("h3_track_pt_track_eta_track_phi_associatedtrack_primary"), track.pt(), track.eta(), track.phi(), eventWeight);
+        registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_associatedtrack_primary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi(), eventWeight);
+        registry.fill(HIST("h2_particle_pt_track_pt_residual_associatedtrack_primary"), jMcParticleFromTrack.pt(), (jMcParticleFromTrack.pt() - track.pt()) / jMcParticleFromTrack.pt(), eventWeight);
+
+        registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_associatedtrack_primary"), track.pt(), track.eta(), track.phi(), eventWeight);
+        registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_associatedtrack_primary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi(), eventWeight);
+        registry.fill(HIST("h2_particle_pt_high_track_pt_high_residual_associatedtrack_primary"), jMcParticleFromTrack.pt(), (jMcParticleFromTrack.pt() - track.pt()) / jMcParticleFromTrack.pt(), eventWeight);
+
+        if (std::find(seenMcParticlesVector.begin(), seenMcParticlesVector.end(), jMcParticleFromTrack.globalIndex()) != seenMcParticlesVector.end()) {
+          registry.fill(HIST("h3_track_pt_track_eta_track_phi_associatedtrack_split_primary"), track.pt(), track.eta(), track.phi(), eventWeight);
+          registry.fill(HIST("h3_particle_pt_particle_eta_particle_phi_associatedtrack_split_primary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi(), eventWeight);
+
+          registry.fill(HIST("h3_track_pt_high_track_eta_track_phi_associatedtrack_split_primary"), track.pt(), track.eta(), track.phi(), eventWeight);
+          registry.fill(HIST("h3_particle_pt_high_particle_eta_particle_phi_associatedtrack_split_primary"), jMcParticleFromTrack.pt(), jMcParticleFromTrack.eta(), jMcParticleFromTrack.phi(), eventWeight);
+        } else {
+          seenMcParticlesVector.push_back(jMcParticleFromTrack.globalIndex());
+        }
+
+        if (abs(jMcParticleFromTrack.eta()) < trackEtaAcceptanceCountQA) { // not actually applied here but it will give an idea of what will be done in the post processing
+          registry.fill(HIST("hTrackCutsCounts"), 4.5);
+        }
+      }
+    }
+  }
+  PROCESS_SWITCH(TrackEfficiencyJets, processEFficiencyPurityWeighted, "Histograms for efficiency and purity quantities for weighted simulations", false);
+
   void processTracks(soa::Filtered<aod::JetCollisions>::iterator const& collision,
                      soa::Filtered<soa::Join<aod::JetTracks, aod::JTrackExtras>> const& tracks)
   {
@@ -499,7 +663,7 @@ struct TrackEfficiencyJets {
       return;
     }
 
-    float eventWeight = mcCollision.weight();
+    float eventWeight = eventWeight;
     registry.fill(HIST("h_mccollisions"), 0.5);
     registry.fill(HIST("h_mccollisions_weighted"), 0.5, eventWeight);