Addition of MC truth process

Author: Preet Pati

PWGCF/TwoParticleCorrelations/Tasks/pidDiHadron.cxx

@@ -183,9 +183,11 @@ struct PidDiHadron {
 
   // make the filters and cuts.
   Filter collisionFilter = (nabs(aod::collision::posZ) < cfgCutVertex);
-  Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtPOIMin) && (aod::track::pt < cfgCutPtPOIMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls) && (nabs(aod::track::dcaZ) < cfgCutDCAz);
+  Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtPOIMin) && (aod::track::pt < cfgCutPtPOIMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls) && (nabs(aod::track::dcaZ) < cfgCutDCAz);
+
   using FilteredCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSel, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::Mults>>;
   using FilteredTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTOFbeta, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr>>;
+  using FilteredMcTracks = soa::Filtered<soa::Join<aod::Tracks, aod::McTrackLabels, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTOFbeta, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr>>;
   using V0TrackCandidate = aod::V0Datas;
 
   Preslice<aod::Tracks> perCollision = aod::track::collisionId;
@@ -356,7 +358,7 @@ struct PidDiHadron {
       massAxisReso = {resoSwitchVals[kMassBins][iPhi], resoCutVals[kMassMin][iPhi], resoCutVals[kMassMax][iPhi], "M_{K^{+}K^{-}} (GeV/c^{2})"};
 
     // Event Counter
-    if ((doprocessSame || doprocessSameReso) && cfgUseAdditionalEventCut) {
+    if ((doprocessSame || doprocessSameReso || doprocessMC) && cfgUseAdditionalEventCut) {
       histos.add("hEventCount", "Number of Events;; Count", {HistType::kTH1D, {{15, -0.5, 14.5}}});
       histos.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(kFilteredEvents + 1, "Filtered event");
       histos.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(kAfterSel8 + 1, "After sel8");
@@ -512,6 +514,15 @@ struct PidDiHadron {
     if (doprocessMixed || doprocessMixedReso) {
       histos.add("deltaEta_deltaPhi_mixed", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}});
     }
+    if (doprocessMC) {
+      histos.add("hNsigmaPionSelected", "hNsigmaPionSelected", {HistType::kTH1D, {axisPt}});
+      histos.add("hNsigmaKaonSelected", "hNsigmaKaonSelected", {HistType::kTH1D, {axisPt}});
+      histos.add("hNsigmaProtonSelected", "hNsigmaProtonSelected", {HistType::kTH1D, {axisPt}});
+
+      histos.add("hNsigmaPionTrue", "hNsigmaPionTrue", {HistType::kTH1D, {axisPt}});
+      histos.add("hNsigmaKaonTrue", "hNsigmaKaonTrue", {HistType::kTH1D, {axisPt}});
+      histos.add("hNsigmaProtonTrue", "hNsigmaProtonTrue", {HistType::kTH1D, {axisPt}});
+    }
 
     histos.add("eventcount", "bin", {HistType::kTH1F, {{4, 0, 4, "bin"}}}); // histogram to see how many events are in the same and mixed event
 
@@ -1467,6 +1478,58 @@ struct PidDiHadron {
     }
   }
   PROCESS_SWITCH(PidDiHadron, processMixedReso, "Process mixed events", true);
+
+  void processMC(FilteredCollisions::iterator const& collision, FilteredMcTracks const& tracksmc, aod::BCsWithTimestamps const&, aod::McParticles const&)
+  {
+    float cent = -1.;
+    if (!cfgCentTableUnavailable)
+      cent = getCentrality(collision);
+    if (cfgUseAdditionalEventCut && !selectionEvent(collision, tracksmc.size(), cent, true))
+      return;
+
+    if (cfgSelCollByNch && (tracksmc.size() < cfgCutMultMin || tracksmc.size() >= cfgCutMultMax)) {
+      return;
+    }
+    if (!cfgSelCollByNch && !cfgCentTableUnavailable && (cent < cfgCutCentMin || cent >= cfgCutCentMax)) {
+      return;
+    }
+
+    // loop over all tracks
+    for (auto const& track : tracksmc) {
+
+      if (!trackSelected(track))
+        continue;
+
+      int pidIndex = getNsigmaPID(track);
+
+      // Check the PDG code for the particles (MC truth) and match with analysed Nsigma PID
+      if (std::abs(track.mcParticle().pdgCode()) == PDG_t::kPiPlus) {
+        histos.fill(HIST("hNsigmaPionTrue"), track.pt());
+
+        if (pidIndex == kPions) {
+          histos.fill(HIST("hNsigmaPionSelected"), track.pt());
+        }
+      } // Pion condition
+
+      if (std::abs(track.mcParticle().pdgCode()) == PDG_t::kKPlus) {
+        histos.fill(HIST("hNsigmaKaonTrue"), track.pt());
+
+        if (pidIndex == kKaons) {
+          histos.fill(HIST("hNsigmaKaonSelected"), track.pt());
+        }
+      } // Kaon condition
+
+      if (std::abs(track.mcParticle().pdgCode()) == PDG_t::kProton) {
+        histos.fill(HIST("hNsigmaProtonTrue"), track.pt());
+
+        if (pidIndex == kProtons) {
+          histos.fill(HIST("hNsigmaProtonSelected"), track.pt());
+        }
+      } // Proton condition
+
+    } // end of tracks MC loop
+  } // end of process MC
+  PROCESS_SWITCH(PidDiHadron, processMC, "Process Monte Carlo", true);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)