[PWGEM] Add option to separate EMCal and DCal cluster pairing for flow task

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/Tasks/emcalPhotonMcTask.cxx

@@ -27,12 +27,14 @@
 
 #include <CCDB/BasicCCDBManager.h>
 #include <CommonConstants/MathConstants.h>
+#include <CommonConstants/PhysicsConstants.h>
 #include <DataFormatsParameters/GRPMagField.h>
 #include <DataFormatsParameters/GRPObject.h>
 #include <Framework/ASoA.h>
 #include <Framework/ASoAHelpers.h>
 #include <Framework/AnalysisHelpers.h>
 #include <Framework/AnalysisTask.h>
+#include <Framework/Concepts.h>
 #include <Framework/Configurable.h>
 #include <Framework/HistogramRegistry.h>
 #include <Framework/HistogramSpec.h>
@@ -51,7 +53,12 @@
 #include <TPDGCode.h>
 #include <TTree.h>
 
+#include <sys/types.h>
+
+#include <algorithm>
 #include <cmath>
+#include <cstddef>
+#include <cstdint>
 #include <random>
 #include <string>
 #include <utility>
@@ -148,14 +155,16 @@ ClusterMcInfo classifyCluster(const TGroup& g, TIter& mcCluster, TIter& mcCluste
     info.isPhoton = true;
   }
 
-  info.photonId = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster, mcParticles, std::vector<int>{PDG_t::kPi0, Pdg::kEta});
+  info.photonId = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster, mcParticles, std::vector<int>{PDG_t::kPi0, Pdg::kEta, Pdg::kOmega, Pdg::kEtaPrime});
   info.isFromPi0 = info.photonId >= 0;
 
   return info;
 }
