[PWGEM] PhotonMeson: Add basic EMCal QA to EMCPhotonCut (#14785)

Author: mhemmer-cern

PWGEM/PhotonMeson/Core/EMCPhotonCut.cxx

@@ -82,6 +82,12 @@ void EMCPhotonCut::SetUseSecondaryTM(bool flag)
   LOG(info) << "EM Photon Cluster Cut, using secondary TM cut is set to : " << mUseTM;
 }
 
+void EMCPhotonCut::SetDoQA(bool flag)
+{
+  mDoQA = flag;
+  LOG(info) << "EM Photon Cluster Cut, QA is set to: " << mUseTM;
+}
+
 void EMCPhotonCut::print() const
 {
   LOG(info) << "EMCal Photon Cut:";

PWGEM/PhotonMeson/Core/EMCPhotonCut.h

@@ -20,6 +20,7 @@
 #include "PWGEM/PhotonMeson/DataModel/gammaTables.h"
 
 #include <Framework/ASoA.h>
+#include <Framework/HistogramRegistry.h>
 
 #include <TNamed.h>
 
@@ -92,7 +93,7 @@ class EMCPhotonCut : public TNamed
   EMCPhotonCut() = default;
   EMCPhotonCut(const char* name, const char* title) : TNamed(name, title) {}
 
-  enum class EMCPhotonCuts : int {
+  enum class EMCPhotonCuts : std::uint8_t {
     // cluster cut
     kDefinition = 0,
     kEnergy,
@@ -105,6 +106,12 @@ class EMCPhotonCut : public TNamed
     kNCuts
   };
 
+  enum class TrackType : std::uint8_t {
+    // cluster cut
+    kPrimary = 0,
+    kSecondary,
+  };
+
   static const char* mCutNames[static_cast<int>(EMCPhotonCuts::kNCuts)];
 
   constexpr auto getClusterId(o2::soa::is_iterator auto const& t) const
@@ -126,7 +133,7 @@ class EMCPhotonCut : public TNamed
   /// \param GetPhiCut lambda to get the phi cut value
   /// \param applyEoverP bool to check if E/p should be checked (for secondaries we do not check this!)
   bool checkTrackMatching(o2::soa::is_iterator auto const& cluster, IsTrackIterator auto& emcmatchedtrack, o2::soa::RowViewSentinel const emcmatchedtrackEnd,
-                          bool applyEoverP, auto GetEtaCut, auto GetPhiCut) const
+                          bool applyEoverP, auto GetEtaCut, auto GetPhiCut, o2::framework::HistogramRegistry* fRegistry = nullptr, TrackType trackType = TrackType::kPrimary) const
   {
     // advance to cluster
     while (emcmatchedtrack != emcmatchedtrackEnd && getClusterId(emcmatchedtrack) < cluster.globalIndex()) {
@@ -150,29 +157,137 @@ class EMCPhotonCut : public TNamed
                   (dPhi > GetPhiCut(trackpt)) ||
                   (applyEoverP && cluster.e() / trackp >= mMinEoverP);
       if (!fail) {
+        if (mDoQA && fRegistry != nullptr) {
+          switch (trackType) {
+            case TrackType::kPrimary:
+              fRegistry->fill(HIST("Cluster/hTrackdEtadPhi"), dEta, dPhi);
+              fRegistry->fill(HIST("Cluster/hTrackdEtaPt"), dEta, trackpt);
+              fRegistry->fill(HIST("Cluster/hTrackdPhiPt"), dPhi, trackpt);
+              break;
+            case TrackType::kSecondary:
+              fRegistry->fill(HIST("Cluster/hSecTrackdEtadPhi"), dEta, dPhi);
+              fRegistry->fill(HIST("Cluster/hSecTrackdEtaPt"), dEta, trackpt);
+              fRegistry->fill(HIST("Cluster/hSecTrackdPhiPt"), dPhi, trackpt);
+              break;
+            default:
+              break;
+          }
+        }
         return false; // cluster got a track matche to it
       }
       ++emcmatchedtrack;
     }
     return true; // all tracks checked, cluster survives
   }
 
+  void fillBeforeClusterHistogram(o2::soa::is_iterator auto const& cluster, o2::framework::HistogramRegistry* fRegistry = nullptr)
+  {
+
+    if (mDoQA == false || fRegistry == nullptr) {
+      return;
+    }
+
+    fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), 0, cluster.e());
+    fRegistry->fill(HIST("Cluster/before/hE"), 0, cluster.e());
+    fRegistry->fill(HIST("Cluster/before/hPt"), 0, cluster.pt());
+    fRegistry->fill(HIST("Cluster/before/hEtaPhi"), 0, cluster.eta(), cluster.phi());
+    fRegistry->fill(HIST("Cluster/before/hNCell"), 0, cluster.nCells(), cluster.e());
+    fRegistry->fill(HIST("Cluster/before/hM02"), 0, cluster.m02(), cluster.e());
+    fRegistry->fill(HIST("Cluster/before/hTime"), 0, cluster.time(), cluster.e());
+  }
+
+  void fillAfterClusterHistogram(o2::soa::is_iterator auto const& cluster, o2::framework::HistogramRegistry* fRegistry = nullptr)
+  {
+
+    if (mDoQA == false || fRegistry == nullptr) {
+      return;
+    }
+
+    fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), 0, cluster.e());
+    fRegistry->fill(HIST("Cluster/before/hE"), 0, cluster.e());
+    fRegistry->fill(HIST("Cluster/before/hPt"), 0, cluster.pt());
+    fRegistry->fill(HIST("Cluster/before/hEtaPhi"), 0, cluster.eta(), cluster.phi());
+    fRegistry->fill(HIST("Cluster/before/hNCell"), 0, cluster.nCells(), cluster.e());
+    fRegistry->fill(HIST("Cluster/before/hM02"), 0, cluster.m02(), cluster.e());
+    fRegistry->fill(HIST("Cluster/before/hTime"), 0, cluster.time(), cluster.e());
+  }
+
   /// \brief check if given clusters survives all cuts
   /// \param flags EMBitFlags where results will be stored
   /// \param cluster cluster table to check
   /// \param matchedTracks matched primary tracks table
   /// \param matchedSecondaries matched secondary tracks table
