[PWGEM/PhotonMeson] Changes to BC wise pi0 task

PWGEM/PhotonMeson/DataModel/bcWiseTables.h

@@ -1,4 +1,4 @@
 // Copyright 2019-2020 CERN and copyright holders of ALICE O2.
 // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
 // All rights not expressly granted are reserved.
 //
@@ -19,6 +19,8 @@
 #ifndef PWGEM_PHOTONMESON_DATAMODEL_BCWISETABLES_H_
 #define PWGEM_PHOTONMESON_DATAMODEL_BCWISETABLES_H_
 
+#include <limits>
+
 #include "Framework/AnalysisDataModel.h"
 
 namespace o2::aod
@@ -65,13 +67,13 @@
 DECLARE_SOA_COLUMN(HasNoTFROFBorder, hasNoTFROFBorder, bool);           //! not in the TF border or ITS ROF border region
 DECLARE_SOA_COLUMN(StoredFT0MAmplitude, storedFT0MAmplitude, uint16_t); //! ft0a+c amplitude
 
 DECLARE_SOA_DYNAMIC_COLUMN(FT0MAmplitude, ft0Amplitude, [](uint16_t storedFT0MAmplitude) -> float { return storedFT0MAmplitude / emdownscaling::downscalingFactors[emdownscaling::kFT0Amp]; }); //! FT0M amplitude
 } // namespace bcwisebc
 DECLARE_SOA_TABLE(BCWiseBCs, "AOD", "BCWISEBC", //! table of bc wise centrality estimation and event selection input
                   o2::soa::Index<>, bcwisebc::HasFT0, bcwisebc::HasTVX, bcwisebc::HaskTVXinEMC, bcwisebc::HasEMCCell, bcwisebc::HasNoTFROFBorder,
                   bcwisebc::StoredFT0MAmplitude, bcwisebc::FT0MAmplitude<bcwisebc::StoredFT0MAmplitude>);
 
 DECLARE_SOA_INDEX_COLUMN(BCWiseBC, bcWiseBC); //! bunch crossing ID used as index
 
 namespace bcwisecollision
 {
@@ -96,14 +98,14 @@
 DECLARE_SOA_COLUMN(StoredTime, storedTime, int16_t);            //! cluster time (10 ps resolution)
 DECLARE_SOA_COLUMN(StoredIsExotic, storedIsExotic, bool);       //! flag to mark cluster as exotic
 
-DECLARE_SOA_DYNAMIC_COLUMN(Definition, definition, [](int8_t storedDefinition) -> int8_t { return storedDefinition; });                                   //! cluster definition, see EMCALClusterDefinition.h
-DECLARE_SOA_DYNAMIC_COLUMN(E, e, [](int16_t storedE) -> float { return storedE / emdownscaling::downscalingFactors[emdownscaling::kEnergy]; });           //! cluster energy (GeV)
-DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, [](int16_t storedEta) -> float { return storedEta / emdownscaling::downscalingFactors[emdownscaling::kEta]; });      //! cluster pseudorapidity
-DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, [](uint16_t storedPhi) -> float { return storedPhi / emdownscaling::downscalingFactors[emdownscaling::kPhi]; });     //! cluster azimuthal angle (0 to 2pi)
-DECLARE_SOA_DYNAMIC_COLUMN(NCells, nCells, [](int16_t storedNCells) -> int16_t { return storedNCells; });                                                 //! number of cells in cluster
-DECLARE_SOA_DYNAMIC_COLUMN(M02, m02, [](int16_t storedM02) -> float { return storedM02 / emdownscaling::downscalingFactors[emdownscaling::kM02]; });      //! shower shape long axis
-DECLARE_SOA_DYNAMIC_COLUMN(Time, time, [](int16_t storedTime) -> float { return storedTime / emdownscaling::downscalingFactors[emdownscaling::kTime]; }); //! cluster time (ns)
-DECLARE_SOA_DYNAMIC_COLUMN(IsExotic, isExotic, [](bool storedIsExotic) -> bool { return storedIsExotic; });                                               //! flag to mark cluster as exotic
+DECLARE_SOA_DYNAMIC_COLUMN(Definition, definition, [](int8_t storedDefinition) -> int8_t { return storedDefinition; });                                                                           //! cluster definition, see EMCALClusterDefinition.h
+DECLARE_SOA_DYNAMIC_COLUMN(E, e, [](int16_t storedE) -> float { return storedE / emdownscaling::downscalingFactors[emdownscaling::kEnergy] + std::numeric_limits<float>::epsilon(); });           //! cluster energy (GeV)
+DECLARE_SOA_DYNAMIC_COLUMN(Eta, eta, [](int16_t storedEta) -> float { return storedEta / emdownscaling::downscalingFactors[emdownscaling::kEta] + std::numeric_limits<float>::epsilon(); });      //! cluster pseudorapidity
+DECLARE_SOA_DYNAMIC_COLUMN(Phi, phi, [](uint16_t storedPhi) -> float { return storedPhi / emdownscaling::downscalingFactors[emdownscaling::kPhi] + std::numeric_limits<float>::epsilon(); });     //! cluster azimuthal angle (0 to 2pi)
+DECLARE_SOA_DYNAMIC_COLUMN(NCells, nCells, [](int16_t storedNCells) -> int16_t { return storedNCells; });                                                                                         //! number of cells in cluster
+DECLARE_SOA_DYNAMIC_COLUMN(M02, m02, [](int16_t storedM02) -> float { return storedM02 / emdownscaling::downscalingFactors[emdownscaling::kM02] + std::numeric_limits<float>::epsilon(); });      //! shower shape long axis
+DECLARE_SOA_DYNAMIC_COLUMN(Time, time, [](int16_t storedTime) -> float { return storedTime / emdownscaling::downscalingFactors[emdownscaling::kTime] + std::numeric_limits<float>::epsilon(); }); //! cluster time (ns)
+DECLARE_SOA_DYNAMIC_COLUMN(IsExotic, isExotic, [](bool storedIsExotic) -> bool { return storedIsExotic; });                                                                                       //! flag to mark cluster as exotic
 
 DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](float storedE, float storedEta) -> float { return storedE / emdownscaling::downscalingFactors[emdownscaling::kEnergy] / std::cosh(storedEta / emdownscaling::downscalingFactors[emdownscaling::kEta]); }); //! cluster pt, assuming m=0 (photons)
 } // namespace bcwisecluster