+
 struct EmcalPhotonMcTask {
   static constexpr float EMCALRadius = 440.f;
   static constexpr float PhiVUndefined = -999.f;
+  static constexpr float Epsilon = 1.e-6f;
 
   Produces<aod::ConvTagCandidates> convTagCandidates;
 
@@ -241,7 +250,7 @@ struct EmcalPhotonMcTask {
   using McColls = o2::soa::Join<o2::aod::EMMCEvents, o2::aod::BinnedGenPts>;
   using McParticles = EMMCParticles;
 
-  PresliceOptional<EMCalPhotons> perCollisionEMC = o2::aod::emccluster::pmeventId;
+  PresliceOptional<aod::EMCEMEventIds> perCollisionEMC = o2::aod::emccluster::pmeventId;
   PresliceOptional<MinMTracks> perEMCClusterMT = o2::aod::mintm::minClusterId;
   PresliceOptional<MinMSTracks> perEMCClusterMS = o2::aod::mintm::minClusterId;
 
@@ -515,7 +524,7 @@ struct EmcalPhotonMcTask {
         ROOT::Math::XYZVector nRef = u.Cross(zAxis); // normal to the plane containing u and z
 
         float phiV = PhiVUndefined; // sentinel for degenerate geometry
-        if (nDecay.R() > 1e-6f && nRef.R() > 1e-6f) {
+        if (nDecay.R() > Epsilon && nRef.R() > Epsilon) {
           float cosPhiV = static_cast<float>(nDecay.Unit().Dot(nRef.Unit()));
           cosPhiV = std::clamp(cosPhiV, -1.f, 1.f);
           phiV = std::acos(cosPhiV);
@@ -627,7 +636,7 @@ struct EmcalPhotonMcTask {
       }
 
       mcCluster1.setCursor(cluster.emmcparticleIds()[0]);
-      int photonid1 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster1, mcParticles, std::vector<int>{PDG_t::kPi0, Pdg::kEta});
+      int photonid1 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcCluster1, mcParticles, std::vector<int>{PDG_t::kPi0, Pdg::kEta, Pdg::kOmega, Pdg::kEtaPrime});
       int motherId = -1;
       if (photonid1 >= 0) {
         mcPhoton1.setCursor(photonid1);

PWGEM/PhotonMeson/Tasks/taskFlowReso.cxx

@@ -50,7 +50,7 @@ using namespace o2::aod::pwgem::photon;
 
 enum QvecEstimator {
   FT0M = 0,
-  FT0A = 1,
+  FT0A,
   FT0C,
   TPCPos,
   TPCNeg,
@@ -60,7 +60,7 @@ enum QvecEstimator {
 
 enum CentralityEstimator {
   None = 0,
-  CFT0A = 1,
+  CFT0A,
   CFT0C,
   CFT0M,
   NCentralityEstimators
@@ -382,7 +382,7 @@ struct TaskFlowReso {
       // no selection on the centrality is applied on purpose to allow for the resolution study in post-processing
       return;
     }
-    if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) {
+    if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange()) || !(collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) {
       return;
     }
     float cent = getCentrality(collision);
@@ -530,7 +530,7 @@ struct TaskFlowReso {
 
 }; // End struct TaskFlowReso
 
-WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+WorkflowSpec defineDataProcessing(ConfigContext const& context)
 {
-  return WorkflowSpec{adaptAnalysisTask<TaskFlowReso>(cfgc)};
+  return WorkflowSpec{adaptAnalysisTask<TaskFlowReso>(context)};
 }

PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx

@@ -38,6 +38,7 @@
 #include <Framework/AnalysisHelpers.h>
 #include <Framework/AnalysisTask.h>
 #include <Framework/BinningPolicy.h>
+#include <Framework/Concepts.h>
 #include <Framework/Configurable.h>
 #include <Framework/Expressions.h>
 #include <Framework/GroupedCombinations.h>
@@ -172,6 +173,7 @@ struct TaskPi0FlowEMC {
     Configurable<bool> cfgEMCUseTM{"cfgEMCUseTM", false, "flag to use EMCal track matching cut or not"};
     Configurable<bool> emcUseSecondaryTM{"emcUseSecondaryTM", false, "flag to use EMCal secondary track matching cut or not"};
     Configurable<bool> cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"};
+    Configurable<bool> separateEMCalDCal{"separateEMCalDCal", false, "flag to only pair EMCal with EMCal and DCal with DCal clusters"};
   } emccuts;
 
   V0PhotonCut fV0PhotonCut;
@@ -251,6 +253,8 @@ struct TaskPi0FlowEMC {
   int runNow = 0;
   int runBefore = -1;
 
+  static constexpr float MaxPhiEMCal = 3.5f;
+
   // Filter clusterFilter = aod::skimmedcluster::time >= emccuts.cfgEMCminTime && aod::skimmedcluster::time <= emccuts.cfgEMCmaxTime && aod::skimmedcluster::m02 >= emccuts.cfgEMCminM02 && aod::skimmedcluster::m02 <= emccuts.cfgEMCmaxM02 && aod::skimmedcluster::e >= emccuts.cfgEMCminE;
   Filter collisionFilter = (nabs(aod::collision::posZ) <= eventcuts.cfgZvtxMax) && (aod::evsel::ft0cOccupancyInTimeRange <= eventcuts.cfgFT0COccupancyMax) && (aod::evsel::ft0cOccupancyInTimeRange >= eventcuts.cfgFT0COccupancyMin);
   // using FilteredEMCalPhotons = soa::Filtered<soa::Join<aod::EMCEMEventIds, aod::MinClusters>>;
@@ -260,10 +264,13 @@ struct TaskPi0FlowEMC {
   using CollsWithQvecs = soa::Join<aod::PMEvents, aod::EMEventsAlias, aod::EMEventsMult_000, aod::EMEventsCent_000, aod::EMEventsQvec_001>;
   using Colls = soa::Join<aod::PMEvents, aod::EMEventsAlias, aod::EMEventsMult_000, aod::EMEventsCent_000>;
 
+  Partition<EMCalPhotons> emcalPhotons = aod::mincluster::storedPhi < std::lround(MaxPhiEMCal * emcdownscaling::downscalingFactors[emcdownscaling::kPhi]);
+  Partition<EMCalPhotons> dcalPhotons = aod::mincluster::storedPhi >= std::lround(MaxPhiEMCal * emcdownscaling::downscalingFactors[emcdownscaling::kPhi]);
+
   static constexpr std::size_t NQVecEntries = 6;
 
-  PresliceOptional<EMCalPhotons> perCollisionEMC = o2::aod::emccluster::pmeventId;
-  PresliceOptional<PCMPhotons> perCollisionPCM = aod::v0photonkf::pmeventId;
+  PresliceOptional<aod::EMCEMEventIds> perCollisionEMC = o2::aod::emccluster::pmeventId;
+  PresliceOptional<aod::V0KFEMEventIds> perCollisionPCM = aod::v0photonkf::pmeventId;
   PresliceOptional<MinMTracks> perEMCClusterMT = o2::aod::mintm::minClusterId;
   PresliceOptional<MinMSTracks> perEMCClusterMS = o2::aod::mintm::minClusterId;
 
@@ -481,6 +488,13 @@ struct TaskPi0FlowEMC {
 
   }; // end init
 
+  /// \brief Check whether a photon (by its phi) falls in the EMCal or DCal acceptance
+  /// \param phi azimuthal angle of the photon
+  static bool isEMCalRegion(float phi)
+  {
+    return phi < MaxPhiEMCal;
+  }
+
   /// Change radians to degree
   /// \param angle in radians
   /// \return angle in degree
@@ -1047,7 +1061,7 @@ struct TaskPi0FlowEMC {
         continue;
       }
       if (rotationConfig.cfgDoRotation.value && nColl % rotationConfig.cfgDownsampling == 0) {
-        rotationBackground<EMCalPhotons>(vMeson, v1, v2, photons1, g1.globalIndex(), g2.globalIndex(), collision);
+        rotationBackground<TPhotons1>(vMeson, v1, v2, photons1, g1.globalIndex(), g2.globalIndex(), collision);
       }
       if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) {
         registry.fill(HIST("hMesonCuts"), 3);
@@ -1093,6 +1107,8 @@ struct TaskPi0FlowEMC {
       }
 
       auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex());
+      auto emcalPhotonsPerCollision = emcalPhotons.sliceBy(perCollisionEMC, collision.globalIndex());
+      auto dcalPhotonsPerCollision = dcalPhotons.sliceBy(perCollisionEMC, collision.globalIndex());
 
       if (emccuts.cfgEnableQA.value) {
         for (const auto& photon : photonsPerCollision) {
@@ -1108,7 +1124,12 @@ struct TaskPi0FlowEMC {
           registry.fill(HIST("clusterQA/hClusterEtaPhiAfter"), photon.phi(), photon.eta()); // after cuts
         }
       }
-      runPairingLoop(collision, photonsPerCollision, photonsPerCollision, flags, flags);
+      if (emccuts.separateEMCalDCal.value) {
+        runPairingLoop(collision, emcalPhotonsPerCollision, emcalPhotonsPerCollision, flags, flags);
+        runPairingLoop(collision, dcalPhotonsPerCollision, dcalPhotonsPerCollision, flags, flags);
+      } else {
+        runPairingLoop(collision, photonsPerCollision, photonsPerCollision, flags, flags);
+      }
       if (rotationConfig.cfgDoRotation.value) {
         if (nColl % rotationConfig.cfgDownsampling == 0) {
           nColl = 1; // reset counter
@@ -1163,6 +1184,10 @@ struct TaskPi0FlowEMC {
         if (!(flags.test(g1.globalIndex())) || !(flags.test(g2.globalIndex()))) {
           continue;
         }
+        if (emccuts.separateEMCalDCal.value && isEMCalRegion(g1.phi()) != isEMCalRegion(g2.phi())) {
+          continue; // only pair EMCal-EMCal or DCal-DCal
+        }
+
         // Cut edge clusters away, similar to rotation method to ensure same acceptance is used
         if (cfgDistanceToEdge.value > 0) {
           if (checkEtaPhi1D(g1.eta(), RecoDecay::constrainAngle(g1.phi())) >= cfgEMCalMapLevelBackground.value) {

PWGEM/PhotonMeson/Tasks/testTaskEmc.cxx

@@ -172,7 +172,7 @@ struct TestTaskEmc {
   PROCESS_SWITCH(TestTaskEmc, processEMCalCalib, "Process EMCal calibration same event", true);
 }; // End struct TestTaskEmc
 
-WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+WorkflowSpec defineDataProcessing(ConfigContext const& context)
 {
-  return WorkflowSpec{adaptAnalysisTask<TestTaskEmc>(cfgc)};
+  return WorkflowSpec{adaptAnalysisTask<TestTaskEmc>(context)};
 }

PWGEM/PhotonMeson/Utils/MCUtilities.h

@@ -17,6 +17,7 @@
 #define PWGEM_PHOTONMESON_UTILS_MCUTILITIES_H_
 
 #include <Framework/ASoA.h>
+#include <Framework/Concepts.h>
 
 #include <TPDGCode.h>