[PWGEM/PhotonMeson] Add MC functionality 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.
 //
@@ -37,7 +37,6 @@
   kFT0MCent,
   kZVtx,
   kFT0Amp,
-  kCellAmpSum,
   kpT,
   nObservables
 };
@@ -53,27 +52,26 @@
   1E2,  // Cluster time
   2E0,  // FT0M centrality
   1E3,  // Z-vertex position
-  1E0,  // FT0M amplitude
-  1E0,  // Cell energy
+  1E-1, // FT0M amplitude
   1E3}; // MC pi0 pt
 } // namespace emdownscaling
 
 namespace bcwisebc
 {
-DECLARE_SOA_COLUMN(HasFT0, hasFT0, bool);                                   //! has_foundFT0()
-DECLARE_SOA_COLUMN(HasTVX, hasTVX, bool);                                   //! has the TVX trigger flag
-DECLARE_SOA_COLUMN(HaskTVXinEMC, haskTVXinEMC, bool);                       //! kTVXinEMC
-DECLARE_SOA_COLUMN(HasEMCCell, hasEMCCell, bool);                           //! at least one EMCal cell in the BC
-DECLARE_SOA_COLUMN(HasNoTFROFBorder, hasNoTFROFBorder, bool);               //! not in the TF border or ITS ROF border region
-DECLARE_SOA_COLUMN(StoredFT0MAmplitude, storedFT0MAmplitude, unsigned int); //! ft0a+c amplitude
-DECLARE_SOA_COLUMN(StoredEMCalnCells, storedEMCalnCells, unsigned int);     //! number of emcal cells
-DECLARE_SOA_COLUMN(StoredEMCalCellEnergy, storedEMCalCellEnergy, float);    //! sum of energy in emcal cells
+DECLARE_SOA_COLUMN(HasFT0, hasFT0, bool);                               //! has_foundFT0()
+DECLARE_SOA_COLUMN(HasTVX, hasTVX, bool);                               //! has the TVX trigger flag
+DECLARE_SOA_COLUMN(HaskTVXinEMC, haskTVXinEMC, bool);                   //! kTVXinEMC
+DECLARE_SOA_COLUMN(HasEMCCell, hasEMCCell, bool);                       //! at least one EMCal cell in the BC
+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::StoredEMCalnCells, bcwisebc::StoredEMCalCellEnergy);
+                  bcwisebc::StoredFT0MAmplitude, bcwisebc::FT0MAmplitude<bcwisebc::StoredFT0MAmplitude>);
 
 DECLARE_SOA_INDEX_COLUMN(BCWiseBC, bcWiseBC); //! bunch crossing ID used as index
 
 namespace bcwisecollision
 {
@@ -117,29 +115,29 @@
 
 namespace bcwisemcpi0s
 {
-DECLARE_SOA_COLUMN(ParticleIdPi0, particleIdPi0, int); //! ID of the pi0 in the MC stack
-DECLARE_SOA_COLUMN(StoredPt, storedPt, uint16_t);      //! Transverse momentum of generated pi0 (10 MeV)
-DECLARE_SOA_COLUMN(IsAccepted, isAccepted, bool);      //! Both decay photons are within the EMCal acceptance
-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_COLUMN(StoredPt, storedPt, uint16_t); //! Transverse momentum of generated pi0 (1 MeV -> Maximum pi0 pT of ~65 GeV)
+DECLARE_SOA_COLUMN(IsAccepted, isAccepted, bool); //! Both decay photons are within the EMCal acceptance
+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)
 } // namespace bcwisemcpi0s
 
 DECLARE_SOA_TABLE(BCWiseMCPi0s, "AOD", "BCWISEMCPI0", //! table of pi0s on MC level
-                  o2::soa::Index<>, BCWiseBCId, bcwisemcpi0s::ParticleIdPi0, bcwisemcpi0s::StoredPt, bcwisemcpi0s::IsAccepted, bcwisemcpi0s::IsPrimary, bcwisemcpi0s::IsFromWD,
+                  o2::soa::Index<>, BCWiseBCId, bcwisemcpi0s::StoredPt, bcwisemcpi0s::IsAccepted, bcwisemcpi0s::IsPrimary, bcwisemcpi0s::IsFromWD,
                   bcwisemcpi0s::Pt<bcwisemcpi0s::StoredPt>);
 
 namespace bcwisemccluster
 {
-DECLARE_SOA_COLUMN(StoredE, storedE, uint16_t); //! energy of cluster inducing particle (1 MeV precision)
+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(E, e, [](uint16_t storedE) -> float { return storedE / 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]; }); //! 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
-                  o2::soa::Index<>, BCWiseBCId, bcwisemccluster::StoredE, bcwisemcpi0s::ParticleIdPi0,
-                  bcwisemccluster::E<bcwisemccluster::StoredE>);
+                  o2::soa::Index<>, BCWiseBCId, bcwisemccluster::Pi0ID, bcwisemccluster::StoredTrueE,
+                  bcwisemccluster::TrueE<bcwisemccluster::StoredTrueE>);
 
 } // namespace o2::aod
 

