[PWGCF] FU: adding y cut for D0s and possibility to run with prompt D0s

PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerMCTruthTask.cxx

@@ -106,6 +106,9 @@ struct FemtoUniverseProducerMCTruthTask {
     Configurable<float> confEtaFilterCut{"confEtaFilterCut", 0.8, "Eta cut for the filtering tracks"};
   } ConfFilteringTracks;
 
+  // D0/D0bar cuts
+  Configurable<float> yD0CandGenMax{"yD0CandGenMax", 0.5, "Rapidity cut for the D0/D0bar mesons"};
+
   FemtoUniverseCollisionSelection colCuts;
   FemtoUniverseTrackSelection trackCuts;
   HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::QAObject};
@@ -162,8 +165,6 @@ struct FemtoUniverseProducerMCTruthTask {
       /// if the most open selection criteria are not fulfilled there is no
       /// point looking further at the track
 
-      if (particle.eta() < -ConfFilteringTracks.confEtaFilterCut || particle.eta() > ConfFilteringTracks.confEtaFilterCut)
-        continue;
       if (particle.pt() < ConfFilteringTracks.confPtLowFilterCut || particle.pt() > ConfFilteringTracks.confPtHighFilterCut)
         continue;
 
@@ -188,6 +189,24 @@ struct FemtoUniverseProducerMCTruthTask {
           continue;
       }
 
+      // check if D0/D0bar mesons pass the rapidity cut
+      // if pass then saving the orgin of D0/D0bar
+      // check if tracks (besides D0/D0bar) pass pseudorapidity cut
+      int8_t origin = -99;
+      if (std::abs(particle.pdgCode()) == Pdg::kD0) {
+        if (std::abs(particle.y()) > yD0CandGenMax) {
+          continue;
+        } else {
+          origin = RecoDecay::getCharmHadronOrigin(tracks, particle);
+        }
+      } else {
+        if (std::abs(particle.eta()) > ConfFilteringTracks.confEtaFilterCut) {
+          continue;
+        } else {
+          origin = -99;
+        }
+      }
+
       /// check if we end-up with the correct final state using MC info
       int8_t sign = 0;
       if (std::abs(pdgCode) == Pdg::kD0 && !RecoDecay::isMatchedMCGen(tracks, particle, Pdg::kD0, std::array{+kPiPlus, -kKPlus}, true, &sign)) {
@@ -216,8 +235,8 @@ struct FemtoUniverseProducerMCTruthTask {
                   pdgCode,
                   pdgCode,
                   childIDs,
-                  0,
-                  0);
+                  origin,
+                  -999.);
     }
   }
 

PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx

@@ -55,7 +55,7 @@
 #include <CCDB/BasicCCDBManager.h>
 
 #include "Math/Vector4D.h"
 #include "TLorentzVector.h"
 #include "TMath.h"
 #include <TPDGCode.h>
 
@@ -323,9 +323,12 @@
 
   // D0/D0bar mesons
   struct : o2::framework::ConfigurableGroup {
-    Configurable<float> confD0D0barCandEtaCut{"confD0D0barCandEtaCut", 0.8, "max. cand. pseudorapidity"};
+    Configurable<float> trackD0CandEtaMax{"trackD0CandEtaMax", 0.8, "max. track/D0 cand. pseudorapidity"};
+    Configurable<double> yD0CandGenMax{"yD0CandGenMax", 0.5, "max. gen. D0 cand. rapidity"};
+    Configurable<double> yD0CandMax{"yD0CandMax", 0.8, "max. D0 cand. rapidity"};
     Configurable<float> trackD0pTGenMin{"trackD0pTGenMin", 0.0, "MC Truth, min. pT for tracks and D0/D0bar cand."};
     Configurable<float> trackD0pTGenMax{"trackD0pTGenMax", 24.0, "MC Truth, max. pT for tracks and D0/D0bar cand."};
+    Configurable<bool> useYCutD0Cand{"useYCutD0Cand", true, "True - apply cut on y of D0 cand./false - apply cut on eta"};
     Configurable<bool> storeD0D0barDoubleMassHypo{"storeD0D0barDoubleMassHypo", false, "Store D0/D0bar cand. which pass selection criteria for both, D0 and D0bar"};
     Configurable<std::vector<int>> classMlD0D0bar{"classMlD0D0bar", {0, 1, 2}, "Indexes of ML scores to be stored. Three indexes max."};
   } ConfD0Selection;