@@ -120,7 +122,7 @@
 DECLARE_SOA_COLUMN(IsPrimary, isPrimary, bool);   //! mcParticle.isPhysicalPrimary() || mcParticle.producedByGenerator()
 DECLARE_SOA_COLUMN(IsFromWD, isFromWD, bool);     //! Pi0 from a weak decay according to pwgem::photonmeson::utils::mcutil::IsFromWD
 
-DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](uint16_t storedpt) -> float { return storedpt / emdownscaling::downscalingFactors[emdownscaling::kpT]; }); //! pT of pi0 (GeV)
+DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, [](uint16_t storedpt) -> float { return storedpt / emdownscaling::downscalingFactors[emdownscaling::kpT] + std::numeric_limits<float>::epsilon(); }); //! pT of pi0 (GeV)
 } // namespace bcwisemcpi0s
 
 DECLARE_SOA_TABLE(BCWiseMCPi0s, "AOD", "BCWISEMCPI0", //! table of pi0s on MC level
@@ -132,7 +134,7 @@
 DECLARE_SOA_COLUMN(Pi0ID, pi0ID, int32_t);              //! Index of the mother pi0 (-1 if not from pi0)
 DECLARE_SOA_COLUMN(StoredTrueE, storedTrueE, uint16_t); //! energy of cluster inducing particle (1 MeV -> Maximum cluster energy of ~65 GeV)
 
-DECLARE_SOA_DYNAMIC_COLUMN(TrueE, trueE, [](uint16_t storedTrueE) -> float { return storedTrueE / emdownscaling::downscalingFactors[emdownscaling::kEnergy]; }); //! energy of cluster inducing particle (GeV)
+DECLARE_SOA_DYNAMIC_COLUMN(TrueE, trueE, [](uint16_t storedTrueE) -> float { return storedTrueE / emdownscaling::downscalingFactors[emdownscaling::kEnergy] + std::numeric_limits<float>::epsilon(); }); //! energy of cluster inducing particle (GeV)
 } // namespace bcwisemccluster
 
 DECLARE_SOA_TABLE(BCWiseMCClusters, "AOD", "BCWISEMCCLS", //! table of MC information for clusters -> To be joined with the cluster table

PWGEM/PhotonMeson/TableProducer/bcWiseClusterSkimmer.cxx

@@ -107,7 +107,7 @@
   template <typename OutputType, typename InputType>
   OutputType convertForStorage(InputType const& valueIn, Observables observable)
   {
-    double valueToBeChecked = valueIn * downscalingFactors[observable];
+    double valueToBeChecked = std::round(valueIn * downscalingFactors[observable]);
     if (valueToBeChecked < std::numeric_limits<OutputType>::lowest()) {
       LOG(warning) << "Value " << valueToBeChecked << " of observable " << observable << " below lowest possible value of " << typeid(OutputType).name() << ": " << static_cast<float>(std::numeric_limits<OutputType>::lowest());
       valueToBeChecked = static_cast<float>(std::numeric_limits<OutputType>::lowest());
@@ -205,10 +205,10 @@
   bool isGammaGammaDecay(TMCParticle mcParticle, TMCParticles mcParticles)
   {
     auto daughtersIds = mcParticle.daughtersIds();
     if (daughtersIds.size() != 2)
       return false;
     for (const auto& daughterId : daughtersIds) {
       if (mcParticles.iteratorAt(daughterId).pdgCode() != 22)
         return false;
     }
     return true;
@@ -218,10 +218,10 @@
   bool isAccepted(TMCParticle mcParticle, TMCParticles mcParticles)
   {
     auto daughtersIds = mcParticle.daughtersIds();
     if (daughtersIds.size() != 2)
       return false;
     for (const auto& daughterId : daughtersIds) {
       if (mcParticles.iteratorAt(daughterId).pdgCode() != 22)
         return false;
       int iCellID = -1;
       try {
@@ -273,7 +273,7 @@
       for (const auto& mcCollision : mcCollisionsBC) {
         auto mcParticlesInColl = mcParticles.sliceBy(perMcCollision, mcCollision.globalIndex());
         for (const auto& mcParticle : mcParticlesInColl) {
           if (mcParticle.pdgCode() != 111 || std::abs(mcParticle.y()) > cfgRapidityCut || !isGammaGammaDecay(mcParticle, mcParticles) || mcParticle.pt() < cfgMinPtGenPi0)
             continue;
           bool isPrimary = mcParticle.isPhysicalPrimary() || mcParticle.producedByGenerator();
           bool isFromWD = (aod::pwgem::photonmeson::utils::mcutil::IsFromWD(mcCollision, mcParticle, mcParticles)) > 0;

PWGEM/PhotonMeson/Tasks/emcalBcWisePi0.cxx

@@ -104,6 +104,7 @@ struct EmcalBcWisePi0 {
       mHistManager.add("Generated/pi0_AllBCs", "pT spectrum of generated pi0s in all BCs;#bf{#it{p}_{T} (GeV/#it{c})};#bf{#it{N}_{#pi^{0}}^{gen}}", HistType::kTH1F, {{200, 0, 20}});
       mHistManager.add("Generated/pi0_TVX", "pT spectrum of generated pi0s in TVX triggered BCs;#bf{#it{p}_{T} (GeV/#it{c})};#bf{#it{N}_{#pi^{0}}^{gen}}", HistType::kTH1F, {{200, 0, 20}});
       mHistManager.add("Generated/pi0_kTVXinEMC", "pT spectrum of generated pi0s in kTVXinEMC triggered BCs;#bf{#it{p}_{T} (GeV/#it{c})};#bf{#it{N}_{#pi^{0}}^{gen}}", HistType::kTH1F, {{200, 0, 20}});
+      mHistManager.add("Accepted/pi0_kTVXinEMC", "pT spectrum of accepted pi0s in kTVXinEMC triggered BCs;#bf{#it{p}_{T} (GeV/#it{c})};#bf{#it{N}_{#pi^{0}}^{acc}}", HistType::kTH1F, {{200, 0, 20}});
     }
   }
 
@@ -273,11 +274,15 @@ struct EmcalBcWisePi0 {
   void fillGeneratedPi0Hists(const auto& mcPi0s, const auto& bc)
   {
     for (const auto& mcPi0 : mcPi0s) {
-      mHistManager.fill(HIST("Generated/pi0_AllBCs"), mcPi0.pt());
-      if (bc.hasTVX())
-        mHistManager.fill(HIST("Generated/pi0_TVX"), mcPi0.pt());
-      if (bc.haskTVXinEMC())
-        mHistManager.fill(HIST("Generated/pi0_kTVXinEMC"), mcPi0.pt());
+      if (mcPi0.isPrimary()) {
+        mHistManager.fill(HIST("Generated/pi0_AllBCs"), mcPi0.pt());
+        if (bc.hasTVX())
+          mHistManager.fill(HIST("Generated/pi0_TVX"), mcPi0.pt());
+        if (bc.haskTVXinEMC())
+          mHistManager.fill(HIST("Generated/pi0_kTVXinEMC"), mcPi0.pt());
+        if (mcPi0.isAccepted() && bc.haskTVXinEMC())
+          mHistManager.fill(HIST("Accepted/pi0_kTVXinEMC"), mcPi0.pt());
+      }
     }
   }