PWGEM/PhotonMeson/TableProducer/bcWiseClusterSkimmer.cxx

@@ -18,6 +18,7 @@
 
 #include <limits>
 #include <vector>
+#include <map>
 
 #include "Framework/runDataProcessing.h"
 #include "Framework/AnalysisTask.h"
@@ -63,6 +64,12 @@
   Configurable<float> cfgMinTime{"cfgMinTime", -25, "Minimum time of selected clusters (ns)"};
   Configurable<float> cfgMaxTime{"cfgMaxTime", 25, "Maximum time of selected clusters (ns)"};
   Configurable<float> cfgRapidityCut{"cfgRapidityCut", 0.8f, "Maximum absolute rapidity of counted generated particles"};
+  // Configurable<float> cfgMinPtGenPi0{"cfgMinPtGenPi0", 0., "Minimum pT for stored generated pi0s (reduce disk space of derived data)"};
+
+  Configurable<bool> cfgRequirekTVXinEMC{"cfgRequirekTVXinEMC", false, "Only store kTVXinEMC triggered BCs"};
+  Configurable<bool> cfgRequireGoodRCTQuality{"cfgRequireGoodRCTQuality", false, "Only store BCs with good quality of T0 and EMC in RCT"};
+
+  aod::rctsel::RCTFlagsChecker isFT0EMCGoodRCTChecker{aod::rctsel::kFT0Bad, aod::rctsel::kEMCBad};
 
   expressions::Filter energyFilter = aod::emcalcluster::energy > cfgMinClusterEnergy;
   expressions::Filter m02Filter = (aod::emcalcluster::nCells == 1 || (aod::emcalcluster::m02 > cfgMinM02 && aod::emcalcluster::m02 < cfgMaxM02));
@@ -71,6 +78,8 @@
 
   HistogramRegistry mHistManager{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false};
 
+  std::map<int32_t, int32_t> fMapPi0Index; // Map to connect the MC index of the pi0 to the one saved in the derived table
+
   void init(framework::InitContext&)
   {
     const int nEventBins = 6;
@@ -79,15 +88,13 @@
     for (int iBin = 0; iBin < nEventBins; iBin++)
       mHistManager.get<TH1>(HIST("nBCs"))->GetXaxis()->SetBinLabel(iBin + 1, binLabels[iBin]);
 
-    mHistManager.add("hAmplitudeVsCentrality", "FT0M AmplitudeVsCentrality;FT0M Centrality;FT0M Amplitude", HistType::kTH2F, {{105, 0, 105}, {400, 0, 200000}});
-    mHistManager.add("sumAmp", "Sum Amplitude", HistType::kTH1F, {{1000, 0., 10.}});
-    mHistManager.add("nCells", "Number of EMCal cells", HistType::kTH1F, {{5000, -0.5, 5000.5}});
-
     LOG(info) << "| Timing cut: " << cfgMinTime << " < t < " << cfgMaxTime;
     LOG(info) << "| M02 cut: " << cfgMinM02 << " < M02 < " << cfgMaxM02;
     LOG(info) << "| E cut: E > " << cfgMinClusterEnergy;
 
     o2::emcal::Geometry::GetInstanceFromRunNumber(300000);
+    if (cfgRequireGoodRCTQuality)
+      isFT0EMCGoodRCTChecker.init({aod::rctsel::kFT0Bad, aod::rctsel::kEMCBad});
   }
 
   /// \brief Process EMCAL clusters (either ambigous or unique)
@@ -129,28 +136,32 @@
   {
     for (const auto& cluster : clusters) {
       float clusterInducerEnergy = 0.;
-      int pi0MCIndex = -1;
+      int32_t pi0MCIndex = -1;
       if (cluster.amplitudeA().size() > 0) {
         int clusterInducerId = cluster.mcParticleIds()[0];
         auto clusterInducer = mcParticles.iteratorAt(clusterInducerId);
         clusterInducerEnergy = clusterInducer.e();
         int daughterId = aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(clusterInducer, mcParticles, std::vector<int>{111});
-        if (daughterId > 0) {
+        if (daughterId > 0)
           pi0MCIndex = mcParticles.iteratorAt(daughterId).mothersIds()[0];
-        }
       }
-      mcclusterTable(bcID, convertForStorage<uint16_t>(clusterInducerEnergy, kEnergy), pi0MCIndex);
+      if (pi0MCIndex > 0)
+        pi0MCIndex = fMapPi0Index[pi0MCIndex];
+      mcclusterTable(bcID, pi0MCIndex, convertForStorage<uint16_t>(clusterInducerEnergy, kEnergy));
     }
   }
 
-  void processEventProperties(const auto& bc, const auto& collisionsInBC, const auto& cellsInBC)
+  bool isBCSelected(const auto& bc)
   {
-    float sumAmp = 0.;
-    for (const auto& cell : cellsInBC)
-      sumAmp += cell.amplitude();
-    mHistManager.fill(HIST("sumAmp"), sumAmp);
-    mHistManager.fill(HIST("nCells"), cellsInBC.size());
+    if (cfgRequirekTVXinEMC && !bc.selection_bit(aod::evsel::kIsTriggerTVX))
+      return false;
+    if (cfgRequireGoodRCTQuality && !isFT0EMCGoodRCTChecker(bc))
+      return false;
+    return true;
+  }
 