@@ -450,11 +453,11 @@
   }
 
   template <typename TrackType>
   aod::femtouniverseparticle::CutContainerType PIDBitmask(const TrackType& track)
   {
     static const o2::track::PID pids[] = {o2::track::PID::Proton, o2::track::PID::Pion, o2::track::PID::Kaon};
     aod::femtouniverseparticle::CutContainerType mask = 0u;
     for (UInt_t i = 0; i < 3; ++i) {
       if (isNSigmaTPC(trackCuts.getNsigmaTPC(track, pids[i])))
         mask |= (1u << i);
       if (isNSigmaTOF(track.p(), trackCuts.getNsigmaTOF(track, pids[i]), track.hasTOF()))
@@ -938,8 +941,8 @@
         phiOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kFake;
       }
 
       TLorentzVector part1Vec;
       TLorentzVector part2Vec;
 
       const auto mMassOne = o2::constants::physics::MassKPlus;  // FIXME: Get from the PDG service of the common header
       const auto mMassTwo = o2::constants::physics::MassKMinus; // FIXME: Get from the PDG service of the common header
@@ -947,7 +950,7 @@
       part1Vec.SetPtEtaPhiM(kaon1MC.pt(), kaon1MC.eta(), kaon1MC.phi(), mMassOne);
       part2Vec.SetPtEtaPhiM(kaon2MC.pt(), kaon2MC.eta(), kaon2MC.phi(), mMassTwo);
 
       TLorentzVector sumVec(part1Vec);
       sumVec += part2Vec;
 
       float phiEta = sumVec.Eta();
@@ -1484,7 +1487,11 @@
         continue;
       }
 
-      if (std::abs(hfCand.eta()) > ConfD0Selection.confD0D0barCandEtaCut) {
+      if (ConfD0Selection.useYCutD0Cand && std::abs(hfHelper.yD0(hfCand)) > ConfD0Selection.yD0CandMax) {
+        continue;
+      }
+
+      if (!(ConfD0Selection.useYCutD0Cand) && std::abs(hfCand.eta()) > ConfD0Selection.trackD0CandEtaMax) {
         continue;
       }
 
@@ -1603,7 +1610,11 @@
         continue;
       }
 
-      if (std::abs(hfCand.eta()) > ConfD0Selection.confD0D0barCandEtaCut) {
+      if (ConfD0Selection.useYCutD0Cand && std::abs(hfHelper.yD0(hfCand)) > ConfD0Selection.yD0CandMax) {
+        continue;
+      }
+
+      if (!(ConfD0Selection.useYCutD0Cand) && std::abs(hfCand.eta()) > ConfD0Selection.trackD0CandEtaMax) {
         continue;
       }
 
@@ -1729,7 +1740,11 @@
         continue;
       }
 
-      if (std::abs(hfCand.eta()) > ConfD0Selection.confD0D0barCandEtaCut) {
+      if (ConfD0Selection.useYCutD0Cand && std::abs(hfHelper.yD0(hfCand)) > ConfD0Selection.yD0CandMax) {
+        continue;
+      }
+
+      if (!(ConfD0Selection.useYCutD0Cand) && std::abs(hfCand.eta()) > ConfD0Selection.trackD0CandEtaMax) {
         continue;
       }
 
@@ -1889,8 +1904,8 @@
         continue;
       }
 
       TLorentzVector part1Vec;
       TLorentzVector part2Vec;
 
       const auto mMassOne = o2::constants::physics::MassKPlus;  // FIXME: Get from the PDG service of the common header
       const auto mMassTwo = o2::constants::physics::MassKMinus; // FIXME: Get from the PDG service of the common header
@@ -1898,7 +1913,7 @@
       part1Vec.SetPtEtaPhiM(p1.pt(), p1.eta(), p1.phi(), mMassOne);
       part2Vec.SetPtEtaPhiM(p2.pt(), p2.eta(), p2.phi(), mMassTwo);
 
       TLorentzVector sumVec(part1Vec);
       sumVec += part2Vec;
 
       float phiEta = sumVec.Eta();
@@ -2118,8 +2133,7 @@
     std::vector<int> tmpIDtrack;
 
     for (const auto& particle : mcParts) {
-      if (particle.eta() < -ConfFilterCuts.confEtaFilterCut || particle.eta() > ConfFilterCuts.confEtaFilterCut)
-        continue;
+
       if (particle.pt() < ConfD0Selection.trackD0pTGenMin || particle.pt() > ConfD0Selection.trackD0pTGenMax)
         continue;
 
@@ -2152,9 +2166,18 @@
         /// check if we have D0(bar) → π± K∓
         continue;
       }
+
       if (std::abs(particle.pdgCode()) == Pdg::kD0) {
-        origin = RecoDecay::getCharmHadronOrigin(mcParts, particle);
-        mcGenFlag = particle.flagMcMatchGen();
+        if (std::abs(particle.y()) > ConfD0Selection.yD0CandGenMax) {
+          continue;
+        } else {
+          origin = RecoDecay::getCharmHadronOrigin(mcParts, particle);
+          mcGenFlag = particle.flagMcMatchGen();
+        }
+      } else {
+        if (std::abs(particle.eta()) > ConfD0Selection.trackD0CandEtaMax) {
+          continue;
+        }
       }
 
       outputParts(outputCollision.lastIndex(),