-  void AreSelectedRunning(EMBitFlags& flags, o2::soa::is_table auto const& clusters, IsTrackContainer auto const& emcmatchedtracks, IsTrackContainer auto const& secondaries) const
+  /// \param fRegistry HistogramRegistry pointer of the main task
+  void AreSelectedRunning(EMBitFlags& flags, o2::soa::is_table auto const& clusters, IsTrackContainer auto const& emcmatchedtracks, IsTrackContainer auto const& secondaries, o2::framework::HistogramRegistry* fRegistry = nullptr)
   {
     auto emcmatchedtrackIter = emcmatchedtracks.begin();
     auto emcmatchedtrackEnd = emcmatchedtracks.end();
     auto secondaryIter = secondaries.begin();
     auto secondaryEnd = secondaries.end();
     size_t iCluster = 0;
+
+    // Check if we want to do QA and provided proper histogram registry
+    if (mDoQA && fRegistry != nullptr) {
+      const o2::framework::AxisSpec thAxisClusterEnergy{500, 0, 50, "#it{E}_{cls} (GeV)"};
+      const o2::framework::AxisSpec thAxisMomentum{250, 0., 25., "#it{p}_{T} (GeV/#it{c})"};
+      const o2::framework::AxisSpec thAxisDEta{200, -0.1, 0.1, "#Delta#eta"};
+      const o2::framework::AxisSpec thAxisDPhi{200, -0.1, 0.1, "#Delta#varphi (rad)"};
+      const o2::framework::AxisSpec thAxisEnergy{500, 0., 50., "#it{E} (GeV)"};
+      const o2::framework::AxisSpec thAxisEta{320, -0.8, 0.8, "#eta"};
+      const o2::framework::AxisSpec thAxisPhi{500, 0, o2::constants::math::TwoPI, "#varphi (rad)"};
+      const o2::framework::AxisSpec thAxisNCell{51, -0.5, 50.5, "#it{N}_{cell}"};
+      const o2::framework::AxisSpec thAxisM02{200, 0, 2.0, "#it{M}_{02}"};
+      const o2::framework::AxisSpec thAxisTime{300, -150, +150, "#it{t}_{cls} (ns)"};
+      const o2::framework::AxisSpec thAxisEoverP{400, 0, 10., "#it{E}_{cls}/#it{p}_{track} (#it{c})"};
+
+      fRegistry->add("Cluster/before/hE", "E_{cluster};#it{E}_{cluster} (GeV);#it{N}_{cluster}", o2::framework::kTH1F, {thAxisClusterEnergy}, true);
+      fRegistry->add("Cluster/before/hPt", "Transverse momenta of clusters;#it{p}_{T} (GeV/c);#it{N}_{cluster}", o2::framework::kTH1F, {thAxisClusterEnergy}, true);
+      fRegistry->add("Cluster/before/hNgamma", "Number of #gamma candidates per collision;#it{N}_{#gamma} per collision;#it{N}_{collisions}", o2::framework::kTH1F, {{1001, -0.5f, 1000.5f}}, true);
+      fRegistry->add("Cluster/before/hEtaPhi", "#eta vs #varphi;#eta;#varphi (rad.)", o2::framework::kTH2F, {thAxisEta, thAxisPhi}, true);
+      fRegistry->add("Cluster/before/hNCell", "#it{N}_{cells};N_{cells} (GeV);#it{E}_{cluster} (GeV)", o2::framework::kTH2F, {thAxisNCell, thAxisClusterEnergy}, true);
+      fRegistry->add("Cluster/before/hM02", "Long ellipse axis;#it{M}_{02} (cm);#it{E}_{cluster} (GeV)", o2::framework::kTH2F, {thAxisM02, thAxisClusterEnergy}, true);
+      fRegistry->add("Cluster/before/hTime", "Cluster time;#it{t}_{cls} (ns);#it{E}_{cluster} (GeV)", o2::framework::kTH2F, {thAxisTime, thAxisClusterEnergy}, true);
+
+      fRegistry->addClone("Cluster/before/", "Cluster/after/");
+
+      auto hClusterQualityCuts = fRegistry->add<TH2>("Cluster/hClusterQualityCuts", "Energy at which clusters are removed by a given cut;;#it{E} (GeV)", o2::framework::kTH2F, {{static_cast<int>(EMCPhotonCut::EMCPhotonCuts::kNCuts) + 2, -0.5, static_cast<double>(EMCPhotonCut::EMCPhotonCuts::kNCuts) + 1.5}, thAxisClusterEnergy}, true);
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(1, "In");
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(2, "Definition");
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(3, "Energy");
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(4, "NCell");
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(5, "M02");
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(6, "Timing");
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(7, "TM");
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(8, "Sec. TM");
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(9, "Exotic");
+      hClusterQualityCuts->GetXaxis()->SetBinLabel(10, "Out");
+
+      fRegistry->add("Cluster/hTrackdEtadPhi", "d#eta vs. d#varphi of matched tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDEta, thAxisDPhi}, true);
+      fRegistry->add("Cluster/hTrackdEtaPt", "d#eta vs. track pT of matched tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDEta, thAxisMomentum}, true);
+      fRegistry->add("Cluster/hTrackdPhiPt", "d#varphi vs. track pT of matched tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDPhi, thAxisMomentum}, true);
+      fRegistry->add("Cluster/hSecTrackdEtadPhi", "d#eta vs. d#varphi of matched secondary tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDEta, thAxisDPhi}, true);
+      fRegistry->add("Cluster/hSecTrackdEtaPt", "d#eta vs. track pT of matched secondary tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDEta, thAxisMomentum}, true);
+      fRegistry->add("Cluster/hSecTrackdPhiPt", "d#varphi vs. track pT of matched secondary tracks;d#eta;d#varphi (rad.)", o2::framework::kTH2F, {thAxisDPhi, thAxisMomentum}, true);
+    }
+
+    nTotClusterPerColl = 0;
+    currentCollID = clusters.iteratorAt(0).emeventId();
     for (const auto& cluster : clusters) {
+      if (currentCollID == cluster.emeventId()) {
+        ++nTotClusterPerColl;
+      } else {
+        fRegistry->fill(HIST("Cluster/before/hNgamma"), nTotClusterPerColl);
+        nTotClusterPerColl = 0;
+      }
+      if (mDoQA == true || fRegistry != nullptr) {
+        fillBeforeClusterHistogram(cluster, fRegistry);
+      }
       if (!IsSelectedRunning(cluster, emcmatchedtrackIter, emcmatchedtrackEnd, secondaryIter, secondaryEnd)) {
         flags.set(iCluster);
+      } else if (mDoQA == true || fRegistry != nullptr) {
+        fillAfterClusterHistogram(cluster, fRegistry);
       }
+      currentCollID = cluster.emeventId();
       ++iCluster;
     }
   }
@@ -184,32 +299,47 @@ class EMCPhotonCut : public TNamed
   /// \param secondaryIter current iterator of matched secondary tracks
   /// \param secondaryEnd end iterator of matched secondary tracks
   /// \return true if cluster survives all cuts else false
-  bool IsSelectedRunning(o2::soa::is_iterator auto const& cluster, IsTrackIterator auto& emcmatchedtrackIter, o2::soa::RowViewSentinel const emcmatchedtrackEnd, IsTrackIterator auto& secondaryIter, o2::soa::RowViewSentinel const secondaryEnd) const
+  bool IsSelectedRunning(o2::soa::is_iterator auto const& cluster, IsTrackIterator auto& emcmatchedtrackIter, o2::soa::RowViewSentinel const emcmatchedtrackEnd, IsTrackIterator auto& secondaryIter, o2::soa::RowViewSentinel const secondaryEnd, o2::framework::HistogramRegistry* fRegistry = nullptr)
   {
     if (!IsSelectedEMCalRunning(EMCPhotonCuts::kDefinition, cluster)) {
+      fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast<int>(EMCPhotonCuts::kDefinition) + 1, cluster.e());
       return false;
     }
     if (!IsSelectedEMCalRunning(EMCPhotonCuts::kEnergy, cluster)) {
+      fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast<int>(EMCPhotonCuts::kEnergy) + 1, cluster.e());
       return false;
     }
     if (!IsSelectedEMCalRunning(EMCPhotonCuts::kNCell, cluster)) {
+      fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast<int>(EMCPhotonCuts::kNCell) + 1, cluster.e());
       return false;
     }
     if (!IsSelectedEMCalRunning(EMCPhotonCuts::kM02, cluster)) {
+      fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast<int>(EMCPhotonCuts::kM02) + 1, cluster.e());
       return false;
     }
     if (!IsSelectedEMCalRunning(EMCPhotonCuts::kTiming, cluster)) {
+      fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast<int>(EMCPhotonCuts::kTiming) + 1, cluster.e());
       return false;
     }
     if (mUseTM && (!IsSelectedEMCalRunning(EMCPhotonCuts::kTM, cluster, emcmatchedtrackIter, emcmatchedtrackEnd))) {
+      fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast<int>(EMCPhotonCuts::kTM) + 1, cluster.e());
       return false;
     }
     if (mUseSecondaryTM && (!IsSelectedEMCalRunning(EMCPhotonCuts::kSecondaryTM, cluster, secondaryIter, secondaryEnd))) {
+      fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast<int>(EMCPhotonCuts::kSecondaryTM) + 1, cluster.e());
       return false;
     }
     if (!IsSelectedEMCalRunning(EMCPhotonCuts::kExotic, cluster)) {
+      fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast<int>(EMCPhotonCuts::kExotic) + 1, cluster.e());
       return false;
     }
