[PWGCF,PWGHF] Add ITS PID to femtodream (#9317)

Author: ariedel-cern

PWGCF/DataModel/FemtoDerived.h

@@ -175,6 +175,14 @@ DECLARE_SOA_COLUMN(TOFNSigmaPr, tofNSigmaPr, float); //! Nsigma separation with
 DECLARE_SOA_COLUMN(TOFNSigmaDe, tofNSigmaDe, float); //! Nsigma separation with the TOF detector for deuteron
 DECLARE_SOA_COLUMN(TOFNSigmaTr, tofNSigmaTr, float); //! Nsigma separation with the TOF detector for triton
 DECLARE_SOA_COLUMN(TOFNSigmaHe, tofNSigmaHe, float); //! Nsigma separation with the TOF detector for helium3
+DECLARE_SOA_COLUMN(ITSSignal, itsSignal, float);
+DECLARE_SOA_COLUMN(ITSNSigmaEl, itsNSigmaEl, float); //! Nsigma separation with the Its detector for electron
+DECLARE_SOA_COLUMN(ITSNSigmaPi, itsNSigmaPi, float); //! Nsigma separation with the Its detector for pion
+DECLARE_SOA_COLUMN(ITSNSigmaKa, itsNSigmaKa, float); //! Nsigma separation with the Its detector for kaon
+DECLARE_SOA_COLUMN(ITSNSigmaPr, itsNSigmaPr, float); //! Nsigma separation with the Its detector for proton
+DECLARE_SOA_COLUMN(ITSNSigmaDe, itsNSigmaDe, float); //! Nsigma separation with the Its detector for deuteron
+DECLARE_SOA_COLUMN(ITSNSigmaTr, itsNSigmaTr, float); //! Nsigma separation with the Its detector for triton
+DECLARE_SOA_COLUMN(ITSNSigmaHe, itsNSigmaHe, float); //! Nsigma separation with the Its detector for helium3
 DECLARE_SOA_COLUMN(DaughDCA, daughDCA, float);       //! DCA between daughters
 DECLARE_SOA_COLUMN(TransRadius, transRadius, float); //! Transverse radius of the decay vertex
 DECLARE_SOA_COLUMN(DecayVtxX, decayVtxX, float);     //! X position of the decay vertex
@@ -359,6 +367,14 @@ DECLARE_SOA_TABLE_STAGED(FDExtParticles, "FDEXTPARTICLE",
                          femtodreamparticle::TOFNSigmaDe,
                          femtodreamparticle::TOFNSigmaTr,
                          femtodreamparticle::TOFNSigmaHe,
+                         femtodreamparticle::ITSSignal,
+                         femtodreamparticle::ITSNSigmaEl,
+                         femtodreamparticle::ITSNSigmaPi,
+                         femtodreamparticle::ITSNSigmaKa,
+                         femtodreamparticle::ITSNSigmaPr,
+                         femtodreamparticle::ITSNSigmaDe,
+                         femtodreamparticle::ITSNSigmaTr,
+                         femtodreamparticle::ITSNSigmaHe,
                          femtodreamparticle::DaughDCA,
                          femtodreamparticle::TransRadius,
                          femtodreamparticle::DecayVtxX,
@@ -457,3 +473,4 @@ using MixingHash = MixingHashes::iterator;
 } // namespace o2::aod
 
 #endif // PWGCF_DATAMODEL_FEMTODERIVED_H_
+       //

PWGCF/FemtoDream/Core/femtoDreamParticleHisto.h

@@ -72,7 +72,7 @@ class FemtoDreamParticleHisto
 
   // comment
   template <o2::aod::femtodreamMCparticle::MCType mc, typename T>
-  void init_debug(std::string folderName, T& multAxis, T& multPercentileAxis, T& pTAxis, T& etaAxis, T& phiAxis, T& tempFitVarAxis, T& dcazAxis, T& NsigmaTPCAxis, T& NsigmaTOFAxis, T& NsigmaTPCTOFAxis, T& /*TPCclustersAxis*/, bool correlatedPlots)
+  void init_debug(std::string folderName, T& multAxis, T& multPercentileAxis, T& pTAxis, T& etaAxis, T& phiAxis, T& tempFitVarAxis, T& dcazAxis, T& NsigmaTPCAxis, T& NsigmaTOFAxis, T& NsigmaTPCTOFAxis, T& NsigmaITSAxis, bool correlatedPlots)
   {
 
     std::string folderSuffix = static_cast<std::string>(o2::aod::femtodreamMCparticle::MCTypeName[mc]).c_str();
@@ -117,6 +117,14 @@ class FemtoDreamParticleHisto
       mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaComb_d").c_str(), "n#sigma_{comb}^{d}", kTH2F, {pTAxis, NsigmaTPCTOFAxis});
       mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaComb_tr").c_str(), "n#sigma_{comb}^{tr}", kTH2F, {pTAxis, NsigmaTPCTOFAxis});
       mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaComb_he3").c_str(), "n#sigma_{comb}^{he3}", kTH2F, {pTAxis, NsigmaTPCTOFAxis});
+      mHistogramRegistry->add((folderName + folderSuffix + "/ITSSignal").c_str(), "<cluster size>x<cos#lambda>", kTH2F, {pTAxis, NsigmaITSAxis});
+      mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_el").c_str(), "n#sigma_{ITS}^{e}", kTH2F, {pTAxis, NsigmaITSAxis});
+      mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_pi").c_str(), "n#sigma_{ITS}^{#pi}", kTH2F, {pTAxis, NsigmaITSAxis});
+      mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_K").c_str(), "n#sigma_{ITS}^{K}", kTH2F, {pTAxis, NsigmaITSAxis});
+      mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_p").c_str(), "n#sigma_{ITS}^{p}", kTH2F, {pTAxis, NsigmaITSAxis});
+      mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_d").c_str(), "n#sigma_{ITS}^{d}", kTH2F, {pTAxis, NsigmaITSAxis});
+      mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_tr").c_str(), "n#sigma_{ITS}^{tr}", kTH2F, {pTAxis, NsigmaITSAxis});
+      mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_he3").c_str(), "n#sigma_{ITS}^{he3}", kTH2F, {pTAxis, NsigmaITSAxis});
       if (correlatedPlots) {
         mHistogramRegistry->add((folderName + folderSuffix + "/HighDcorrelator").c_str(), "", kTHnSparseF, {multAxis, multPercentileAxis, pTAxis, etaAxis, phiAxis, tempFitVarAxis, dcazAxis, NsigmaTPCAxis, NsigmaTOFAxis});
       }
@@ -217,7 +225,7 @@ class FemtoDreamParticleHisto
   /// \param tempFitVarBins binning of the tempFitVar (DCA_xy in case of tracks, CPA in case of V0s, etc.)
   /// \param isMC add Monte Carlo truth histograms to the output file
   template <typename T>
-  void init(HistogramRegistry* registry, T& MultBins, T& PercentileBins, T& pTBins, T& etaBins, T& phiBins, T& tempFitVarBins, T& NsigmaTPCBins, T& NsigmaTOFBins, T& NsigmaTPCTOFBins, T& TPCclustersBins, T& InvMassBins, bool isMC, int pdgCode, bool isDebug = false, bool correlatedPlots = false)
+  void init(HistogramRegistry* registry, T& MultBins, T& PercentileBins, T& pTBins, T& etaBins, T& phiBins, T& tempFitVarBins, T& NsigmaTPCBins, T& NsigmaTOFBins, T& NsigmaTPCTOFBins, T& NsigmaITSBins, T& InvMassBins, bool isMC, int pdgCode, bool isDebug = false, bool correlatedPlots = false)
   {
     mPDG = pdgCode;
     if (registry) {
@@ -247,15 +255,15 @@ class FemtoDreamParticleHisto
       framework::AxisSpec NsigmaTPCAxis = {NsigmaTPCBins, "n#sigma_{TPC}"};
       framework::AxisSpec NsigmaTOFAxis = {NsigmaTOFBins, "n#sigma_{TOF}"};
       framework::AxisSpec NsigmaTPCTOFAxis = {NsigmaTPCTOFBins, "n#sigma_{TPC+TOF}"};
-      framework::AxisSpec TPCclustersAxis = {TPCclustersBins, "TPC found clusters"};
+      framework::AxisSpec NsigmaITSAxis = {NsigmaITSBins, "n#sigma_{ITS}"};
       framework::AxisSpec InvMassAxis = {InvMassBins, "M_{inv} (GeV/#it{c}^{2})"};
 
       std::string folderName = (static_cast<std::string>(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]).c_str() + static_cast<std::string>(mFolderSuffix[mFolderSuffixType])).c_str();
 
       // Fill here the actual histogramms by calling init_base and init_MC
       init_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(folderName, tempFitVarAxisTitle, pTAxis, tempFitVarAxis, InvMassAxis, multAxis);
       if (isDebug) {
-        init_debug<o2::aod::femtodreamMCparticle::MCType::kRecon>(folderName, multAxis, multPercentileAxis, pTAxis, etaAxis, phiAxis, tempFitVarAxis, dcazAxis, NsigmaTPCAxis, NsigmaTOFAxis, NsigmaTPCTOFAxis, TPCclustersAxis, correlatedPlots);
+        init_debug<o2::aod::femtodreamMCparticle::MCType::kRecon>(folderName, multAxis, multPercentileAxis, pTAxis, etaAxis, phiAxis, tempFitVarAxis, dcazAxis, NsigmaTPCAxis, NsigmaTOFAxis, NsigmaTPCTOFAxis, NsigmaITSAxis, correlatedPlots);
       }
       if (isMC) {
         init_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(folderName, tempFitVarAxisTitle, pTAxis, tempFitVarAxis, InvMassAxis, multAxis);
@@ -350,6 +358,14 @@ class FemtoDreamParticleHisto
       mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaComb_d"), momentum, std::sqrt(part.tpcNSigmaDe() * part.tpcNSigmaDe() + part.tofNSigmaDe() * part.tofNSigmaDe()));
       mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaComb_tr"), momentum, std::sqrt(part.tpcNSigmaTr() * part.tpcNSigmaTr() + part.tofNSigmaTr() * part.tofNSigmaTr()));
       mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaComb_he3"), momentum, std::sqrt(part.tpcNSigmaHe() * part.tpcNSigmaHe() + part.tofNSigmaHe() * part.tofNSigmaHe()));
+      mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/ITSSignal"), momentum, part.itsSignal());
+      mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaITS_el"), momentum, part.itsNSigmaEl());
+      mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaITS_pi"), momentum, part.itsNSigmaPi());
+      mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaITS_K"), momentum, part.itsNSigmaKa());
+      mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaITS_p"), momentum, part.itsNSigmaPr());
+      mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaITS_d"), momentum, part.itsNSigmaDe());
+      mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaITS_tr"), momentum, part.itsNSigmaTr());
+      mHistogramRegistry->fill(HIST(o2::aod::femtodreamparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("/nSigmaITS_he3"), momentum, part.itsNSigmaHe());
 
       if (correlatedPlots) {
 

PWGCF/FemtoDream/Core/femtoDreamTrackSelection.h

@@ -24,6 +24,8 @@
 
 #include "PWGCF/DataModel/FemtoDerived.h"
 #include "Common/DataModel/TrackSelectionTables.h"
+#include "Common/DataModel/PIDResponse.h"
+#include "Common/DataModel/PIDResponseITS.h"
 #include "Common/Core/TrackSelection.h"
 #include "Common/Core/TrackSelectionDefaults.h"
 #include "PWGCF/FemtoDream/Core/femtoDreamObjectSelection.h"
@@ -129,6 +131,14 @@ class FemtoDreamTrackSelection : public FemtoDreamObjectSelection<float, femtoDr
   template <typename T>
   auto getNsigmaTOF(T const& track, o2::track::PID pid);
 
+  /// Computes the n_sigma for a track and a particle-type hypothesis in the ITS
+  /// \tparam T Data type of the track
+  /// \param track Track for which PID is evaluated
+  /// \param pid Particle species for which PID is evaluated
+  /// \return Value of n_{sigma, ITS}
+  template <typename T>
+  auto getNsigmaITS(T const& track, o2::track::PID pid);
+
   /// Checks whether the most open combination of all selection criteria is fulfilled
   /// \tparam T Data type of the track
   /// \param track Track
@@ -146,7 +156,7 @@ class FemtoDreamTrackSelection : public FemtoDreamObjectSelection<float, femtoDr
   /// \param Eta eta of the track
   /// \param Dca dca of the track with respect to primary vertex
   /// \return The bit-wise container for the selections, separately with all selection criteria, and the PID
-  template <typename cutContainerType, typename T, typename R>
+  template <bool useItsPid = false, typename cutContainerType, typename T, typename R>
   std::array<cutContainerType, 2> getCutContainer(T const& track, R Pt, R Eta, R Dcaxy);
 
   /// Some basic QA histograms
@@ -291,7 +301,7 @@ class FemtoDreamTrackSelection : public FemtoDreamObjectSelection<float, femtoDr
                                                                           "Maximal DCA_z (cm)",
                                                                           "Minimal DCA (cm)",
                                                                           "Maximal PID (nSigma)"}; ///< Helper information for the different selections
-};                                                                                                 // namespace femtoDream
+}; // namespace femtoDream
 
 template <o2::aod::femtodreamparticle::ParticleType part, o2::aod::femtodreamparticle::TrackType tracktype, typename cutContainerType>
 void FemtoDreamTrackSelection::init(HistogramRegistry* QAregistry, HistogramRegistry* Registry)
@@ -384,6 +394,28 @@ auto FemtoDreamTrackSelection::getNsigmaTOF(T const& track, o2::track::PID pid)
   return o2::aod::pidutils::tofNSigma(pid, track);
 }
 
+template <typename T>
+auto FemtoDreamTrackSelection::getNsigmaITS(T const& track, o2::track::PID pid)
+{
+  if (pid == o2::track::PID::Electron) {
+    return track.itsNSigmaEl();
+  } else if (pid == o2::track::PID::Pion) {
+    return track.itsNSigmaPi();
+  } else if (pid == o2::track::PID::Kaon) {
+    return track.itsNSigmaKa();
+  } else if (pid == o2::track::PID::Proton) {
+    return track.itsNSigmaPr();
+  } else if (pid == o2::track::PID::Deuteron) {
+    return track.itsNSigmaDe();
+  } else if (pid == o2::track::PID::Triton) {
+    return track.itsNSigmaTr();
+  } else if (pid == o2::track::PID::Helium3) {
+    return track.itsNSigmaHe();
+  }
+  // if nothing matched, return default value
+  return -999.f;
+}
+
 template <typename T>
 bool FemtoDreamTrackSelection::isSelectedMinimal(T const& track)
 {
@@ -460,7 +492,7 @@ bool FemtoDreamTrackSelection::isSelectedMinimal(T const& track)
   return true;
 }
 
-template <typename cutContainerType, typename T, typename R>
+template <bool useItsPid, typename cutContainerType, typename T, typename R>
 std::array<cutContainerType, 2> FemtoDreamTrackSelection::getCutContainer(T const& track, R Pt, R Eta, R Dca)
 {
   cutContainerType output = 0;
@@ -479,23 +511,30 @@ std::array<cutContainerType, 2> FemtoDreamTrackSelection::getCutContainer(T cons
   const auto dcaZ = track.dcaZ();
   const auto dca = Dca;
 
-  std::vector<float> pidTPC, pidTOF;
+  std::vector<float> pidTPC, pidTOF, pidITS;
   for (auto it : mPIDspecies) {
     pidTPC.push_back(getNsigmaTPC(track, it));
     pidTOF.push_back(getNsigmaTOF(track, it));
+    if constexpr (useItsPid) {
+      pidITS.push_back(getNsigmaITS(track, it));
+    }
   }
 
   float observable = 0.;
   for (auto& sel : mSelections) {
     const auto selVariable = sel.getSelectionVariable();
     if (selVariable == femtoDreamTrackSelection::kPIDnSigmaMax) {
-      /// PID needs to be handled a bit differently since we may need more than one species
+      /// PID needsgetNsigmaITSto be handled a bit differently since we may need more than one species
       for (size_t i = 0; i < mPIDspecies.size(); ++i) {
         auto pidTPCVal = pidTPC.at(i) - nSigmaPIDOffsetTPC;
         auto pidTOFVal = pidTOF.at(i) - nSigmaPIDOffsetTOF;
         auto pidComb = std::sqrt(pidTPCVal * pidTPCVal + pidTOFVal * pidTOFVal);
         sel.checkSelectionSetBitPID(pidTPCVal, outputPID);
         sel.checkSelectionSetBitPID(pidComb, outputPID);
+        if constexpr (useItsPid) {
+          auto pidITSVal = pidITS.at(i);
+          sel.checkSelectionSetBitPID(pidITSVal, outputPID);
+        }
       }
     } else {
       /// for the rest it's all the same

PWGCF/FemtoDream/Core/femtoDreamUtils.h

@@ -118,5 +118,22 @@ inline bool containsNameValuePair(const std::vector<T>& myVector, const std::str
   return false; // No match found
 }
 
+template <typename T>
+float itsSignal(T const& track)
+{
+  uint32_t clsizeflag = track.itsClusterSizes();
+  auto clSizeLayer0 = (clsizeflag >> (0 * 4)) & 0xf;
+  auto clSizeLayer1 = (clsizeflag >> (1 * 4)) & 0xf;
+  auto clSizeLayer2 = (clsizeflag >> (2 * 4)) & 0xf;
+  auto clSizeLayer3 = (clsizeflag >> (3 * 4)) & 0xf;
+  auto clSizeLayer4 = (clsizeflag >> (4 * 4)) & 0xf;
+  auto clSizeLayer5 = (clsizeflag >> (5 * 4)) & 0xf;
+  auto clSizeLayer6 = (clsizeflag >> (6 * 4)) & 0xf;
+  int numLayers = 7;
+  int sumClusterSizes = clSizeLayer1 + clSizeLayer2 + clSizeLayer3 + clSizeLayer4 + clSizeLayer5 + clSizeLayer6 + clSizeLayer0;
+  float cosLamnda = 1. / std::cosh(track.eta());
+  return (static_cast<float>(sumClusterSizes) / numLayers) * cosLamnda;
+};
+
 } // namespace o2::analysis::femtoDream
 #endif // PWGCF_FEMTODREAM_CORE_FEMTODREAMUTILS_H_

PWGCF/FemtoDream/Core/femtoDreamV0Selection.h

@@ -550,8 +550,8 @@ template <typename cutContainerType, typename C, typename V, typename T>
 std::array<cutContainerType, 5>
   FemtoDreamV0Selection::getCutContainer(C const& /*col*/, V const& v0, T const& posTrack, T const& negTrack)
 {
-  auto outputPosTrack = PosDaughTrack.getCutContainer<cutContainerType>(posTrack, v0.positivept(), v0.positiveeta(), v0.dcapostopv());
-  auto outputNegTrack = NegDaughTrack.getCutContainer<cutContainerType>(negTrack, v0.negativept(), v0.negativeeta(), v0.dcanegtopv());
+  auto outputPosTrack = PosDaughTrack.getCutContainer<false, cutContainerType>(posTrack, v0.positivept(), v0.positiveeta(), v0.dcapostopv());
+  auto outputNegTrack = NegDaughTrack.getCutContainer<false, cutContainerType>(negTrack, v0.negativept(), v0.negativeeta(), v0.dcanegtopv());
   cutContainerType output = 0;
   size_t counter = 0;