Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mhemmer-cern wants to merge 2 commits into
base: master
Choose a base branch
from
Open

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
53eaaa7
[PWGEM] PhotonMeson: Add basic EMCal QA to EMCPhotonCut
mhemmer-cern Feb 3, 2026
243024c
fix missing argument in calibETaskEmc.cxx
mhemmer-cern Feb 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions PWGEM/PhotonMeson/Core/EMCPhotonCut.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:";
Expand Down
151 changes: 145 additions & 6 deletions PWGEM/PhotonMeson/Core/EMCPhotonCut.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "PWGEM/PhotonMeson/DataModel/gammaTables.h"

#include <Framework/ASoA.h>
#include <Framework/HistogramRegistry.h>

#include <TNamed.h>

Expand Down Expand Up @@ -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,
Expand All @@ -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
Expand All @@ -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()) {
Expand All @@ -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;
}
}
Expand All @@ -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;
}

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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};
Expand Down
8 changes: 4 additions & 4 deletions PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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))) {
Expand Down
5 changes: 3 additions & 2 deletions PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
Loading