+  void processEventProperties(const auto& bc, const auto& collisionsInBC, const auto& cellsInBC)
+  {
     bool hasFT0 = bc.has_foundFT0();
     bool hasTVX = bc.selection_bit(aod::evsel::kIsTriggerTVX);
     bool haskTVXinEMC = bc.alias_bit(kTVXinEMC);
@@ -170,25 +181,21 @@
 
     float ft0Amp = hasFT0 ? bc.foundFT0().sumAmpA() + bc.foundFT0().sumAmpC() : 0.;
 
-    bcTable(hasFT0, hasTVX, haskTVXinEMC, hasEMCCell, hasNoTFROFBorder, convertForStorage<unsigned int>(ft0Amp, kFT0Amp), convertForStorage<unsigned int>(cellsInBC.size(), kNCells), convertForStorage<float>(sumAmp, kCellAmpSum));
+    bcTable(hasFT0, hasTVX, haskTVXinEMC, hasEMCCell, hasNoTFROFBorder, convertForStorage<uint16_t>(ft0Amp, kFT0Amp));
 
     for (const auto& collision : collisionsInBC) {
       collisionTable(bcTable.lastIndex(), convertForStorage<uint8_t>(collision.centFT0M(), kFT0MCent), convertForStorage<int16_t>(collision.posZ(), kZVtx));
-      mHistManager.fill(HIST("hAmplitudeVsCentrality"), collision.centFT0M(), ft0Amp);
     }
-
-    if (collisionsInBC.size() == 0)
-      mHistManager.fill(HIST("hAmplitudeVsCentrality"), 103, ft0Amp);
   }
 
   template <typename TMCParticle, typename TMCParticles>
   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;
@@ -198,10 +205,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 {
@@ -218,6 +225,8 @@
   void processData(MyBCs const& bcs, MyCollisions const& collisions, aod::FT0s const&, SelectedCells const& cells, SelectedUniqueClusters const& uClusters, SelectedAmbiguousClusters const& aClusters)
   {
     for (const auto& bc : bcs) {
+      if (!isBCSelected(bc))
+        continue;
       auto collisionsInBC = collisions.sliceBy(perFoundBC, bc.globalIndex());
       auto cellsInBC = cells.sliceBy(cellsPerBC, bc.globalIndex());
 
@@ -240,6 +249,8 @@
   void processMC(MyBCs const& bcs, MyMCCollisions const& collisions, aod::McCollisions const& mcCollisions, aod::FT0s const&, SelectedCells const& cells, SelectedUniqueMCClusters const& uClusters, SelectedAmbiguousMCClusters const& aClusters, aod::McParticles const& mcParticles)
   {
     for (const auto& bc : bcs) {
+      if (!isBCSelected(bc))
+        continue;
       auto collisionsInBC = collisions.sliceBy(perFoundBC, bc.globalIndex());
       auto cellsInBC = cells.sliceBy(cellsPerBC, bc.globalIndex());
 
@@ -249,11 +260,12 @@
       for (const auto& mcCollision : mcCollisionsBC) {
         auto mcParticlesInColl = mcParticles.sliceBy(perMcCollision, mcCollision.globalIndex());
         for (const auto& mcParticle : mcParticlesInColl) {
-          if (mcParticle.pdgCode() != 111 || fabs(mcParticle.y()) > cfgRapidityCut || !isGammaGammaDecay(mcParticle, mcParticles))
+          if (mcParticle.pdgCode() != 111 || std::abs(mcParticle.y()) > cfgRapidityCut || !isGammaGammaDecay(mcParticle, mcParticles))
             continue;
           bool isPrimary = mcParticle.isPhysicalPrimary() || mcParticle.producedByGenerator();
           bool isFromWD = (aod::pwgem::photonmeson::utils::mcutil::IsFromWD(mcCollision, mcParticle, mcParticles)) > 0;
-          mcpi0Table(bc.globalIndex(), mcParticle.globalIndex(), convertForStorage<uint16_t>(mcParticle.pt(), kpT), isAccepted(mcParticle, mcParticles), isPrimary, isFromWD);
+          mcpi0Table(bc.globalIndex(), convertForStorage<uint16_t>(mcParticle.pt(), kpT), isAccepted(mcParticle, mcParticles), isPrimary, isFromWD);
+          fMapPi0Index[mcParticle.globalIndex()] = static_cast<int32_t>(mcpi0Table.lastIndex());
         }
       }
 
@@ -266,6 +278,7 @@
         processClusters(clustersInBC, bcTable.lastIndex());
         processClusterMCInfo(clustersInBC, bc.globalIndex(), mcParticles);
       }
+      fMapPi0Index.clear();
     }
   }
   PROCESS_SWITCH(bcWiseClusterSkimmer, processMC, "Run skimming for MC", false);