+    fRegistry->fill(HIST("Cluster/hClusterQualityCuts"), static_cast<int>(EMCPhotonCuts::kNCuts) + 1, cluster.e());
+    if (currentCollID == cluster.emeventId()) {
+      ++nAccClusterPerColl;
+    } else {
+      fRegistry->fill(HIST("Cluster/before/hNgamma"), nAccClusterPerColl);
+      nAccClusterPerColl = 0;
+    }
     return true;
   }
 
@@ -439,6 +569,10 @@ class EMCPhotonCut : public TNamed
   /// \param flag flag to use secondary track matching
   void SetUseSecondaryTM(bool flag = false);
 
+  /// \brief Set flag to do QA
+  /// \param flag flag to do QA
+  void SetDoQA(bool flag = false);
+
   /// \brief Set parameters for track matching delta eta = a + (pT + b)^c
   /// \param a a in a + (pT + b)^c
   /// \param b b in a + (pT + b)^c
@@ -517,8 +651,13 @@ class EMCPhotonCut : public TNamed
   float mMaxTime{25.f};        ///< maximum cluster timing
   float mMinEoverP{1.75f};     ///< minimum cluster energy over track momentum ratio needed for the pair to be considered matched
   bool mUseExoticCut{true};    ///< flag to decide if the exotic cluster cut is to be checked or not
