Select MC charged particles

Author: omvazque

PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx

@@ -69,6 +69,7 @@ struct UccZdc {
 
   static constexpr float kCollEnergy{2.68};
   static constexpr float kZero{0.};
+  static constexpr float minCharge{3.f};
 
   // Configurables Event Selection
   Configurable<bool> isNoCollInTimeRangeStrict{"isNoCollInTimeRangeStrict", true, "use isNoCollInTimeRangeStrict?"};
@@ -819,6 +820,7 @@ struct UccZdc {
 
   // Preslice<aod::McParticles> perMCCollision = aod::mcparticle::mcCollisionId;
   Preslice<TheFilteredSimTracks> perCollision = aod::track::collisionId;
+  Service<o2::framework::O2DatabasePDG> pdg;
   TRandom* randPointer = new TRandom();
   void processMCclosure(aod::McCollisions::iterator const& mccollision, soa::SmallGroups<o2::aod::SimCollisions> const& collisions, o2::aod::BCsRun3 const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::McParticles const& mcParticles, TheFilteredSimTracks const& simTracks)
   {
@@ -879,9 +881,10 @@ struct UccZdc {
           return;
         }
 
-        std::vector<float> pTs;
-        std::vector<float> vecFD;
-        std::vector<float> vecOneOverEff;
+        std::vector<double> pTs;
+        std::vector<double> vecFD;
+        std::vector<double> vecEff;
+        std::vector<double> vecOneOverEffXFD;
         // std::vector<float> wIs;
         const auto& groupedTracks{simTracks.sliceBy(perCollision, collision.globalIndex())};
 
@@ -901,6 +904,8 @@ struct UccZdc {
         }
 
         // Calculates the event weight, W_k
+        const int foundNchBin{efficiency->GetXaxis()->FindBin(nchRaw)};
+
         for (const auto& track : groupedTracks) {
           // Track Selection
           if (track.eta() < minEta || track.eta() > maxEta) {
@@ -912,32 +917,51 @@ struct UccZdc {
           if (!track.isGlobalTrack()) {
             continue;
           }
+          if (!track.has_mcParticle()) {
+            continue;
+          }
+          const auto& particle{track.mcParticle()};
+
+          auto charge{0.};
+          // Get the MC particle
+          auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
+          if (pdgParticle != nullptr) {
+            charge = pdgParticle->Charge();
+          } else {
+            continue;
+          }
 
-          float pt{track.pt()};
-          int foundNchBin{efficiency->GetXaxis()->FindBin(nchRaw)};
-          int foundPtBin{efficiency->GetYaxis()->FindBin(pt)};
-          float effValue{1.};
-          float fdValue{1.};
+          // Is it a charged particle?
+          if (std::abs(charge) < minCharge) {
+            continue;
+          }
+          // Is it a primary particle?
+          // if (!particle.isPhysicalPrimary()) { continue; }
+
+          const double pt{static_cast<double>(track.pt())};
+          const int foundPtBin{efficiency->GetYaxis()->FindBin(pt)};
+          double effValue{1.};
+          double fdValue{1.};
 
           if (applyEff) {
             effValue = efficiency->GetBinContent(foundNchBin, foundPtBin);
             fdValue = fd->GetBinContent(fd->FindBin(pt));
           }
           if ((effValue > 0.) && (fdValue > 0.)) {
             pTs.emplace_back(pt);
-            vecOneOverEff.emplace_back(1. / effValue);
+            vecEff.emplace_back(effValue);
             vecFD.emplace_back(fdValue);
+            vecOneOverEffXFD.emplace_back(fdValue / effValue);
           }
         }
-
-        nchMult = std::accumulate(vecOneOverEff.begin(), vecOneOverEff.end(), 0);
+        nchMult = std::accumulate(vecOneOverEffXFD.begin(), vecOneOverEffXFD.end(), 0);
         if (nchMult < minNchSel) {
           return;
         }
 
         double p1, p2, p3, p4, w1, w2, w3, w4;
         p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0;
-        getPTpowers(pTs, vecOneOverEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4);
+        getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4);
 
         const double denTwoParCorr{std::pow(w1, 2.) - w2};
         const double numTwoParCorr{std::pow(p1, 2.) - p2};
@@ -971,6 +995,21 @@ struct UccZdc {
           if (particle.pt() < minPt || particle.pt() > maxPt) {
             continue;
           }
+
+          auto charge{0.};
+          // Get the MC particle
+          auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
+          if (pdgParticle != nullptr) {
+            charge = pdgParticle->Charge();
+          } else {
+            continue;
+          }
+
+          // Is it a charged particle?
+          if (std::abs(charge) < minCharge) {
+            continue;
+          }
+          // Is it a primary particle?
           if (!particle.isPhysicalPrimary()) {
             continue;
           }
@@ -1044,9 +1083,23 @@ struct UccZdc {
           if (!track.has_mcParticle()) {
             continue;
           }
-          registry.fill(HIST("Pt_all_ch"), nchRaw, track.pt());
-
+          // Get the MC particle
           const auto& particle{track.mcParticle()};
+          auto charge{0.};
+          auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
+          if (pdgParticle != nullptr) {
+            charge = pdgParticle->Charge();
+          } else {
+            continue;
+          }
+
+          // Is it a charged particle?
+          if (std::abs(charge) < minCharge) {
+            continue;
+          }
+          // All charged particles
+          registry.fill(HIST("Pt_all_ch"), nchRaw, track.pt());
+          // Is it a primary particle?
           if (!particle.isPhysicalPrimary()) {
             continue;
           }
@@ -1075,6 +1128,21 @@ struct UccZdc {
           if (particle.pt() < minPt || particle.pt() > maxPt) {
             continue;
           }
+
+          auto charge{0.};
+          // Get the MC particle
+          auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
+          if (pdgParticle != nullptr) {
+            charge = pdgParticle->Charge();
+          } else {
+            continue;
+          }
+
+          // Is it a charged particle?
+          if (std::abs(charge) < minCharge) {
+            continue;
+          }
+          // Is it a primary particle?
           if (!particle.isPhysicalPrimary()) {
             continue;
           }
@@ -1101,14 +1169,14 @@ struct UccZdc {
   PROCESS_SWITCH(UccZdc, processMCclosure, "Process MC closure", false);
 
   template <typename T, typename U>
-  void getPTpowers(const T& pTs, const T& vecOneOverEff, const T& vecFD, U& pOne, U& wOne, U& pTwo, U& wTwo, U& pThree, U& wThree, U& pFour, U& wFour)
+  void getPTpowers(const T& pTs, const T& vecEff, const T& vecFD, U& pOne, U& wOne, U& pTwo, U& wTwo, U& pThree, U& wThree, U& pFour, U& wFour)
   {
     pOne = wOne = pTwo = wTwo = pThree = wThree = pFour = wFour = 0.;
     for (std::size_t i = 0; i < pTs.size(); ++i) {
-      const float pTi{pTs.at(i)};
-      const float eFFi{vecOneOverEff.at(i)};
-      const float fDi{vecFD.at(i)};
-      const float wEighti{eFFi * fDi};
+      const double pTi{pTs.at(i)};
+      const double eFFi{vecEff.at(i)};
+      const double fDi{vecFD.at(i)};
+      const double wEighti{pow(eFFi, -1.) * fDi};
       pOne += wEighti * pTi;
       wOne += wEighti;
       pTwo += std::pow(wEighti * pTi, 2.);