[PWGEM] Add option to separate EMCal and DCal cluster pairing for flow task (#17246)

Author: mhemmer-cern

PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h

@@ -47,7 +47,7 @@ enum Observable {
 };
 
 // Values in tables are stored in downscaled format to save disk space
-const std::array<float, nObservables> downscalingFactors{
+constexpr std::array<float, nObservables> downscalingFactors{
   1E0, // Cluster definition
   1E3, // Cluster energy
   1E4, // Cluster eta

PWGEM/PhotonMeson/DataModel/gammaTables.h

@@ -102,23 +102,31 @@ using V0LegMCLabel = V0LegMCLabels::iterator;
 // 3. EMEMCClusterMCLabels: Vector of EM MC particle ID of largest contributor to cluster
 namespace emcclustermclabel
 {
-DECLARE_SOA_ARRAY_INDEX_COLUMN(McParticle, emmcparticle) //!
+DECLARE_SOA_ARRAY_INDEX_COLUMN(McParticle, mcParticle)        //! Array of indice of McParticles
+DECLARE_SOA_COLUMN(Amplitude, amplitude, std::vector<float>); //! vector of aplitudes stored in cluster. Ordering is identical to the ordering of EMMCParticleId
 } // namespace emcclustermclabel
 
 // NOTE: MC labels. This table has one vector of global mc particle ids for each reconstructed emc cluster (joinable with emccluster table)
 DECLARE_SOA_TABLE(EMCClusterMCLabels, "AOD", "EMCClsMCLABEL", //!
                   emcclustermclabel::McParticleIds);
-using EMCClusterMCLabel = EMCClusterMCLabels::iterator;
+
+DECLARE_SOA_TABLE_VERSIONED(EMCClusterMCLabels_001, "AOD", "EMCClsMCLABEL", 1, //!
+                            emcclustermclabel::McParticleIds, emcclustermclabel::Amplitude);
+using EMCClusterMCLabel = EMCClusterMCLabels_001::iterator;
 
 namespace ememcclustermclabel
 {
-DECLARE_SOA_ARRAY_INDEX_COLUMN(EMMCParticle, emmcparticle); //!
+DECLARE_SOA_ARRAY_INDEX_COLUMN(EMMCParticle, emmcparticle);   //! Array of indice of EmMcParticles
+DECLARE_SOA_COLUMN(Amplitude, amplitude, std::vector<float>); //! vector of aplitudes stored in cluster. Ordering is identical to the ordering of EMMCParticleId
 } // namespace ememcclustermclabel
 
 // NOTE: MC labels. This table has a vector of entries for each reconstructed emc cluster (joinable with emccluster table)
 DECLARE_SOA_TABLE(EMEMCClusterMCLabels, "AOD", "EMEMCClsMCLABEL", //!
                   ememcclustermclabel::EMMCParticleIds);
-using EMEMCClusterMCLabel = EMEMCClusterMCLabels::iterator;
+
+DECLARE_SOA_TABLE_VERSIONED(EMEMCClusterMCLabels_001, "AOD", "EMEMCClsMCLABEL", 1, //!
+                            ememcclustermclabel::EMMCParticleIds, ememcclustermclabel::Amplitude);
+using EMEMCClusterMCLabel = EMEMCClusterMCLabels_001::iterator;
 
 namespace v0leg
 {
@@ -146,9 +154,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITS, meanClusterSizeITS, [](uint32_t i
   }
   if (nl > 0) {
     return static_cast<float>(total_cluster_size) / static_cast<float>(nl);
-  } else {
-    return 0;
   }
+  return 0;
 });
 DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSib, meanClusterSizeITSib, [](uint32_t itsClusterSizes) -> float {
   int total_cluster_size = 0, nl = 0;
@@ -161,9 +168,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSib, meanClusterSizeITSib, [](uint32
   }
   if (nl > 0) {
     return static_cast<float>(total_cluster_size) / static_cast<float>(nl);
-  } else {
-    return 0;
   }
+  return 0;
 });
 DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSob, meanClusterSizeITSob, [](uint32_t itsClusterSizes) -> float {
   int total_cluster_size = 0, nl = 0;
@@ -176,9 +182,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(MeanClusterSizeITSob, meanClusterSizeITSob, [](uint32
   }
   if (nl > 0) {
     return static_cast<float>(total_cluster_size) / static_cast<float>(nl);
-  } else {
-    return 0;
   }
+  return 0;
 });
 } // namespace v0leg
 DECLARE_SOA_TABLE(V0Legs_000, "AOD", "V0LEG", //!

PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

@@ -10,8 +10,8 @@
 // or submit itself to any jurisdiction.
 
 /// \file associateMCinfoPhoton.cxx
-/// \brief This code produces reduced events for photon analyses
-/// \author Daiki Sekihata (daiki.sekihata@cern.ch)
+/// \brief This code produces EmMc tables were the McParticleIds get reshuffled due to not storing all McParticles
+/// \author Daiki Sekihata (daiki.sekihata@cern.ch), Marvin Hemmer (marvin.hemmer@cern.ch), Nicolas Strangmann (nicolas.strangmann@cern.ch)
 
 #include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h"
 #include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
@@ -38,6 +38,7 @@
 #include <sys/types.h>
 
 #include <algorithm>
+#include <array>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
@@ -57,7 +58,7 @@ using namespace o2::constants::physics;
 using MyCollisionsMC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::EMEvSels>;
 using TracksMC = soa::Join<aod::TracksIU, aod::McTrackLabels>;
 using FwdTracksMC = soa::Join<aod::FwdTracks, aod::McFwdTrackLabels>;
-using MyEMCClusters = soa::Join<aod::MinClusters, aod::EMCClusterMCLabels>;
+using MyEMCClusters = soa::Join<aod::MinClusters, aod::EMCClusterMCLabels_001>;
 
 struct Counter {
   int particles{0};
@@ -77,7 +78,7 @@ struct AssociateMCInfoPhoton {
   Produces<o2::aod::EMMCParticles> emmcparticles;
   Produces<o2::aod::V0LegMCLabels> v0legmclabels;
   Produces<o2::aod::EMPrimaryElectronMCLabels> emprimaryelectronmclabels;
-  Produces<o2::aod::EMEMCClusterMCLabels> ememcclustermclabels;
+  Produces<o2::aod::EMEMCClusterMCLabels_001> ememcclustermclabels;
 
   Produces<o2::aod::BinnedGenPts> binnedGenPt;
 
@@ -100,6 +101,7 @@ struct AssociateMCInfoPhoton {
 
     // !!Don't change pt,eta,y binning. These binnings have to be consistent with binned data at analysis.!!
     std::vector<double> ptbins;
+    ptbins.reserve(72);
     for (int i = 0; i < 2; i++) {                // o2-linter: disable=magic-number (just numbers for binning)
       ptbins.emplace_back(0.05 * (i - 0) + 0.0); // from 0 to 0.05 GeV/c, every 0.05 GeV/c
     }
@@ -116,7 +118,7 @@ struct AssociateMCInfoPhoton {
     const AxisSpec axisRapidity{{0.0, +0.8, +0.9}, "rapidity |y|"};
 
     static constexpr uint NParticleNames = 9;
-    static constexpr std::string_view ParticleNames[NParticleNames] = {
+    static constexpr std::array<std::string_view, NParticleNames> ParticleNames = {
       "Gamma", "Pi0", "Eta", "Omega", "Phi",
       "ChargedPion", "ChargedKaon", "K0S", "Lambda"};
 
@@ -267,17 +269,19 @@ struct AssociateMCInfoPhoton {
               genGamma[binNumber]++;
               break;
             case PDG_t::kPi0:
-              if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles))
+              if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles)) {
                 continue;
+              }
               registry.fill(HIST("Generated/h2PtY_Pi0"), mcParticle.pt(), std::fabs(mcParticle.y()));
               genPi0[binNumber]++;
               if (isMesonAccepted) {
                 registry.fill(HIST("Generated/h2PtY_Accepted_Pi0"), mcParticle.pt(), std::fabs(mcParticle.y()));
               }
               break;
             case Pdg::kEta:
-              if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles))
+              if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles)) {
                 continue;
+              }
               registry.fill(HIST("Generated/h2PtY_Eta"), mcParticle.pt(), std::fabs(mcParticle.y()));
               genEta[binNumber]++;
               if (isMesonAccepted) {
@@ -291,7 +295,7 @@ struct AssociateMCInfoPhoton {
       } // end of mc track loop
 
       // make an entry for this MC event only if it was not already added to the table
-      if (!(fEventLabels.find(mcCollisionIter.globalIndex()) != fEventLabels.end())) {
+      if (!fEventLabels.contains(mcCollisionIter.globalIndex())) {
         mcevents(mcCollisionIter.globalIndex(), mcCollisionIter.generatorsID(), mcCollisionIter.posX(), mcCollisionIter.posY(), mcCollisionIter.posZ(), mcCollisionIter.impactParameter(), mcCollisionIter.eventPlaneAngle());
         fEventLabels[mcCollisionIter.globalIndex()] = fCounter.events;
         fCounter.events++;
@@ -438,15 +442,17 @@ struct AssociateMCInfoPhoton {
         mcCollisionIter.setCursor(collisionIter.mcCollisionId());
 
         // TODO: test
-        if (emccluster.emmcparticleIds().size() <= 0) {
+        if (emccluster.mcParticleIds().size() <= 0) {
           continue;
         }
         std::vector<int32_t> vEmcMcParticleIds;
+        std::vector<float> vAmplitudes;
 
-        vEmcMcParticleIds.reserve(emccluster.emmcparticleIds().size());
+        vEmcMcParticleIds.reserve(emccluster.mcParticleIds().size());
+        vAmplitudes.reserve(emccluster.mcParticleIds().size());
 
-        for (const auto& emcParticleId : emccluster.emmcparticleIds()) {
-          mcPhoton.setCursor(emcParticleId);
+        for (size_t iCont = 0; iCont < emccluster.mcParticleIds().size(); iCont++) {
+          mcPhoton.setCursor(emccluster.mcParticleIds()[iCont]);
 
           // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack
           auto [iter, isNew] = fNewLabels.try_emplace(mcPhoton.globalIndex(), fCounter.particles);
@@ -456,6 +462,7 @@ struct AssociateMCInfoPhoton {
             fCounter.particles++;
           }
           vEmcMcParticleIds.emplace_back(iter->second);
+          vAmplitudes.emplace_back(emccluster.amplitude()[iCont]);
           // ememcclustermclabels(fNewLabels.find(mcPhoton.index())->second);
 
           // Next, store mother-chain of this reconstructed track.
@@ -465,7 +472,7 @@ struct AssociateMCInfoPhoton {
           }
           selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter);
         } // end of loop over mc particles of the current emc cluster
-        ememcclustermclabels(vEmcMcParticleIds);
+        ememcclustermclabels(vEmcMcParticleIds, vAmplitudes);
 
       } // end of em emc cluster loop
     }
@@ -625,8 +632,8 @@ struct AssociateMCInfoPhoton {
   PROCESS_SWITCH(AssociateMCInfoPhoton, processDummy, "processDummy", true);
 };
 
-WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+WorkflowSpec defineDataProcessing(ConfigContext const& context)
 {
   return WorkflowSpec{
-    adaptAnalysisTask<AssociateMCInfoPhoton>(cfgc, TaskName{"associate-mc-info-photon"})};
+    adaptAnalysisTask<AssociateMCInfoPhoton>(context, TaskName{"associate-mc-info-photon"})};
 }

PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx

@@ -53,7 +53,7 @@ struct SkimmerGammaCalo {
   Preslice<o2::aod::EMCMatchSecs> psMSperCluster = o2::aod::emcalclustercell::emcalclusterId;
 
   Produces<aod::SkimEMCClusters_001> tableGammaEMCReco;
-  Produces<aod::EMCClusterMCLabels> tableEMCClusterMCLabels;
+  Produces<aod::EMCClusterMCLabels_001> tableEMCClusterMCLabels;
   Produces<aod::SkimEMCCells> tableCellEMCReco;
 
   Produces<aod::EmEmcClusters_000> tableEmEmcClusters;
@@ -331,15 +331,16 @@ struct SkimmerGammaCalo {
         continue;
       }
       std::vector<int32_t> mcLabels;
+      std::vector<float> amplitudes;
+      mcLabels.reserve(emccluster.amplitudeA().size());
+      amplitudes.reserve(emccluster.amplitudeA().size());
       for (size_t iCont = 0; iCont < emccluster.amplitudeA().size(); iCont++) {
         mcLabels.push_back(emccluster.mcParticleIds()[iCont]);
+        amplitudes.push_back(emccluster.amplitudeA()[iCont]);
       }
-      // LOGF(info, "---- New Cluster ---");
-      // for (unsigned long int iCont = 0; iCont < mcLabels.size(); iCont++) {
-      //   LOGF(info, "iCont = %d, mcParticle = %d, amplitudeA = %.5f", iCont, mcLabels.at(iCont), emccluster.amplitudeA()[iCont]);
-      // }
-      tableEMCClusterMCLabels(mcLabels);
+      tableEMCClusterMCLabels(mcLabels, amplitudes);
       mcLabels.clear();
+      amplitudes.clear();
     }
   }
   PROCESS_SWITCH(SkimmerGammaCalo, processMC, "process MC info", false); // Run this in addition to processRec for MCs to copy the cluster mc labels from the EMCALMCClusters to the skimmed EMCClusterMCLabels table