-  bool mUseTM{true};           ///< flag to decide if track matching cut is to be checek or not
-  bool mUseSecondaryTM{false}; ///< flag to decide if seconary track matching cut is to be checek or not
+  bool mUseTM{true};           ///< flag to decide if track matching cut is to be check or not
+  bool mUseSecondaryTM{false}; ///< flag to decide if seconary track matching cut is to be check or not
+  bool mDoQA{false};           ///< flag to decide if QA should be done or not
+
+  uint nTotClusterPerColl{0}; ///< running number of all cluster per collision used for QA
+  uint nAccClusterPerColl{0}; ///< running number of accepted cluster per collision used for QA
+  int currentCollID{-1};      ///< running collision ID of clusters used for QA
 
   TrackMatchingParams mTrackMatchingEtaParams = {-1.f, 0.f, 0.f};
   TrackMatchingParams mTrackMatchingPhiParams = {-1.f, 0.f, 0.f};

PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx

@@ -702,7 +702,7 @@ struct CalibTaskEmc {
     int nColl = 1;
 
     EMBitFlags emcFlags(clusters.size());
-    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds);
+    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, &registry);
 
     for (const auto& collision : collisions) {
       auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex());
@@ -812,7 +812,7 @@ struct CalibTaskEmc {
     int nColl = 1;
 
     EMBitFlags emcFlags(clusters.size());
-    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds);
+    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, &registry);
 
     for (const auto& collision : collisions) {
       o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(&registry, collision);
@@ -917,7 +917,7 @@ struct CalibTaskEmc {
   {
     float energyCorrectionFactor = 1.f;
     EMBitFlags emcFlags(clusters.size());
-    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds);
+    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, &registry);
 
     SameKindPair<Colls, EMCalPhotons, BinningType> pair{binningOnPositions, mixingConfig.cfgMixingDepth, -1, &cache}; // indicates that 5 events should be mixed and under/overflow (-1) to be ignored
 
@@ -998,7 +998,7 @@ struct CalibTaskEmc {
   {
     float energyCorrectionFactor = 1.f;
     EMBitFlags emcFlags(clusters.size());
-    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds);
+    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, &registry);
 
     for (const auto& [c1, photonEMC, c2, photonPCM] : pairPCMEMC) {
       if (!(fEMEventCut.IsSelected(c1)) || !(fEMEventCut.IsSelected(c2))) {

PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx

@@ -348,6 +348,7 @@ struct TaskPi0FlowEMC {
     fEMCCut.SetClusterizer(emccuts.clusterDefinition);
     fEMCCut.SetUseTM(emccuts.cfgEMCUseTM.value);                // disables or enables TM
     fEMCCut.SetUseSecondaryTM(emccuts.emcUseSecondaryTM.value); // disables or enables secondary TM
+    fEMCCut.SetDoQA(emccuts.cfgEnableQA.value);
   }
 
   void definePCMCut()
@@ -1095,7 +1096,7 @@ struct TaskPi0FlowEMC {
   void processEMCal(CollsWithQvecs const& collisions, EMCalPhotons const& clusters, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds)
   {
     EMBitFlags flags(clusters.size());
-    fEMCCut.AreSelectedRunning(flags, clusters, matchedPrims, matchedSeconds);
+    fEMCCut.AreSelectedRunning(flags, clusters, matchedPrims, matchedSeconds, &registry);
 
     if (cfgDoReverseScaling.value) {
       energyCorrectionFactor = 1.0505f;
@@ -1238,7 +1239,7 @@ struct TaskPi0FlowEMC {
   void processEMCalPCMC(CollsWithQvecs const& collisions, EMCalPhotons const& clusters, PCMPhotons const& photons, aod::V0Legs const&, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds)
   {
     EMBitFlags emcFlags(clusters.size());
-    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds);
+    fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, &registry);
 
     if (cfgDoReverseScaling.value) {
       energyCorrectionFactor = 1.0505f;