From dbe4476edc5da829f919c3e91d3c7c1a0b7d8770 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 13 Nov 2026 07:16:48 +0000 Subject: [PATCH] Please consider the following formatting changes --- PWGJE/Tasks/promptphotonQA.cxx | 1181 +++++++++++++++++++------------- 1 file changed, 689 insertions(+), 492 deletions(-) diff --git a/PWGJE/Tasks/promptphotonQA.cxx b/PWGJE/Tasks/promptphotonQA.cxx index 2992336b4c3..0f16e5f2b12 100644 --- a/PWGJE/Tasks/promptphotonQA.cxx +++ b/PWGJE/Tasks/promptphotonQA.cxx @@ -17,12 +17,12 @@ #include "PWGJE/DataModel/EMCALClusters.h" #include "PWGJE/DataModel/Jet.h" #include "PWGJE/DataModel/JetReducedData.h" -#include "Common/DataModel/EventSelection.h" + #include "Common/Core/RecoDecay.h" -#include "CommonConstants/MathConstants.h" -#include +#include "Common/DataModel/EventSelection.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "Framework/O2DatabasePDGPlugin.h" + +#include #include #include #include @@ -30,16 +30,20 @@ #include #include #include +#include #include #include + +#include + #include #include #include #include #include #include -#include #include +#include using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -47,9 +51,10 @@ using namespace o2::framework::expressions; #include #include #include + +#include #include #include -#include #include struct PromptphotonQA { @@ -72,8 +77,13 @@ struct PromptphotonQA { static constexpr double MaxExtrapolationAbsEta = 10.; static constexpr double MinMatchingExtrapolationPhi = -10.; - //basic enum for cateogrization of prompt photons - enum class Origin : int { Unclassified = 0, Direct = 1, Fragmentation = 2, Decay = 3, Excluded = 4, Mixed = 5 }; + // basic enum for cateogrization of prompt photons + enum class Origin : int { Unclassified = 0, + Direct = 1, + Fragmentation = 2, + Decay = 3, + Excluded = 4, + Mixed = 5 }; struct OriginalResult { Origin origin = Origin::Unclassified; bool leptonBranch = false; @@ -144,8 +154,7 @@ struct PromptphotonQA { (aod::jcluster::nlm <= cfgMaxNLM) && (aod::jcluster::isExotic == cfgExoticContribution); - - //AOD merged tables. We of course require the correct task for all of these tables. + // AOD merged tables. We of course require the correct task for all of these tables. using TrackCandidates = soa::Join; using RecoTracks = soa::Join; using DataCollisions = soa::Join; @@ -155,61 +164,80 @@ struct PromptphotonQA { using FilteredDataClusters = soa::Filtered; using FilteredMCClusters = soa::Filtered; using BcCandidates = soa::Join; - //init presciles and filters + // init presciles and filters PresliceUnsorted emcTrackPerTrack = aod::jemctrack::trackId; Preslice clustersPerCollision = aod::jcluster::collisionId; Preslice tracksPerCollision = aod::jtrack::collisionId; int trackFilter = -1; std::vector triggerMaskBits; - - //WRAPPERS to fill histograms WITH and WITHOUT the pThat weight from the MC. QA1 is for 1D histos, QA2 is for 2d histos. - //by having these wrappers we ensure that the weighted and unweighted histos are always filed simultaneously and identically. + // WRAPPERS to fill histograms WITH and WITHOUT the pThat weight from the MC. QA1 is for 1D histos, QA2 is for 2d histos. + // by having these wrappers we ensure that the weighted and unweighted histos are always filed simultaneously and identically. struct QA1 { std::shared_ptr weighted, unweighted; - void fill(double x, double weight) { - weighted->Fill(x,weight); - if (unweighted) unweighted->Fill(x,1.0); + void fill(double x, double weight) + { + weighted->Fill(x, weight); + if (unweighted) + unweighted->Fill(x, 1.0); } - void setBinLabel(int bin, const char* label) { - weighted->GetXaxis()->SetBinLabel(bin,label); - if (unweighted) unweighted->GetXaxis()->SetBinLabel(bin,label); + void setBinLabel(int bin, const char* label) + { + weighted->GetXaxis()->SetBinLabel(bin, label); + if (unweighted) + unweighted->GetXaxis()->SetBinLabel(bin, label); } }; struct QA2 { std::shared_ptr weighted, unweighted; - void fill(double x, double y, double weight) { - weighted->Fill(x,y,weight); - if (unweighted) unweighted->Fill(x,y,1.0); + void fill(double x, double y, double weight) + { + weighted->Fill(x, y, weight); + if (unweighted) + unweighted->Fill(x, y, 1.0); } }; // Histogram storage and selection stages, as well as several helper structs - enum PerformanceStage { BeforeM02, AfterM02, Neutral, Isolated, IsolatedEOverP, - StrictNeutral, MatchRecovery, MatchRecoveryIsolated, nPerformanceStages }; + enum PerformanceStage { BeforeM02, + AfterM02, + Neutral, + Isolated, + IsolatedEOverP, + StrictNeutral, + MatchRecovery, + MatchRecoveryIsolated, + nPerformanceStages }; const std::array performanceStages{ "beforeM02", "afterM02", "neutral", "isolated", "isolatedEOverP", "strictNeutral", "matchedEOverPRecovery", "matchedEOverPRecoveryIsolated"}; // Exclusive truth categories for energy stacks. Inclusive prompt is an overlay, // since direct, fragmentation and mixed prompt are already separate stack entries. - enum class ClusterOrigin : int { Direct, Fragmentation, MixedPrompt, Decay, Charged, - Other, NoDominant, UnknownTruth, Count }; + enum class ClusterOrigin : int { Direct, + Fragmentation, + MixedPrompt, + Decay, + Charged, + Other, + NoDominant, + UnknownTruth, + Count }; static constexpr int NClusterOrigins = static_cast(ClusterOrigin::Count); - const std::array clusterOrigins{ + const std::array clusterOrigins{ "direct", "fragmentation", "mixedPrompt", "decay", "charged", "other", "noDominant", "unknownTruth"}; - const std::array analysisEnergyStages{ + const std::array analysisEnergyStages{ "beforeM02", "afterM02", "validMatching", "neutral", "coneEligible", "isolated", "isolatedEOverP"}; struct ClusterEnergyQA { std::shared_ptr total, prompt; - std::array,NClusterOrigins> exclusive; + std::array, NClusterOrigins> exclusive; }; struct AnalysisCutQA { std::shared_ptr candidates, goodSignal; std::shared_ptr candidatesVsPt, goodSignalVsPt; - std::array energy; + std::array energy; }; struct ConeParticle { @@ -241,10 +269,10 @@ struct PromptphotonQA { std::shared_ptr energy, pt, eta, phi, m02, m20, cells, time, nlm; std::shared_ptr m02Energy, m20Energy, etaPhi; }; - + struct QAState { - std::array,2> analysisCuts; - std::array,2> clusterEnergy; + std::array, 2> analysisCuts; + std::array, 2> clusterEnergy; // Neutral, cone-eligible neutral, and isolated neutral candidate samples. std::array signalValidation; std::array tracks; @@ -260,10 +288,10 @@ struct PromptphotonQA { // Efficiency sets are reset once per dataframe, not per reconstructed collision. // This also deduplicates photons across split reconstructed collisions. - std::array,nPerformanceStages>,2> efficiencyMatched; - std::array,2> efficiencyAll, efficiencyAccepted, coneAll, coneAccepted; - std::array,nPerformanceStages>,2> efficiencyNumerator, purityNumerator; - std::array,nPerformanceStages> purityDenominator, validLabels; + std::array, nPerformanceStages>, 2> efficiencyMatched; + std::array, 2> efficiencyAll, efficiencyAccepted, coneAll, coneAccepted; + std::array, nPerformanceStages>, 2> efficiencyNumerator, purityNumerator; + std::array, nPerformanceStages> purityDenominator, validLabels; bool recoEventAccepted = false; // Stages: quality+acceptance before M02; after M02; charged-veto pass; veto fail. @@ -277,88 +305,95 @@ struct PromptphotonQA { std::shared_ptr genPromptEta, genPromptPhi; std::array, 2> genPhi, genEta, genPt; std::shared_ptr genEvents, ancestryStatus, directBeforeCalo, directInCalo; - std::array,2> selectedEventTruth, selectedConeTruth; - std::array,nPerformanceStages> directFiducialTagged, fiducialCandidates; + std::array, 2> selectedEventTruth, selectedConeTruth; + std::array, nPerformanceStages> directFiducialTagged, fiducialCandidates; }; QAState qaState; //************************************************** - //HELPER STRUCTS DONE////////////////////////////// + // HELPER STRUCTS DONE////////////////////////////// //************************************************** - - //NOW FOR HELPER FUNCTIONS FOR HISTOGRAM BOOKING + // NOW FOR HELPER FUNCTIONS FOR HISTOGRAM BOOKING //////////////////////////////////////////////////// - //checks whether isolation is required + // checks whether isolation is required static bool isolatedStage(int stage) { return stage == Isolated || stage == IsolatedEOverP || stage == MatchRecoveryIsolated; } - //simple phi02pi wrapper - static double phi02pi(double phi) { + // simple phi02pi wrapper + static double phi02pi(double phi) + { return RecoDecay::constrainAngle(phi); } - //Checks whether origin is prompt for a cluster - static bool isPrompt(Origin origin) { + // Checks whether origin is prompt for a cluster + static bool isPrompt(Origin origin) + { return origin == Origin::Direct || origin == Origin::Fragmentation || origin == Origin::Mixed; } - + // Each Fill makes the weighted and unit-weight fill together, with identical selection. - static std::string unweightedTitle(std::string title) { + static std::string unweightedTitle(std::string title) + { for (auto const& word : {std::string("Weighted "), std::string("weighted ")}) { std::size_t pos; - while ((pos = title.find(word)) != std::string::npos) title.erase(pos, word.size()); + while ((pos = title.find(word)) != std::string::npos) + title.erase(pos, word.size()); } auto pos = title.find("Sum weights"); - if (pos != std::string::npos) title.replace(pos, 11, "Events"); + if (pos != std::string::npos) + title.replace(pos, 11, "Events"); return title; } - //folder bookkeeping in the root file - static bool isMCPath(std::string const& name) { - return name.rfind("REC/",0)==0 || name.rfind("GEN/",0)==0; + // folder bookkeeping in the root file + static bool isMCPath(std::string const& name) + { + return name.rfind("REC/", 0) == 0 || name.rfind("GEN/", 0) == 0; } - static std::string weightPath(std::string const& name, bool weighted) { + static std::string weightPath(std::string const& name, bool weighted) + { const auto slash = name.find('/'); - return name.substr(0,slash+1)+(weighted ? "weighted/" : "unweighted/")+name.substr(slash+1); + return name.substr(0, slash + 1) + (weighted ? "weighted/" : "unweighted/") + name.substr(slash + 1); } - - //NOW WE BOOKKEEP ALL OF THE HISTOGRAMS - std::shared_ptr h1(std::string const& name, std::string const& title, int n, double lo, double hi, bool extended = false) { + // NOW WE BOOKKEEP ALL OF THE HISTOGRAMS + std::shared_ptr h1(std::string const& name, std::string const& title, int n, double lo, double hi, bool extended = false) + { auto& registry = extended ? extendedHistos : histos; auto result = std::make_shared(); const bool paired = isMCPath(name); const auto mainTitle = paired ? title : unweightedTitle(title); - result->weighted = registry.add(paired ? weightPath(name,true) : name, mainTitle.c_str(), kTH1D, {{n,lo,hi}}, true); - if (paired) result->unweighted = registry.add(weightPath(name,false), unweightedTitle(title).c_str(), kTH1D, {{n,lo,hi}}, true); + result->weighted = registry.add(paired ? weightPath(name, true) : name, mainTitle.c_str(), kTH1D, {{n, lo, hi}}, true); + if (paired) + result->unweighted = registry.add(weightPath(name, false), unweightedTitle(title).c_str(), kTH1D, {{n, lo, hi}}, true); return result; } - - std::shared_ptr h2(std::string const& name, std::string const& title, AxisSpec x, AxisSpec y, bool extended = false) { + + std::shared_ptr h2(std::string const& name, std::string const& title, AxisSpec x, AxisSpec y, bool extended = false) + { auto& registry = extended ? extendedHistos : histos; auto result = std::make_shared(); const bool paired = isMCPath(name); - result->weighted = registry.add(paired ? weightPath(name,true) : name, title.c_str(), kTH2D, {x,y}, true); - if (paired) result->unweighted = registry.add(weightPath(name,false), unweightedTitle(title).c_str(), kTH2D, {x,y}, true); + result->weighted = registry.add(paired ? weightPath(name, true) : name, title.c_str(), kTH2D, {x, y}, true); + if (paired) + result->unweighted = registry.add(weightPath(name, false), unweightedTitle(title).c_str(), kTH2D, {x, y}, true); return result; } - + std::shared_ptr bookMap(std::string const& path) { - return h2(path, ";#eta;#varphi (rad)", {100,-1,1}, {144,0,o2::constants::math::TwoPI}, true); + return h2(path, ";#eta;#varphi (rad)", {100, -1, 1}, {144, 0, o2::constants::math::TwoPI}, true); } IsolationQA bookIsolation(std::string const& path) { IsolationQA qa; - qa.sumPt = h2(path + "/coneSumPt", ";Candidate pT (GeV/c);Charged cone sum pT (GeV/c)", {250,0,500}, {250,0,50}, true); - qa.eOverP = h2(path + "/coneEOverP", ";Candidate pT (GeV/c);E / cone sum |p|", {250,0,500}, {200,0,10}, true); - const std::array*>,4> entries{{ - {"/eligible", &qa.eligible}, {"/isolated", &qa.isolated}, - {"/eOverPPass", &qa.eOverPPass}, {"/combined", &qa.combined}}}; + qa.sumPt = h2(path + "/coneSumPt", ";Candidate pT (GeV/c);Charged cone sum pT (GeV/c)", {250, 0, 500}, {250, 0, 50}, true); + qa.eOverP = h2(path + "/coneEOverP", ";Candidate pT (GeV/c);E / cone sum |p|", {250, 0, 500}, {200, 0, 10}, true); + const std::array*>, 4> entries{{{"/eligible", &qa.eligible}, {"/isolated", &qa.isolated}, {"/eOverPPass", &qa.eOverPPass}, {"/combined", &qa.combined}}}; for (auto const& entry : entries) { *entry.second = h1(path + entry.first, ";Candidate pT (GeV/c);Candidates", 250, 0, 500, true); } @@ -383,106 +418,121 @@ struct PromptphotonQA { q.etaPhi->fill(eta, phi02pi(phi), w); } - void bookCluster(ClusterQA& q, std::string const& p) { + void bookCluster(ClusterQA& q, std::string const& p) + { q.etaPhi = bookMap(p + "/etaPhi"); - q.energy = h1(p+"/energy", ";E_{cluster} (GeV);Weighted clusters", 250, 0, 500); - q.pt = h1(p+"/pt", ";E/cosh(eta) (GeV/c);Weighted clusters", 250, 0, 500); - q.eta = h1(p+"/eta", ";eta;Weighted clusters", 80, -0.8, 0.8); - q.phi = h1(p+"/phi", ";phi;Weighted clusters", 90, 0, o2::constants::math::TwoPI); - q.m02 = h1(p+"/m02", ";M02;Weighted clusters", 100, 0, 2); - q.m20 = h1(p+"/m20", ";M20;Weighted clusters", 100, 0, 2); - q.cells = h1(p+"/nCells", ";N cells;Weighted clusters", 100, -0.5, 99.5); - q.time = h1(p+"/time", ";Time (ns);Weighted clusters", 200, -100, 100); - q.nlm = h1(p+"/nlm", ";N local maxima;Weighted clusters", 20, -0.5, 19.5); - q.m02Energy = h2(p+"/m02VsEnergy", ";E (GeV);M02", {250,0,500}, {100,0,2}); - q.m20Energy = h2(p+"/m20VsEnergy", ";E (GeV);M20", {250,0,500}, {100,0,2}); + q.energy = h1(p + "/energy", ";E_{cluster} (GeV);Weighted clusters", 250, 0, 500); + q.pt = h1(p + "/pt", ";E/cosh(eta) (GeV/c);Weighted clusters", 250, 0, 500); + q.eta = h1(p + "/eta", ";eta;Weighted clusters", 80, -0.8, 0.8); + q.phi = h1(p + "/phi", ";phi;Weighted clusters", 90, 0, o2::constants::math::TwoPI); + q.m02 = h1(p + "/m02", ";M02;Weighted clusters", 100, 0, 2); + q.m20 = h1(p + "/m20", ";M20;Weighted clusters", 100, 0, 2); + q.cells = h1(p + "/nCells", ";N cells;Weighted clusters", 100, -0.5, 99.5); + q.time = h1(p + "/time", ";Time (ns);Weighted clusters", 200, -100, 100); + q.nlm = h1(p + "/nlm", ";N local maxima;Weighted clusters", 20, -0.5, 19.5); + q.m02Energy = h2(p + "/m02VsEnergy", ";E (GeV);M02", {250, 0, 500}, {100, 0, 2}); + q.m20Energy = h2(p + "/m20VsEnergy", ";E (GeV);M20", {250, 0, 500}, {100, 0, 2}); } - - template void fillCluster(ClusterQA& q, C const& c, double w) { + + template + void fillCluster(ClusterQA& q, C const& c, double w) + { q.etaPhi->fill(c.eta(), phi02pi(c.phi()), w); - q.energy->fill(c.energy(),w); q.pt->fill(c.energy()/std::cosh(c.eta()),w); - q.eta->fill(c.eta(),w); q.phi->fill(phi02pi(c.phi()),w); - q.m02->fill(c.m02(),w); q.m20->fill(c.m20(),w); - q.cells->fill(c.nCells(),w); q.time->fill(c.time(),w); q.nlm->fill(c.nlm(),w); - q.m02Energy->fill(c.energy(),c.m02(),w); q.m20Energy->fill(c.energy(),c.m20(),w); + q.energy->fill(c.energy(), w); + q.pt->fill(c.energy() / std::cosh(c.eta()), w); + q.eta->fill(c.eta(), w); + q.phi->fill(phi02pi(c.phi()), w); + q.m02->fill(c.m02(), w); + q.m20->fill(c.m20(), w); + q.cells->fill(c.nCells(), w); + q.time->fill(c.time(), w); + q.nlm->fill(c.nlm(), w); + q.m02Energy->fill(c.energy(), c.m02(), w); + q.m20Energy->fill(c.energy(), c.m20(), w); } //********************************************* - //DONE WITH HISTO BOOKING. NOW FOR THE MAIN INIT + // DONE WITH HISTO BOOKING. NOW FOR THE MAIN INIT //********************************************* std::shared_ptr energyH1(std::string const& name, std::string const& title) { - const bool paired=isMCPath(name); - auto& registry=paired ? clusterEnergyHistos : dataClusterEnergyHistos; - auto result=std::make_shared(); - result->weighted=registry.add(paired ? weightPath(name,true) : name, - title.c_str(),kTH1D,{{250,0,500}},true); - if (paired) result->unweighted=registry.add(weightPath(name,false), - unweightedTitle(title).c_str(),kTH1D,{{250,0,500}},true); + const bool paired = isMCPath(name); + auto& registry = paired ? clusterEnergyHistos : dataClusterEnergyHistos; + auto result = std::make_shared(); + result->weighted = registry.add(paired ? weightPath(name, true) : name, + title.c_str(), kTH1D, {{250, 0, 500}}, true); + if (paired) + result->unweighted = registry.add(weightPath(name, false), + unweightedTitle(title).c_str(), kTH1D, {{250, 0, 500}}, true); return result; } ClusterEnergyQA bookClusterEnergy(std::string const& path, bool truth) { ClusterEnergyQA q; - q.total=energyH1(path+"/total", ";Cluster E (GeV);Clusters"); + q.total = energyH1(path + "/total", ";Cluster E (GeV);Clusters"); if (truth) { - q.prompt=energyH1(path+"/prompt", ";Cluster E (GeV);Prompt-dominated clusters (inclusive overlay)"); - for (int i=0; i branches{"neutral", "strictNeutral", "matchedEOverPRecovery"}; - const std::array labels{ + const std::array branches{"neutral", "strictNeutral", "matchedEOverPRecovery"}; + const std::array labels{ "quality+energy+configured acceptance", "M02", "valid track matching", "neutral branch selection", "cone eligible", "charged isolation", "cone E/sum|p|"}; - for (int sample=0; sample<2; ++sample) for (int branch=0; branch<3; ++branch) { - auto& q=qaState.analysisCuts[sample][branch]; - const std::string base=(sample==0 ? "REC/analysisCuts/" : "DATA/analysisCuts/")+branches[branch]+"/"; - for (int stage=0; stage<7; ++stage) - q.energy[stage]=bookClusterEnergy(base+"energy/"+analysisEnergyStages[stage],sample==0); - q.candidates=h1(base+"candidates", ";Cumulative selection;Candidates",7,-0.5,6.5,true); - q.candidatesVsPt=h2(base+"candidatesVsPt", ";Cumulative selection;Cluster pT (GeV/c)", - {7,-0.5,6.5},{250,0,500},true); - if (sample==0) { - q.goodSignal=h1(base+"goodSignal", ";Cumulative selection;Prompt-dominated candidates",7,-0.5,6.5,true); - q.goodSignalVsPt=h2(base+"goodSignalVsPt", ";Cumulative selection;Cluster pT (GeV/c)", - {7,-0.5,6.5},{250,0,500},true); - } - for (int bin=1; bin<=7; ++bin) { - q.candidates->setBinLabel(bin,labels[bin-1]); - if (q.goodSignal) q.goodSignal->setBinLabel(bin,labels[bin-1]); - for (auto const& hist : {q.candidatesVsPt,q.goodSignalVsPt}) { - if (!hist) continue; - hist->weighted->GetXaxis()->SetBinLabel(bin,labels[bin-1]); - if (hist->unweighted) hist->unweighted->GetXaxis()->SetBinLabel(bin,labels[bin-1]); + for (int sample = 0; sample < 2; ++sample) + for (int branch = 0; branch < 3; ++branch) { + auto& q = qaState.analysisCuts[sample][branch]; + const std::string base = (sample == 0 ? "REC/analysisCuts/" : "DATA/analysisCuts/") + branches[branch] + "/"; + for (int stage = 0; stage < 7; ++stage) + q.energy[stage] = bookClusterEnergy(base + "energy/" + analysisEnergyStages[stage], sample == 0); + q.candidates = h1(base + "candidates", ";Cumulative selection;Candidates", 7, -0.5, 6.5, true); + q.candidatesVsPt = h2(base + "candidatesVsPt", ";Cumulative selection;Cluster pT (GeV/c)", + {7, -0.5, 6.5}, {250, 0, 500}, true); + if (sample == 0) { + q.goodSignal = h1(base + "goodSignal", ";Cumulative selection;Prompt-dominated candidates", 7, -0.5, 6.5, true); + q.goodSignalVsPt = h2(base + "goodSignalVsPt", ";Cumulative selection;Cluster pT (GeV/c)", + {7, -0.5, 6.5}, {250, 0, 500}, true); + } + for (int bin = 1; bin <= 7; ++bin) { + q.candidates->setBinLabel(bin, labels[bin - 1]); + if (q.goodSignal) + q.goodSignal->setBinLabel(bin, labels[bin - 1]); + for (auto const& hist : {q.candidatesVsPt, q.goodSignalVsPt}) { + if (!hist) + continue; + hist->weighted->GetXaxis()->SetBinLabel(bin, labels[bin - 1]); + if (hist->unweighted) + hist->unweighted->GetXaxis()->SetBinLabel(bin, labels[bin - 1]); + } } } - } } - void init(InitContext const&) { + void init(InitContext const&) + { bookAnalysisCuts(); - + if (doprocessMCGen && doprocessMCGenJE) throw std::runtime_error("Select only one GEN input: processMCGen (truth only) or processMCGenJE (JE truth and collision labels)"); if (doprocessMCGen && cfgGenReqRec) throw std::runtime_error("cfgGenReqRec requires JE labels: disable processMCGen and enable processMCGenJE, or run REC alone"); if (cfgMinPromptFraction <= 0.f || cfgMinPromptFraction > 1.f) throw std::runtime_error("cfgMinPromptFraction must be in (0,1]"); - if (cfgLowM02 > cfgHighM02 || cfgMinTime > cfgMaxTime || cfgLowClusterE > cfgHighClusterE) + if (cfgLowM02 > cfgHighM02 || cfgMinTime > cfgMaxTime || cfgLowClusterE > cfgHighClusterE) throw std::runtime_error("Invalid cluster selection bounds"); if (!(cfgIsoRadius > 0.f && cfgIsoMaxSumPt >= 0.f && cfgIsoMinEOverP >= 0.f)) throw std::runtime_error("Invalid isolation radius or threshold"); if (!(cfgMatchMinEOverP > 0.f && cfgDCGapEta >= 0.f && cfgDCGapEta < cfgEMChighPSeta && cfgDClowPSphi < cfgDCGapEndPhi && cfgDCGapEndPhi <= cfgDChighPSphi)) throw std::runtime_error("Invalid matched E/P or DCal geometry"); - - for (int sample = 0; sample < NSamples; ++sample) { + + for (int sample = 0; sample < NSamples; ++sample) { auto& q = qaState.tracks[sample]; const std::string base = sample == 0 ? "REC/tracking/" : "DATA/tracking/"; q.global = bookTrackKinematics(base + "global"); @@ -492,141 +542,144 @@ struct PromptphotonQA { q.surfaceStatus->setBinLabel(1, "missing/invalid extrapolation"); q.surfaceStatus->setBinLabel(2, "valid outside calo"); q.surfaceStatus->setBinLabel(3, "valid inside calo"); - q.etaShift = h2(base + "etaShift", ";Global pT;eta(surface)-eta(global)", {250,0,50}, {200,-1,1}, true); - q.phiShift = h2(base + "phiShift", ";Global pT;phi(surface)-phi(global)", {250,0,50}, {200,-o2::constants::math::PI,o2::constants::math::PI}, true); - + q.etaShift = h2(base + "etaShift", ";Global pT;eta(surface)-eta(global)", {250, 0, 50}, {200, -1, 1}, true); + q.phiShift = h2(base + "phiShift", ";Global pT;phi(surface)-phi(global)", {250, 0, 50}, {200, -o2::constants::math::PI, o2::constants::math::PI}, true); } - - const std::array trackTruthNames{ + + const std::array trackTruthNames{ "allGenerated", "acceptedEvents", "globalMatched", "surfaceExtrapolatedMatched", "surfaceInCaloMatched"}; - + for (int i = 0; i < NTrackingTruthSets; ++i) qaState.trackingTruth[i] = bookTrackKinematics("REC/tracking/performance/" + trackTruthNames[i]); - + qaState.trackingLabelStatus = h1("REC/tracking/labelStatus", ";Status;Selected tracks in accepted events", 5, -0.5, 4.5, true); const char* trackLabelNames[]{"missing label", "invalid label", "different MC collision", "not selected charged primary", "matched charged primary"}; - for (int i = 0; i < NTrackingTruthSets; ++i) qaState.trackingLabelStatus->setBinLabel(i+1, trackLabelNames[i]); - qaState.ancestryStatus = h1("GEN/ancestry", ";Classification;Photons in fiducial acceptance", 5,-0.5,4.5,true); + for (int i = 0; i < NTrackingTruthSets; ++i) + qaState.trackingLabelStatus->setBinLabel(i + 1, trackLabelNames[i]); + qaState.ancestryStatus = h1("GEN/ancestry", ";Classification;Photons in fiducial acceptance", 5, -0.5, 4.5, true); const char* ancestryLabels[]{"other", "direct", "fragmentation", "mixed", "lepton ancestry only or unresolved"}; - for (int i=0;i<5;++i) qaState.ancestryStatus->setBinLabel(i+1, ancestryLabels[i]); + for (int i = 0; i < 5; ++i) + qaState.ancestryStatus->setBinLabel(i + 1, ancestryLabels[i]); qaState.genIsolation = bookIsolation("GEN/isolation/prompt"); qaState.allPrimaryMap = bookMap("GEN/allPrimaryParticleEtaPhi"); qaState.allPhotonMap = bookMap("GEN/allPrimaryPhotonEtaPhi"); qaState.genParticleMap = bookMap("GEN/chargedPrimaryEtaPhi"); triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(cfgTriggerMasks); trackFilter = jetderiveddatautilities::initialiseTrackSelection(static_cast(cfgTrackFilter)); - const std::array stages{"beforeM02","afterM02","neutral","charged"}; - for (int sample=0; sample<2; ++sample) for (int stage=0; stage<4; ++stage) - qaState.clusterEnergy[sample][stage]=bookClusterEnergy( - (sample==0 ? "REC/clusterEnergy/" : "DATA/clusterEnergy/")+stages[stage],sample==0); - - for (int sample=0;sample<2;++sample) { + const std::array stages{"beforeM02", "afterM02", "neutral", "charged"}; + for (int sample = 0; sample < 2; ++sample) + for (int stage = 0; stage < 4; ++stage) + qaState.clusterEnergy[sample][stage] = bookClusterEnergy( + (sample == 0 ? "REC/clusterEnergy/" : "DATA/clusterEnergy/") + stages[stage], sample == 0); + + for (int sample = 0; sample < 2; ++sample) { const std::string p = sample == 0 ? "REC" : "DATA"; qaState.recoIsolation[sample] = bookIsolation(p + "/isolation/neutral"); qaState.allClusterMap[sample] = bookMap(p + "/allClusterEtaPhi"); qaState.trackMap[sample] = bookMap(p + "/selectedTrackEtaPhi"); - qaState.matchedEOverP[sample] = h2(p + "/association/eOverP", ";Track p (GeV/c);Cluster E / track p", {250,0,250}, {200,0,10}, true); - qaState.nLinked[sample] = h1(p + "/association/nLinked", ";Linked tracks;Clusters after M02", 32,-0.5,31.5, true); - qaState.nLinkedValid[sample] = h1(p+"/association/nLinkedValid", ";Linked tracks;Clusters with valid matching",32,-0.5,31.5,true); - qaState.invalidMatchingPt[sample] = h1(p+"/association/invalidMatchingPt", ";Cluster pT;Clusters with unusable matching",250,0,500,true); - qaState.recoveredIsoPt[sample] = h1(p+"/association/recoveryIsolatedPt", ";Cluster pT;Isolated unmatched or recovered clusters",250,0,500,true); - qaState.strictPt[sample] = h1(p+"/association/strictNeutralPt", ";Cluster pT;Unmatched clusters",250,0,500,true); - qaState.recoveryPt[sample] = h1(p+"/association/recoveryPt", ";Cluster pT;Unmatched or recovered clusters",250,0,500,true); - qaState.nMatched[sample] = h1(p + "/association/nMatched", ";Tracks passing matching window;Clusters with valid matching", 32,-0.5,31.5, true); - qaState.eventCuts[sample] = h1(p+"/eventCutFlow", ";Selection;Weighted events", 7, -0.5, 6.5); - const char* labels[]{"valid weight","vertex","sel8","EMC readout","trigger","event track","missing EMC extrapolation"}; - for(int b=0;b<7;++b) qaState.eventCuts[sample]->setBinLabel(b+1,labels[b]); - qaState.clusterCuts[sample] = h1(p+"/clusterCutFlow", ";Selection;Weighted clusters", 9, -0.5, 8.5); - const char* clabels[]{"all","definition","time","minimum energy","cells","NLM","exotic","energy+acceptance","M02"}; - for(int b=0;b<9;++b) qaState.clusterCuts[sample]->setBinLabel(b+1,clabels[b]); - qaState.sumPt[sample] = h1(p+"/matchedSumPt", ";Sum matched track pT (GeV/c);Weighted clusters",250,0,250); - qaState.eOverSumPt[sample] = h1(p+"/energyOverSumPt", ";E / sum pT;Weighted matched clusters",200,0,10); - qaState.matchResidual[sample] = h2(p+"/trackMatchResidual", ";Delta phi at EMCal;Delta eta at EMCal",{200,-0.1,0.1},{200,-0.1,0.1}); - for(int st=0;st<4;++st) bookCluster(qaState.blindQA[sample][st],p+"/"+stages[st]); + qaState.matchedEOverP[sample] = h2(p + "/association/eOverP", ";Track p (GeV/c);Cluster E / track p", {250, 0, 250}, {200, 0, 10}, true); + qaState.nLinked[sample] = h1(p + "/association/nLinked", ";Linked tracks;Clusters after M02", 32, -0.5, 31.5, true); + qaState.nLinkedValid[sample] = h1(p + "/association/nLinkedValid", ";Linked tracks;Clusters with valid matching", 32, -0.5, 31.5, true); + qaState.invalidMatchingPt[sample] = h1(p + "/association/invalidMatchingPt", ";Cluster pT;Clusters with unusable matching", 250, 0, 500, true); + qaState.recoveredIsoPt[sample] = h1(p + "/association/recoveryIsolatedPt", ";Cluster pT;Isolated unmatched or recovered clusters", 250, 0, 500, true); + qaState.strictPt[sample] = h1(p + "/association/strictNeutralPt", ";Cluster pT;Unmatched clusters", 250, 0, 500, true); + qaState.recoveryPt[sample] = h1(p + "/association/recoveryPt", ";Cluster pT;Unmatched or recovered clusters", 250, 0, 500, true); + qaState.nMatched[sample] = h1(p + "/association/nMatched", ";Tracks passing matching window;Clusters with valid matching", 32, -0.5, 31.5, true); + qaState.eventCuts[sample] = h1(p + "/eventCutFlow", ";Selection;Weighted events", 7, -0.5, 6.5); + const char* labels[]{"valid weight", "vertex", "sel8", "EMC readout", "trigger", "event track", "missing EMC extrapolation"}; + for (int b = 0; b < 7; ++b) + qaState.eventCuts[sample]->setBinLabel(b + 1, labels[b]); + qaState.clusterCuts[sample] = h1(p + "/clusterCutFlow", ";Selection;Weighted clusters", 9, -0.5, 8.5); + const char* clabels[]{"all", "definition", "time", "minimum energy", "cells", "NLM", "exotic", "energy+acceptance", "M02"}; + for (int b = 0; b < 9; ++b) + qaState.clusterCuts[sample]->setBinLabel(b + 1, clabels[b]); + qaState.sumPt[sample] = h1(p + "/matchedSumPt", ";Sum matched track pT (GeV/c);Weighted clusters", 250, 0, 250); + qaState.eOverSumPt[sample] = h1(p + "/energyOverSumPt", ";E / sum pT;Weighted matched clusters", 200, 0, 10); + qaState.matchResidual[sample] = h2(p + "/trackMatchResidual", ";Delta phi at EMCal;Delta eta at EMCal", {200, -0.1, 0.1}, {200, -0.1, 0.1}); + for (int st = 0; st < 4; ++st) + bookCluster(qaState.blindQA[sample][st], p + "/" + stages[st]); } - - const std::array truthNames{"prompt","direct","fragmentation"}; - for(int g=0;g<3;++g) for(int st=0;st<4;++st) - bookCluster(qaState.truthQA[g][st],"REC/truth_"+truthNames[g]+"/"+stages[st]); - - qaState.truthStatus = h1("REC/truthStatus", ";Status;Weighted clusters before M02",4,-0.5,3.5); - const char* tlabels[]{"valid labels","no labels","size mismatch","invalid label/fraction"}; - for(int b=0;b<4;++b) qaState.truthStatus->setBinLabel(b+1,tlabels[b]); - qaState.maxPromptFraction = h1("REC/maxPromptFraction", ";Largest single prompt fraction;Weighted clusters",101,0,1.01); - qaState.totalPromptFraction = h1("REC/totalPromptFraction", ";Sum prompt energy fractions;Weighted clusters",101,0,1.01); - qaState.promptFractionEnergy = h2("REC/maxPromptFractionVsEnergy", ";Cluster E (GeV);Largest prompt fraction",{250,0,500},{101,0,1.01}); - const std::array performanceGroups{"prompt","direct"}; - - for (int st=0;st= Isolated); - - for (int g=0;g<2;++g) { - const auto base = "REC/performance/"+performanceGroups[g]; - qaState.coneAll[g] = h1(base+"/efficiency/coneEligibleGenerated", ";True photon pT;Cone-eligible photons",250,0,500,true); - qaState.coneAccepted[g] = h1(base+"/efficiency/coneEligibleAcceptedEvents", ";True photon pT;Cone-eligible photons in accepted events",250,0,500,true); - qaState.selectedEventTruth[g] = h1(base+"/efficiency/selectedEventPopulation", ";True photon pT;Fiducial photons in configured GEN event population",250,0,500,true); - qaState.selectedConeTruth[g] = h1(base+"/efficiency/selectedConePopulation", ";True photon pT;Cone-eligible photons in configured GEN event population",250,0,500,true); - qaState.efficiencyAll[g] = h1(base+"/efficiency/allGenerated", ";True photon pT (GeV/c);Weighted photons",250,0,500); - qaState.efficiencyAccepted[g] = h1(base+"/efficiency/acceptedEvents", ";True photon pT (GeV/c);Weighted photons",250,0,500); - for (int st=0;st= Isolated); - qaState.purityNumerator[g][st] = h1(base+"/purity/"+performanceStages[st], ";Reconstructed cluster pT (GeV/c);Weighted tagged clusters",250,0,500, st >= Isolated); + + const std::array truthNames{"prompt", "direct", "fragmentation"}; + for (int g = 0; g < 3; ++g) + for (int st = 0; st < 4; ++st) + bookCluster(qaState.truthQA[g][st], "REC/truth_" + truthNames[g] + "/" + stages[st]); + + qaState.truthStatus = h1("REC/truthStatus", ";Status;Weighted clusters before M02", 4, -0.5, 3.5); + const char* tlabels[]{"valid labels", "no labels", "size mismatch", "invalid label/fraction"}; + for (int b = 0; b < 4; ++b) + qaState.truthStatus->setBinLabel(b + 1, tlabels[b]); + qaState.maxPromptFraction = h1("REC/maxPromptFraction", ";Largest single prompt fraction;Weighted clusters", 101, 0, 1.01); + qaState.totalPromptFraction = h1("REC/totalPromptFraction", ";Sum prompt energy fractions;Weighted clusters", 101, 0, 1.01); + qaState.promptFractionEnergy = h2("REC/maxPromptFractionVsEnergy", ";Cluster E (GeV);Largest prompt fraction", {250, 0, 500}, {101, 0, 1.01}); + const std::array performanceGroups{"prompt", "direct"}; + + for (int st = 0; st < nPerformanceStages; ++st) + qaState.validLabels[st] = h1("REC/performance/purity/" + performanceStages[st] + "/validLabels", + ";Cluster pT;Candidates with valid contributor labels", 250, 0, 500, true); + + for (int st = 0; st < nPerformanceStages; ++st) + qaState.purityDenominator[st] = h1("REC/performance/purity/" + performanceStages[st] + "/allCandidates", + ";Reconstructed cluster pT (GeV/c);Weighted clusters", 250, 0, 500, st >= Isolated); + + for (int g = 0; g < 2; ++g) { + const auto base = "REC/performance/" + performanceGroups[g]; + qaState.coneAll[g] = h1(base + "/efficiency/coneEligibleGenerated", ";True photon pT;Cone-eligible photons", 250, 0, 500, true); + qaState.coneAccepted[g] = h1(base + "/efficiency/coneEligibleAcceptedEvents", ";True photon pT;Cone-eligible photons in accepted events", 250, 0, 500, true); + qaState.selectedEventTruth[g] = h1(base + "/efficiency/selectedEventPopulation", ";True photon pT;Fiducial photons in configured GEN event population", 250, 0, 500, true); + qaState.selectedConeTruth[g] = h1(base + "/efficiency/selectedConePopulation", ";True photon pT;Cone-eligible photons in configured GEN event population", 250, 0, 500, true); + qaState.efficiencyAll[g] = h1(base + "/efficiency/allGenerated", ";True photon pT (GeV/c);Weighted photons", 250, 0, 500); + qaState.efficiencyAccepted[g] = h1(base + "/efficiency/acceptedEvents", ";True photon pT (GeV/c);Weighted photons", 250, 0, 500); + for (int st = 0; st < nPerformanceStages; ++st) { + qaState.efficiencyNumerator[g][st] = h1(base + "/efficiency/" + performanceStages[st], ";True photon pT (GeV/c);Weighted matched photons", 250, 0, 500, st >= Isolated); + qaState.purityNumerator[g][st] = h1(base + "/purity/" + performanceStages[st], ";Reconstructed cluster pT (GeV/c);Weighted tagged clusters", 250, 0, 500, st >= Isolated); } } // Use one event ensemble on both sides of the acceptance ratio. const std::string acceptance = "REC/performance/direct/acceptance/"; - qaState.directBeforeCalo = h1(acceptance+"beforeCalo", ";True photon pT;Direct photons before geometry",250,0,500,true); - qaState.directInCalo = h1(acceptance+"inCalo", ";True photon pT;Direct photons within geometry",250,0,500,true); - for (int st=0; st signalSamples{"neutral", "neutralConeEligible", "isolated"}; for (int i = 0; i < NSignalSamples; ++i) { auto& q = qaState.signalValidation[i]; const std::string base = "REC/signalValidation/" + signalSamples[i] + "/"; - const std::array*>,8> entries{{ - {"allCandidates", &q.allCandidates}, {"validTruth", &q.validTruth}, - {"anyPrompt", &q.anyPrompt}, {"goodSignal", &q.goodSignal}, - {"badSignal", &q.badSignal}, {"noPrompt", &q.noPrompt}, - {"promptBelowThreshold", &q.promptBelowThreshold}, {"unknownTruth", &q.unknownTruth}}}; + const std::array*>, 8> entries{{{"allCandidates", &q.allCandidates}, {"validTruth", &q.validTruth}, {"anyPrompt", &q.anyPrompt}, {"goodSignal", &q.goodSignal}, {"badSignal", &q.badSignal}, {"noPrompt", &q.noPrompt}, {"promptBelowThreshold", &q.promptBelowThreshold}, {"unknownTruth", &q.unknownTruth}}}; for (auto const& entry : entries) { *entry.second = h1(base + entry.first, ";Cluster pT (GeV/c);Candidates", 250, 0, 500, true); } q.maxPromptFraction = h2(base + "maxPromptFractionVsPt", - ";Cluster pT (GeV/c);Largest single prompt deposited-energy fraction", - {250, 0, 500}, {101, 0, 1.01}, true); + ";Cluster pT (GeV/c);Largest single prompt deposited-energy fraction", + {250, 0, 500}, {101, 0, 1.01}, true); q.totalPromptFraction = h2(base + "totalPromptFractionVsPt", - ";Cluster pT (GeV/c);Sum of prompt deposited-energy fractions", - {250, 0, 500}, {101, 0, 1.01}, true); + ";Cluster pT (GeV/c);Sum of prompt deposited-energy fractions", + {250, 0, 500}, {101, 0, 1.01}, true); } // Match the reconstructed truth-tagged cluster angular binning. qaState.genPromptEta = h1("GEN/prompt/eta", ";eta;Prompt photons", 80, -0.8, 0.8); qaState.genPromptPhi = h1("GEN/prompt/phi", ";phi;Prompt photons", 90, 0, o2::constants::math::TwoPI); - qaState.genEvents = h1("GEN/nEvents", ";MC events;Sum weights",1,0,1); - for(int g=0;g<2;++g) { - std::string p = g==0 ? "GEN/direct" : "GEN/frag"; + qaState.genEvents = h1("GEN/nEvents", ";MC events;Sum weights", 1, 0, 1); + for (int g = 0; g < 2; ++g) { + std::string p = g == 0 ? "GEN/direct" : "GEN/frag"; qaState.genPhotonMap[g] = bookMap(p + "_etaPhiQA"); - qaState.genPhi[g]=h1(p+"_phiQA",";phi;Weighted photons",320,0,o2::constants::math::TwoPI); - qaState.genEta[g]=h1(p+"_etaQA",";eta;Weighted photons",50,-1,1); - qaState.genPt[g]=h1(p+"_ptQA",";pT (GeV/c);Weighted photons",250,0,500); + qaState.genPhi[g] = h1(p + "_phiQA", ";phi;Weighted photons", 320, 0, o2::constants::math::TwoPI); + qaState.genEta[g] = h1(p + "_etaQA", ";eta;Weighted photons", 50, -1, 1); + qaState.genPt[g] = h1(p + "_ptQA", ";pT (GeV/c);Weighted photons", 250, 0, 500); } } - //********************************************* - //INIT COMPLETE - //NOW TIME FOR ANALYSIS-LEVEL HELPER FUNCTIONS //********************************************* - + // INIT COMPLETE + // NOW TIME FOR ANALYSIS-LEVEL HELPER FUNCTIONS + //********************************************* + // Truth classification: follow photon copies, then only quark/gluon ancestors // of a parton emission. Stop each branch at its nearest hard endpoint. - template + template OriginalResult classifyOriginal(Particle const& photon, Table const& particles) const { OriginalResult result; @@ -644,7 +697,8 @@ struct PromptphotonQA { pending.pop_back(); const auto local = id - static_cast(particles.offset()); if (local < 0 || local >= static_cast(particles.size()) || - !visited.emplace(id, partonBranch).second) continue; + !visited.emplace(id, partonBranch).second) + continue; auto mother = particles.rawIteratorAt(local); const int code = std::abs(mother.pdgCode()); const bool parton = (code >= kDown && code <= kTop) || code == kGluon; @@ -667,7 +721,8 @@ struct PromptphotonQA { } // Track, calorimeter and cone selections - template static bool finiteTrack(T const& t) + template + static bool finiteTrack(T const& t) { return std::isfinite(t.pt()) && std::isfinite(t.eta()) && std::isfinite(t.phi()); } @@ -677,7 +732,8 @@ struct PromptphotonQA { { if (!finiteTrack(track) || !std::isfinite(track.dcaXY()) || !std::isfinite(track.dcaZ()) || !std::isfinite(track.tpcCrossedRowsOverFindableCls()) || !std::isfinite(track.tpcChi2NCl()) || - !std::isfinite(track.itsChi2NCl())) return false; + !std::isfinite(track.itsChi2NCl())) + return false; // Basic track cuts if (track.pt() < cfgtrkMinPt) return false; @@ -715,14 +771,18 @@ struct PromptphotonQA { return true; }; // end of track selection - template bool selectedTrack(T const& track) { - if (!finiteTrack(track) || track.pt() < cfgtrkMinPt || std::abs(track.eta()) > cfgtrkMaxEta) return false; - if (cfgJETracks) return jetderiveddatautilities::selectTrack(track, trackFilter); + template + bool selectedTrack(T const& track) + { + if (!finiteTrack(track) || track.pt() < cfgtrkMinPt || std::abs(track.eta()) > cfgtrkMaxEta) + return false; + if (cfgJETracks) + return jetderiveddatautilities::selectTrack(track, trackFilter); auto original = track.template track_as(); return trackSelection(original) && original.isGlobalTrack() && (!cfgGlobalWoDCATrack || original.isGlobalTrackWoDCA()); } - + bool coneEligible(double eta) { return std::abs(eta) <= cfgtrkMaxEta && @@ -731,9 +791,11 @@ struct PromptphotonQA { bool inCalo(double eta, double phi) { - if (!(std::abs(eta) < cfgEMChighPSeta)) return false; + if (!(std::abs(eta) < cfgEMChighPSeta)) + return false; phi = phi02pi(phi); - if (phi > cfgEMClowPSphi && phi < cfgEMChighPSphi) return true; + if (phi > cfgEMClowPSphi && phi < cfgEMChighPSphi) + return true; return phi > cfgDClowPSphi && phi < cfgDChighPSphi && (phi >= cfgDCGapEndPhi || std::abs(eta) > cfgDCGapEta); } @@ -744,12 +806,14 @@ struct PromptphotonQA { std::pair measureIsolation(IsolationQA* qa, double pt, double eta, double phi, double energy, std::vector const& particles, double w) { - if (!coneEligible(eta)) return {false, false}; + if (!coneEligible(eta)) + return {false, false}; double sumPt = 0., sumP = 0.; for (auto const& particle : particles) { const double dEta = particle.eta - eta; const double dPhi = TVector2::Phi_mpi_pi(particle.phi - phi); - if (dEta*dEta + dPhi*dPhi >= cfgIsoRadius*cfgIsoRadius) continue; + if (dEta * dEta + dPhi * dPhi >= cfgIsoRadius * cfgIsoRadius) + continue; sumPt += particle.pt; sumP += particle.pt * std::cosh(particle.eta); } @@ -758,51 +822,61 @@ struct PromptphotonQA { if (qa) { qa->eligible->fill(pt, w); qa->sumPt->fill(pt, sumPt, w); - if (sumP > 0.) qa->eOverP->fill(pt, energy / sumP, w); - if (isolated) qa->isolated->fill(pt, w); - if (passesRatio) qa->eOverPPass->fill(pt, w); - if (isolated && passesRatio) qa->combined->fill(pt, w); + if (sumP > 0.) + qa->eOverP->fill(pt, energy / sumP, w); + if (isolated) + qa->isolated->fill(pt, w); + if (passesRatio) + qa->eOverPPass->fill(pt, w); + if (isolated && passesRatio) + qa->combined->fill(pt, w); } return {isolated, passesRatio}; } // Cluster truth tagging and performance counting - template + template bool isDecayPhoton(Particle const& photon, aod::JMcParticles const& particles) { - const int64_t first=static_cast(particles.offset()); - const int64_t end=first+static_cast(particles.size()); + const int64_t first = static_cast(particles.offset()); + const int64_t end = first + static_cast(particles.size()); std::vector pending; std::set visited; for (auto const& mother : photon.template mothers_as()) pending.push_back(mother.globalIndex()); while (!pending.empty()) { - const int64_t id=pending.back(); + const int64_t id = pending.back(); pending.pop_back(); - if (id=end || !visited.insert(id).second) continue; - const auto mother=particles.rawIteratorAt(id-first); - const int code=std::abs(mother.pdgCode()); - if (code>100 && code<1000000000) return true; - if (code!=kGamma) continue; + if (id < first || id >= end || !visited.insert(id).second) + continue; + const auto mother = particles.rawIteratorAt(id - first); + const int code = std::abs(mother.pdgCode()); + if (code > 100 && code < 1000000000) + return true; + if (code != kGamma) + continue; for (auto const& ancestor : mother.template mothers_as()) pending.push_back(ancestor.globalIndex()); } return false; } - template + template ClusterOrigin contributorOrigin(Particle const& particle, aod::JMcParticles const& particles) { - if (particle.pdgCode()==kGamma) { - const auto result=classifyOriginal(particle,particles); - if (result.origin==Origin::Direct) return ClusterOrigin::Direct; - if (result.origin==Origin::Fragmentation) return ClusterOrigin::Fragmentation; - if (result.origin==Origin::Mixed) return ClusterOrigin::MixedPrompt; - return isDecayPhoton(particle,particles) ? ClusterOrigin::Decay : ClusterOrigin::Other; + if (particle.pdgCode() == kGamma) { + const auto result = classifyOriginal(particle, particles); + if (result.origin == Origin::Direct) + return ClusterOrigin::Direct; + if (result.origin == Origin::Fragmentation) + return ClusterOrigin::Fragmentation; + if (result.origin == Origin::Mixed) + return ClusterOrigin::MixedPrompt; + return isDecayPhoton(particle, particles) ? ClusterOrigin::Decay : ClusterOrigin::Other; } // Classify the labelled contributor itself. Do not double-count its ancestors. - auto* species=pdg->GetParticle(particle.pdgCode()); - return species && species->Charge()!=0. ? ClusterOrigin::Charged : ClusterOrigin::Other; + auto* species = pdg->GetParticle(particle.pdgCode()); + return species && species->Charge() != 0. ? ClusterOrigin::Charged : ClusterOrigin::Other; } struct TruthTag { @@ -812,64 +886,80 @@ struct PromptphotonQA { bool anyPrompt = false; double maxPromptFraction = 0., totalPromptFraction = 0.; bool prompt = false, direct = false, fragmentation = false; - std::array,2> efficiencyIds; + std::array, 2> efficiencyIds; }; - - template TruthTag tagCluster(C const& c, aod::JMcParticles const& particles, double w) { + + template + TruthTag tagCluster(C const& c, aod::JMcParticles const& particles, double w) + { TruthTag tag; const auto ids = c.mcParticlesIds(); const auto fractions = c.amplitudeA(); // EMCAL producer stores getEnergyFractions(), NOT energies. - if (ids.size() != fractions.size()) { qaState.truthStatus->fill(2,w); return tag; } - if (ids.empty()) { qaState.truthStatus->fill(1,w); return tag; } + if (ids.size() != fractions.size()) { + qaState.truthStatus->fill(2, w); + return tag; + } + if (ids.empty()) { + qaState.truthStatus->fill(1, w); + return tag; + } // Validate the complete association before assigning any tag. Never index invalid labels. const int64_t first = particles.offset(); const int64_t end = first + particles.size(); double fractionSum = 0.; - for (std::size_t i=0;i= end || !std::isfinite(fractions[i]) || fractions[i] < 0.f || fractions[i] > MaxContributorFraction) { - qaState.truthStatus->fill(3,w); return tag; + qaState.truthStatus->fill(3, w); + return tag; } fractionSum += fractions[i]; } - if (fractionSum <= 0. || fractionSum > MaxContributorFractionSum) { qaState.truthStatus->fill(3,w); return tag; } + if (fractionSum <= 0. || fractionSum > MaxContributorFractionSum) { + qaState.truthStatus->fill(3, w); + return tag; + } tag.labelsValid = true; tag.stackOrigin = ClusterOrigin::NoDominant; double leadingFraction = -1.; int64_t leadingId = std::numeric_limits::max(); - qaState.truthStatus->fill(0.0,w); + qaState.truthStatus->fill(0.0, w); double maximum = 0., total = 0.; std::size_t i = 0; for (auto const& particle : c.template mcParticles_as()) { const double fraction = fractions[i++]; // The largest qualifying contributor gets exactly one stack category. // Below a 50% threshold several contributors can qualify; ties use the lowest ID. - if (fraction>cfgMinPromptFraction && - (fraction>leadingFraction || (fraction==leadingFraction && particle.globalIndex() cfgMinPromptFraction && + (fraction > leadingFraction || (fraction == leadingFraction && particle.globalIndex() < leadingId))) { + leadingFraction = fraction; + leadingId = particle.globalIndex(); + tag.stackOrigin = contributorOrigin(particle, particles); } - if (particle.pdgCode() != kGamma) continue; + if (particle.pdgCode() != kGamma) + continue; const auto result = classifyOriginal(particle, particles); - if (!isPrompt(result.origin)) continue; + if (!isPrompt(result.origin)) + continue; tag.anyPrompt |= fraction > 0.; maximum = std::max(maximum, fraction); total += fraction; - if (fraction <= cfgMinPromptFraction) continue; + if (fraction <= cfgMinPromptFraction) + continue; tag.efficiencyIds[0].push_back(particle.globalIndex()); - if (!result.directHardIndices.empty()) tag.efficiencyIds[1].push_back(particle.globalIndex()); + if (!result.directHardIndices.empty()) + tag.efficiencyIds[1].push_back(particle.globalIndex()); tag.prompt = true; tag.direct |= !result.directHardIndices.empty(); tag.fragmentation |= !result.fragmentationHardIndices.empty(); } tag.maxPromptFraction = maximum; tag.totalPromptFraction = total; - qaState.maxPromptFraction->fill(maximum,w); - qaState.totalPromptFraction->fill(total,w); - qaState.promptFractionEnergy->fill(c.energy(),maximum,w); + qaState.maxPromptFraction->fill(maximum, w); + qaState.totalPromptFraction->fill(total, w); + qaState.promptFractionEnergy->fill(c.energy(), maximum, w); return tag; } - + void fillSignalValidation(int sample, double pt, double w, TruthTag const& tag) { auto& q = qaState.signalValidation[sample]; @@ -883,86 +973,109 @@ struct PromptphotonQA { q.validTruth->fill(pt, w); q.maxPromptFraction->fill(pt, tag.maxPromptFraction, w); q.totalPromptFraction->fill(pt, tag.totalPromptFraction, w); - if (tag.anyPrompt) q.anyPrompt->fill(pt, w); + if (tag.anyPrompt) + q.anyPrompt->fill(pt, w); if (tag.prompt) { q.goodSignal->fill(pt, w); } else { q.badSignal->fill(pt, w); - if (tag.anyPrompt) q.promptBelowThreshold->fill(pt, w); - else q.noPrompt->fill(pt, w); + if (tag.anyPrompt) + q.promptBelowThreshold->fill(pt, w); + else + q.noPrompt->fill(pt, w); } } void fillClusterEnergy(ClusterEnergyQA& q, double energy, double w, TruthTag const& tag) { - q.total->fill(energy,w); - if (!q.prompt) return; // Data has no truth decomposition. - const int category=static_cast(tag.stackOrigin); - q.exclusive[category]->fill(energy,w); - if (tag.stackOrigin==ClusterOrigin::Direct || tag.stackOrigin==ClusterOrigin::Fragmentation || - tag.stackOrigin==ClusterOrigin::MixedPrompt) q.prompt->fill(energy,w); + q.total->fill(energy, w); + if (!q.prompt) + return; // Data has no truth decomposition. + const int category = static_cast(tag.stackOrigin); + q.exclusive[category]->fill(energy, w); + if (tag.stackOrigin == ClusterOrigin::Direct || tag.stackOrigin == ClusterOrigin::Fragmentation || + tag.stackOrigin == ClusterOrigin::MixedPrompt) + q.prompt->fill(energy, w); } void fillAnalysisCut(int sample, int branch, int stage, double pt, double energy, double w, TruthTag const& tag) { - auto& q=qaState.analysisCuts[sample][branch]; - q.candidates->fill(stage,w); - q.candidatesVsPt->fill(stage,pt,w); - fillClusterEnergy(q.energy[stage],energy,w,tag); + auto& q = qaState.analysisCuts[sample][branch]; + q.candidates->fill(stage, w); + q.candidatesVsPt->fill(stage, pt, w); + fillClusterEnergy(q.energy[stage], energy, w, tag); if (tag.prompt && q.goodSignal) { - q.goodSignal->fill(stage,w); - q.goodSignalVsPt->fill(stage,pt,w); + q.goodSignal->fill(stage, w); + q.goodSignalVsPt->fill(stage, pt, w); } } void fillAnalysisCommonCut(int sample, int stage, double pt, double energy, double w, TruthTag const& tag) { - for (int branch=0; branch<3; ++branch) fillAnalysisCut(sample,branch,stage,pt,energy,w,tag); + for (int branch = 0; branch < 3; ++branch) + fillAnalysisCut(sample, branch, stage, pt, energy, w, tag); } void fillAnalysisFinalCuts(int sample, int branch, double pt, double energy, double w, TruthTag const& tag, bool neutral, bool eligible, bool isolated, bool passesRatio) { - if (!neutral) return; - fillAnalysisCut(sample,branch,3,pt,energy,w,tag); - if (!eligible) return; - fillAnalysisCut(sample,branch,4,pt,energy,w,tag); - if (!isolated) return; - fillAnalysisCut(sample,branch,5,pt,energy,w,tag); - if (passesRatio) fillAnalysisCut(sample,branch,6,pt,energy,w,tag); + if (!neutral) + return; + fillAnalysisCut(sample, branch, 3, pt, energy, w, tag); + if (!eligible) + return; + fillAnalysisCut(sample, branch, 4, pt, energy, w, tag); + if (!isolated) + return; + fillAnalysisCut(sample, branch, 5, pt, energy, w, tag); + if (passesRatio) + fillAnalysisCut(sample, branch, 6, pt, energy, w, tag); } - template void fillPerformance(int stage, C const& c, double w, TruthTag const& tag) { - const double pt = c.energy()/std::cosh(c.eta()); - qaState.purityDenominator[stage]->fill(pt,w); - if (stage == Neutral) { - fillSignalValidation(0, pt, w, tag); - if (coneEligible(c.eta())) fillSignalValidation(1, pt, w, tag); - } else if (stage == Isolated) { - fillSignalValidation(2, pt, w, tag); - } - // Only matches inside the reconstructed fiducial region enter efficiency. - // Inclusive purity and blind QA can still cover the larger REC sample. - if (inCalo(c.eta(), c.phi()) && (!isolatedStage(stage) || coneEligible(c.eta()))) { - for (int g=0;g<2;++g) - qaState.efficiencyMatched[g][stage].insert(tag.efficiencyIds[g].begin(),tag.efficiencyIds[g].end()); - qaState.fiducialCandidates[stage]->fill(pt,w); - if (isolatedStage(stage) ? tag.coneFiducialDirect : tag.fiducialDirect) - qaState.directFiducialTagged[stage]->fill(pt,w); - } - if (tag.labelsValid) qaState.validLabels[stage]->fill(pt,w); - if (tag.prompt) qaState.purityNumerator[0][stage]->fill(pt,w); - if (tag.direct) qaState.purityNumerator[1][stage]->fill(pt,w); + template + void fillPerformance(int stage, C const& c, double w, TruthTag const& tag) + { + const double pt = c.energy() / std::cosh(c.eta()); + qaState.purityDenominator[stage]->fill(pt, w); + if (stage == Neutral) { + fillSignalValidation(0, pt, w, tag); + if (coneEligible(c.eta())) + fillSignalValidation(1, pt, w, tag); + } else if (stage == Isolated) { + fillSignalValidation(2, pt, w, tag); + } + // Only matches inside the reconstructed fiducial region enter efficiency. + // Inclusive purity and blind QA can still cover the larger REC sample. + if (inCalo(c.eta(), c.phi()) && (!isolatedStage(stage) || coneEligible(c.eta()))) { + for (int g = 0; g < 2; ++g) + qaState.efficiencyMatched[g][stage].insert(tag.efficiencyIds[g].begin(), tag.efficiencyIds[g].end()); + qaState.fiducialCandidates[stage]->fill(pt, w); + if (isolatedStage(stage) ? tag.coneFiducialDirect : tag.fiducialDirect) + qaState.directFiducialTagged[stage]->fill(pt, w); + } + if (tag.labelsValid) + qaState.validLabels[stage]->fill(pt, w); + if (tag.prompt) + qaState.purityNumerator[0][stage]->fill(pt, w); + if (tag.direct) + qaState.purityNumerator[1][stage]->fill(pt, w); } - template void fillStage(int sample, int stage, C const& c, double w, TruthTag const& tag) { - fillCluster(qaState.blindQA[sample][stage],c,w); - fillClusterEnergy(qaState.clusterEnergy[sample][stage],c.energy(),w,tag); - if (sample != 0) return; - if (stage <= Neutral) fillPerformance(stage, c, w, tag); - if (tag.prompt) fillCluster(qaState.truthQA[0][stage],c,w); - if (tag.direct) fillCluster(qaState.truthQA[1][stage],c,w); - if (tag.fragmentation) fillCluster(qaState.truthQA[2][stage],c,w); + template + void fillStage(int sample, int stage, C const& c, double w, TruthTag const& tag) + { + fillCluster(qaState.blindQA[sample][stage], c, w); + fillClusterEnergy(qaState.clusterEnergy[sample][stage], c.energy(), w, tag); + if (sample != 0) + return; + if (stage <= Neutral) + fillPerformance(stage, c, w, tag); + if (tag.prompt) + fillCluster(qaState.truthQA[0][stage], c, w); + if (tag.direct) + fillCluster(qaState.truthQA[1][stage], c, w); + if (tag.fragmentation) + fillCluster(qaState.truthQA[2][stage], c, w); } // The first six cut-flow bins include clusters rejected by the framework @@ -970,29 +1083,35 @@ struct PromptphotonQA { void fillBaselineCutFlow(aod::JClusters const& rawClusters, int sample, double w) { for (auto const& c : rawClusters) { - qaState.clusterCuts[sample]->fill(0.0,w); - if (c.definition() != cfgClusterDefinition) continue; + qaState.clusterCuts[sample]->fill(0.0, w); + if (c.definition() != cfgClusterDefinition) + continue; if (std::isfinite(c.eta()) && std::isfinite(c.phi())) qaState.allClusterMap[sample]->fill(c.eta(), phi02pi(c.phi()), w); - qaState.clusterCuts[sample]->fill(1,w); - if (!std::isfinite(c.time()) || c.time() < cfgMinTime || c.time() > cfgMaxTime) continue; - qaState.clusterCuts[sample]->fill(2,w); - if (!std::isfinite(c.energy()) || c.energy() <= cfgMinClusterEnergy) continue; - qaState.clusterCuts[sample]->fill(3,w); - if(c.nCells() < cfgMinNCells) continue; - qaState.clusterCuts[sample]->fill(4,w); - if(c.nlm() > cfgMaxNLM) continue; - qaState.clusterCuts[sample]->fill(5,w); + qaState.clusterCuts[sample]->fill(1, w); + if (!std::isfinite(c.time()) || c.time() < cfgMinTime || c.time() > cfgMaxTime) + continue; + qaState.clusterCuts[sample]->fill(2, w); + if (!std::isfinite(c.energy()) || c.energy() <= cfgMinClusterEnergy) + continue; + qaState.clusterCuts[sample]->fill(3, w); + if (c.nCells() < cfgMinNCells) + continue; + qaState.clusterCuts[sample]->fill(4, w); + if (c.nlm() > cfgMaxNLM) + continue; + qaState.clusterCuts[sample]->fill(5, w); } } // Use the same extrapolation source as the cluster-track matching branch. - template + template bool trackSurface(Track const& track, aod::JEMCTracks const& emctracks, double& eta, double& phi) { if (cfgJETracks) { auto rows = emctracks.sliceBy(emcTrackPerTrack, track.globalIndex()); - if (rows.size() == 0) return false; + if (rows.size() == 0) + return false; eta = rows.iteratorAt(0).etaEmcal(); phi = rows.iteratorAt(0).phiEmcal(); } else { @@ -1002,47 +1121,60 @@ struct PromptphotonQA { } // Reject extrapolation sentinels before wrapping phi. if (!std::isfinite(eta) || !std::isfinite(phi) || std::abs(eta) > MaxExtrapolationAbsEta || - phi < -o2::constants::math::PI || phi > o2::constants::math::TwoPI) return false; + phi < -o2::constants::math::PI || phi > o2::constants::math::TwoPI) + return false; phi = phi02pi(phi); return true; } - template + template bool trackingTruthSelected(Particle const& p) { if (!p.isPhysicalPrimary() || !finiteTrack(p) || p.pt() < cfgtrkMinPt || - std::abs(p.eta()) > cfgtrkMaxEta) return false; + std::abs(p.eta()) > cfgtrkMaxEta) + return false; auto* species = pdg->GetParticle(p.pdgCode()); return species && species->Charge() != 0.; } - template + template void collectTrackingMatches(Tracks const& tracks, aod::JMcTrackLbs const& labels, - aod::JMcParticles const& particles, int64_t collisionId, - aod::JEMCTracks const& emctracks, double w, - std::array,3>& matched) + aod::JMcParticles const& particles, int64_t collisionId, + aod::JEMCTracks const& emctracks, double w, + std::array, 3>& matched) { - const int64_t first=static_cast(particles.offset()); - const int64_t end=first+static_cast(particles.size()); + const int64_t first = static_cast(particles.offset()); + const int64_t end = first + static_cast(particles.size()); for (auto const& track : tracks) { - if (!selectedTrack(track)) continue; + if (!selectedTrack(track)) + continue; auto label = labels.rawIteratorAt(track.globalIndex() - labels.offset()); - if (!label.has_mcParticle()) { qaState.trackingLabelStatus->fill(0.,w); continue; } + if (!label.has_mcParticle()) { + qaState.trackingLabelStatus->fill(0., w); + continue; + } const int64_t id = label.mcParticleId(); if (id < first || id >= end) { - qaState.trackingLabelStatus->fill(1.,w); continue; + qaState.trackingLabelStatus->fill(1., w); + continue; } - auto truth = particles.rawIteratorAt(id-first); + auto truth = particles.rawIteratorAt(id - first); if (!truth.has_mcCollision() || truth.mcCollisionId() != collisionId) { - qaState.trackingLabelStatus->fill(2.,w); continue; + qaState.trackingLabelStatus->fill(2., w); + continue; + } + if (!trackingTruthSelected(truth)) { + qaState.trackingLabelStatus->fill(3., w); + continue; } - if (!trackingTruthSelected(truth)) { qaState.trackingLabelStatus->fill(3.,w); continue; } - qaState.trackingLabelStatus->fill(4.,w); + qaState.trackingLabelStatus->fill(4., w); matched[0].insert(id); double eta, phi; - if (!trackSurface(track, emctracks, eta, phi)) continue; + if (!trackSurface(track, emctracks, eta, phi)) + continue; matched[1].insert(id); - if (inCalo(eta, phi)) matched[2].insert(id); + if (inCalo(eta, phi)) + matched[2].insert(id); } } @@ -1051,7 +1183,7 @@ struct PromptphotonQA { double sumPt, sumP; }; - template + template ClusterMatch matchCluster(Cluster const& c, aod::JEMCTracks const& emctracks, int sample, double w, double clusterPt, bool& missingExtrapolation) { @@ -1061,71 +1193,96 @@ struct PromptphotonQA { qaState.nLinked[sample]->fill(linkedTracks.size(), w); int matchedCount = 0; for (auto const& track : linkedTracks) { - if (!finiteTrack(track)) { validMatching = false; continue; } - if(!selectedTrack(track)) continue; + if (!finiteTrack(track)) { + validMatching = false; + continue; + } + if (!selectedTrack(track)) + continue; double etaT, phiT; - if(cfgJETracks) { - auto extrapolated = emctracks.sliceBy(emcTrackPerTrack,track.globalIndex()); - if(extrapolated.size() == 0) { missingExtrapolation = true; validMatching = false; continue; } + if (cfgJETracks) { + auto extrapolated = emctracks.sliceBy(emcTrackPerTrack, track.globalIndex()); + if (extrapolated.size() == 0) { + missingExtrapolation = true; + validMatching = false; + continue; + } etaT = extrapolated.iteratorAt(0).etaEmcal(); phiT = extrapolated.iteratorAt(0).phiEmcal(); } else { auto original = track.template track_as(); - etaT = original.trackEtaEmcal(); phiT = original.trackPhiEmcal(); + etaT = original.trackEtaEmcal(); + phiT = original.trackPhiEmcal(); } - if(!std::isfinite(etaT) || !std::isfinite(phiT) || std::abs(etaT)>MaxExtrapolationAbsEta || phiT < MinMatchingExtrapolationPhi) { - missingExtrapolation = true; validMatching = false; continue; + if (!std::isfinite(etaT) || !std::isfinite(phiT) || std::abs(etaT) > MaxExtrapolationAbsEta || phiT < MinMatchingExtrapolationPhi) { + missingExtrapolation = true; + validMatching = false; + continue; } - const double dEta = etaT-c.eta(), dPhi = TVector2::Phi_mpi_pi(phiT-c.phi()); - qaState.matchResidual[sample]->fill(dPhi,dEta,w); - const double etaWindow = cfgPtClusterCut ? 0.010+std::pow(track.pt()+4.07,-2.5) : 0.05; - const double phiWindow = cfgPtClusterCut ? 0.015+std::pow(track.pt()+3.65,-2.0) : 0.05; - if(std::abs(dEta) < etaWindow && std::abs(dPhi) < phiWindow) { + const double dEta = etaT - c.eta(), dPhi = TVector2::Phi_mpi_pi(phiT - c.phi()); + qaState.matchResidual[sample]->fill(dPhi, dEta, w); + const double etaWindow = cfgPtClusterCut ? 0.010 + std::pow(track.pt() + 4.07, -2.5) : 0.05; + const double phiWindow = cfgPtClusterCut ? 0.015 + std::pow(track.pt() + 3.65, -2.0) : 0.05; + if (std::abs(dEta) < etaWindow && std::abs(dPhi) < phiWindow) { ++matchedCount; const double momentum = track.pt() * std::cosh(track.eta()); - if (momentum > 0.) qaState.matchedEOverP[sample]->fill(momentum, c.energy() / momentum, w); + if (momentum > 0.) + qaState.matchedEOverP[sample]->fill(momentum, c.energy() / momentum, w); matchedP += momentum; - matched = true; matchedPt += track.pt(); + matched = true; + matchedPt += track.pt(); } } - + // Incomplete matching cannot establish that a cluster is neutral. if (!validMatching) { - qaState.invalidMatchingPt[sample]->fill(clusterPt,w); + qaState.invalidMatchingPt[sample]->fill(clusterPt, w); return {false, matched, matchedPt, matchedP}; } - qaState.nLinkedValid[sample]->fill(linkedTracks.size(),w); + qaState.nLinkedValid[sample]->fill(linkedTracks.size(), w); qaState.nMatched[sample]->fill(matchedCount, w); return {true, matched, matchedPt, matchedP}; } // Shared reconstructed-event processing for MC and data. - template + template void processReco(Collision const& collision, Clusters const& clusters, aod::JClusters const& rawClusters, RecoTracks const& tracks, - aod::JEMCTracks const& emctracks, int sample, double w, Tagger tagger) { + aod::JEMCTracks const& emctracks, int sample, double w, Tagger tagger) + { qaState.recoEventAccepted = false; - if (!std::isfinite(w)) throw std::runtime_error("Non-finite MC event weight"); - qaState.eventCuts[sample]->fill(0.0,w); - if (!std::isfinite(collision.posZ()) || std::abs(collision.posZ()) >= cfgVtxCut) return; - qaState.eventCuts[sample]->fill(1,w); - if (!collision.sel8()) return; - qaState.eventCuts[sample]->fill(2,w); - if (cfgEmcTrigger && !collision.isEmcalReadout()) return; - qaState.eventCuts[sample]->fill(3,w); - if (!jetderiveddatautilities::selectTrigger(collision,triggerMaskBits)) return; - qaState.eventCuts[sample]->fill(4,w); + if (!std::isfinite(w)) + throw std::runtime_error("Non-finite MC event weight"); + qaState.eventCuts[sample]->fill(0.0, w); + if (!std::isfinite(collision.posZ()) || std::abs(collision.posZ()) >= cfgVtxCut) + return; + qaState.eventCuts[sample]->fill(1, w); + if (!collision.sel8()) + return; + qaState.eventCuts[sample]->fill(2, w); + if (cfgEmcTrigger && !collision.isEmcalReadout()) + return; + qaState.eventCuts[sample]->fill(3, w); + if (!jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) + return; + qaState.eventCuts[sample]->fill(4, w); if (cfgRequireEventTrack) { bool found = false; - for (auto const& track : tracks) if (selectedTrack(track)) { found = true; break; } - if (!found) return; + for (auto const& track : tracks) + if (selectedTrack(track)) { + found = true; + break; + } + if (!found) + return; } - qaState.eventCuts[sample]->fill(5,w); + qaState.eventCuts[sample]->fill(5, w); qaState.recoEventAccepted = true; bool missingExtrapolation = false; std::vector constituents; for (auto const& track : tracks) { - if (!selectedTrack(track)) continue; + if (!selectedTrack(track)) + continue; constituents.push_back({track.pt(), track.eta(), track.phi()}); qaState.trackMap[sample]->fill(track.eta(), phi02pi(track.phi()), w); auto& tq = qaState.tracks[sample]; @@ -1136,97 +1293,109 @@ struct PromptphotonQA { } else { // pT remains the global-track pT; only eta/phi are evaluated at the surface. fillTrackKinematics(tq.extrapolated, track.pt(), surfaceEta, surfacePhi, w); - tq.etaShift->fill(track.pt(), surfaceEta-track.eta(), w); - tq.phiShift->fill(track.pt(), TVector2::Phi_mpi_pi(surfacePhi-track.phi()), w); + tq.etaShift->fill(track.pt(), surfaceEta - track.eta(), w); + tq.phiShift->fill(track.pt(), TVector2::Phi_mpi_pi(surfacePhi - track.phi()), w); if (inCalo(surfaceEta, surfacePhi)) { fillTrackKinematics(tq.inCalo, track.pt(), surfaceEta, surfacePhi, w); tq.surfaceStatus->fill(2., w); } else { - tq.surfaceStatus->fill(1., w); - } + tq.surfaceStatus->fill(1., w); + } } } fillBaselineCutFlow(rawClusters, sample, w); for (auto const& c : clusters) { - qaState.clusterCuts[sample]->fill(6,w); - if(!std::isfinite(c.eta()) || !std::isfinite(c.phi()) || !std::isfinite(c.m02()) || !std::isfinite(c.m20())) continue; - if(c.energy() < cfgLowClusterE || c.energy() > cfgHighClusterE || (cfgReqRecPsRec && !inCalo(c.eta(), c.phi()))) continue; - qaState.clusterCuts[sample]->fill(7,w); + qaState.clusterCuts[sample]->fill(6, w); + if (!std::isfinite(c.eta()) || !std::isfinite(c.phi()) || !std::isfinite(c.m02()) || !std::isfinite(c.m20())) + continue; + if (c.energy() < cfgLowClusterE || c.energy() > cfgHighClusterE || (cfgReqRecPsRec && !inCalo(c.eta(), c.phi()))) + continue; + qaState.clusterCuts[sample]->fill(7, w); auto tag = tagger(c); - fillStage(sample,0,c,w,tag); - const double analysisPt=c.energy()/std::cosh(c.eta()); - fillAnalysisCommonCut(sample,0,analysisPt,c.energy(),w,tag); - if(c.m02() < cfgLowM02 || c.m02() > cfgHighM02) continue; - qaState.clusterCuts[sample]->fill(8,w); - fillStage(sample,1,c,w,tag); - fillAnalysisCommonCut(sample,1,analysisPt,c.energy(),w,tag); - const double clusterPt = c.energy()/std::cosh(c.eta()); + fillStage(sample, 0, c, w, tag); + const double analysisPt = c.energy() / std::cosh(c.eta()); + fillAnalysisCommonCut(sample, 0, analysisPt, c.energy(), w, tag); + if (c.m02() < cfgLowM02 || c.m02() > cfgHighM02) + continue; + qaState.clusterCuts[sample]->fill(8, w); + fillStage(sample, 1, c, w, tag); + fillAnalysisCommonCut(sample, 1, analysisPt, c.energy(), w, tag); + const double clusterPt = c.energy() / std::cosh(c.eta()); const auto [validMatching, matched, matchedPt, matchedP] = matchCluster(c, emctracks, sample, w, clusterPt, missingExtrapolation); - if (!validMatching) continue; - fillAnalysisCommonCut(sample,2,clusterPt,c.energy(),w,tag); - qaState.sumPt[sample]->fill(matchedPt,w); + if (!validMatching) + continue; + fillAnalysisCommonCut(sample, 2, clusterPt, c.energy(), w, tag); + qaState.sumPt[sample]->fill(matchedPt, w); bool charged = matched; - if(matched && matchedPt > 0.) { - const double ratio = c.energy()/matchedPt; - qaState.eOverSumPt[sample]->fill(ratio,w); - if(cfgRecoverHighEOverSumPt && ratio >= cfgEOverSumPtMin) charged = false; + if (matched && matchedPt > 0.) { + const double ratio = c.energy() / matchedPt; + qaState.eOverSumPt[sample]->fill(ratio, w); + if (cfgRecoverHighEOverSumPt && ratio >= cfgEOverSumPtMin) + charged = false; } // Compare strict matching and E/sum|p| recovery on the same after-M02 population. // The legacy neutral branch above uses its independently configurable E/sumPt rule. const bool recoveryPass = !matched || (matchedP > 0. && c.energy() >= cfgMatchMinEOverP * matchedP); if (!matched) { - qaState.strictPt[sample]->fill(clusterPt,w); - if (sample == 0) fillPerformance(StrictNeutral,c,w,tag); + qaState.strictPt[sample]->fill(clusterPt, w); + if (sample == 0) + fillPerformance(StrictNeutral, c, w, tag); } if (recoveryPass) { - qaState.recoveryPt[sample]->fill(clusterPt,w); - if (sample == 0) fillPerformance(MatchRecovery,c,w,tag); + qaState.recoveryPt[sample]->fill(clusterPt, w); + if (sample == 0) + fillPerformance(MatchRecovery, c, w, tag); } // Cumulative diagnostics compare the branches without changing existing selections. - const auto [analysisIsolated,analysisRatio]=measureIsolation(nullptr, - clusterPt,c.eta(),c.phi(),c.energy(),constituents,w); - const bool eligible=coneEligible(c.eta()); - fillAnalysisFinalCuts(sample,0,clusterPt,c.energy(),w,tag,!charged,eligible,analysisIsolated,analysisRatio); - fillAnalysisFinalCuts(sample,1,clusterPt,c.energy(),w,tag,!matched,eligible,analysisIsolated,analysisRatio); - fillAnalysisFinalCuts(sample,2,clusterPt,c.energy(),w,tag,recoveryPass,eligible,analysisIsolated,analysisRatio); - fillStage(sample,charged ? 3 : 2,c,w,tag); + const auto [analysisIsolated, analysisRatio] = measureIsolation(nullptr, + clusterPt, c.eta(), c.phi(), c.energy(), constituents, w); + const bool eligible = coneEligible(c.eta()); + fillAnalysisFinalCuts(sample, 0, clusterPt, c.energy(), w, tag, !charged, eligible, analysisIsolated, analysisRatio); + fillAnalysisFinalCuts(sample, 1, clusterPt, c.energy(), w, tag, !matched, eligible, analysisIsolated, analysisRatio); + fillAnalysisFinalCuts(sample, 2, clusterPt, c.energy(), w, tag, recoveryPass, eligible, analysisIsolated, analysisRatio); + fillStage(sample, charged ? 3 : 2, c, w, tag); if (!charged || recoveryPass) { const auto [isolated, passesRatio] = measureIsolation(charged ? nullptr : &qaState.recoIsolation[sample], - c.energy()/std::cosh(c.eta()), c.eta(), c.phi(), c.energy(), constituents, w); + c.energy() / std::cosh(c.eta()), c.eta(), c.phi(), c.energy(), constituents, w); if (isolated && recoveryPass) { - qaState.recoveredIsoPt[sample]->fill(clusterPt,w); - if (sample == 0) fillPerformance(MatchRecoveryIsolated,c,w,tag); + qaState.recoveredIsoPt[sample]->fill(clusterPt, w); + if (sample == 0) + fillPerformance(MatchRecoveryIsolated, c, w, tag); } if (sample == 0 && isolated && !charged) { fillPerformance(Isolated, c, w, tag); - if (passesRatio) fillPerformance(IsolatedEOverP, c, w, tag); + if (passesRatio) + fillPerformance(IsolatedEOverP, c, w, tag); } } } - if(missingExtrapolation) qaState.eventCuts[sample]->fill(6,w); + if (missingExtrapolation) + qaState.eventCuts[sample]->fill(6, w); } // MC bookkeeping helpers: keep truth classification and efficiency counting together. - template + template TruthTag tagClusterForCollision(Cluster const& c, aod::JMcParticles const& particles, int64_t mcCollisionId, double weight) { - auto tag = tagCluster(c,particles,weight); + auto tag = tagCluster(c, particles, weight); // Purity retains all contributors. Efficiency matches only photons from // this reconstructed collision's associated MC collision. for (auto idsIt = tag.efficiencyIds.begin(); idsIt != tag.efficiencyIds.end(); ++idsIt) { auto& ids = *idsIt; std::vector sameCollision; for (auto const& id : ids) { - const auto photon = particles.rawIteratorAt(id-particles.offset()); - if (photon.mcCollisionId() == mcCollisionId) sameCollision.push_back(id); + const auto photon = particles.rawIteratorAt(id - particles.offset()); + if (photon.mcCollisionId() == mcCollisionId) + sameCollision.push_back(id); } ids = std::move(sameCollision); } // efficiencyIds[1] now contains only direct contributors from this MC event. for (auto const& id : tag.efficiencyIds[1]) { - const auto direct = particles.rawIteratorAt(id-particles.offset()); - if (!std::isfinite(direct.pt()) || !inCalo(direct.eta(), direct.phi())) continue; + const auto direct = particles.rawIteratorAt(id - particles.offset()); + if (!std::isfinite(direct.pt()) || !inCalo(direct.eta(), direct.phi())) + continue; tag.fiducialDirect = true; tag.coneFiducialDirect |= coneEligible(direct.eta()); } @@ -1235,19 +1404,21 @@ struct PromptphotonQA { void fillTrackingPerformance(aod::JMcParticles const& particles, std::set const& acceptedMC, - std::array,3> const& trackingMatched) + std::array, 3> const& trackingMatched) { for (auto const& truth : particles) { - if (!truth.has_mcCollision() || !trackingTruthSelected(truth)) continue; + if (!truth.has_mcCollision() || !trackingTruthSelected(truth)) + continue; const auto mcCollision = truth.mcCollision_as(); const double w = mcCollision.weight(); - if (!std::isfinite(w)) throw std::runtime_error("Non-finite tracking truth weight"); + if (!std::isfinite(w)) + throw std::runtime_error("Non-finite tracking truth weight"); fillTrackKinematics(qaState.trackingTruth[0], truth.pt(), truth.eta(), truth.phi(), w); if (acceptedMC.count(mcCollision.globalIndex())) fillTrackKinematics(qaState.trackingTruth[1], truth.pt(), truth.eta(), truth.phi(), w); for (int i = 0; i < NTrackingMatches; ++i) if (trackingMatched[i].count(truth.globalIndex())) - fillTrackKinematics(qaState.trackingTruth[i+2], truth.pt(), truth.eta(), truth.phi(), w); + fillTrackKinematics(qaState.trackingTruth[i + 2], truth.pt(), truth.eta(), truth.phi(), w); } } @@ -1258,35 +1429,46 @@ struct PromptphotonQA { // Truth fiducial region is explicit and shared by numerator/denominator. // No reconstructed energy or shower-shape requirement is applied to truth. for (auto const& photon : particles) { - if (photon.pdgCode()!=kGamma || !std::isfinite(photon.pt()) || !std::isfinite(photon.eta()) || !std::isfinite(photon.phi())) continue; - if (!photon.has_mcCollision()) continue; + if (photon.pdgCode() != kGamma || !std::isfinite(photon.pt()) || !std::isfinite(photon.eta()) || !std::isfinite(photon.phi())) + continue; + if (!photon.has_mcCollision()) + continue; const auto mcCollision = photon.mcCollision_as(); const double weight = mcCollision.weight(); - if (!std::isfinite(weight)) throw std::runtime_error("Non-finite efficiency MC weight"); - const auto result = classifyOriginal(photon,particles); + if (!std::isfinite(weight)) + throw std::runtime_error("Non-finite efficiency MC weight"); + const auto result = classifyOriginal(photon, particles); const bool selectedEvent = !cfgGenReqRec || reconstructedMC.count(mcCollision.globalIndex()); const bool direct = !result.directHardIndices.empty(); - if (direct && selectedEvent) qaState.directBeforeCalo->fill(photon.pt(),weight); - if (!inCalo(photon.eta(), photon.phi())) continue; - if (direct && selectedEvent) qaState.directInCalo->fill(photon.pt(),weight); - const std::array belongs{isPrompt(result.origin),!result.directHardIndices.empty()}; - for (int g=0;g<2;++g) { - if (!belongs[g]) continue; + if (direct && selectedEvent) + qaState.directBeforeCalo->fill(photon.pt(), weight); + if (!inCalo(photon.eta(), photon.phi())) + continue; + if (direct && selectedEvent) + qaState.directInCalo->fill(photon.pt(), weight); + const std::array belongs{isPrompt(result.origin), !result.directHardIndices.empty()}; + for (int g = 0; g < 2; ++g) { + if (!belongs[g]) + continue; const bool eligible = coneEligible(photon.eta()); if (selectedEvent) { - qaState.selectedEventTruth[g]->fill(photon.pt(),weight); - if (eligible) qaState.selectedConeTruth[g]->fill(photon.pt(),weight); + qaState.selectedEventTruth[g]->fill(photon.pt(), weight); + if (eligible) + qaState.selectedConeTruth[g]->fill(photon.pt(), weight); } if (eligible) { - qaState.coneAll[g]->fill(photon.pt(),weight); - if (acceptedMC.count(mcCollision.globalIndex())) qaState.coneAccepted[g]->fill(photon.pt(),weight); + qaState.coneAll[g]->fill(photon.pt(), weight); + if (acceptedMC.count(mcCollision.globalIndex())) + qaState.coneAccepted[g]->fill(photon.pt(), weight); } - qaState.efficiencyAll[g]->fill(photon.pt(),weight); - if (acceptedMC.count(mcCollision.globalIndex())) qaState.efficiencyAccepted[g]->fill(photon.pt(),weight); - for (int st=0;stfill(photon.pt(), weight); + if (acceptedMC.count(mcCollision.globalIndex())) + qaState.efficiencyAccepted[g]->fill(photon.pt(), weight); + for (int st = 0; st < nPerformanceStages; ++st) // Pair isolated numerators with the truth cone-eligible denominators. if ((eligible || !isolatedStage(st)) && - qaState.efficiencyMatched[g][st].count(photon.globalIndex())) qaState.efficiencyNumerator[g][st]->fill(photon.pt(),weight); + qaState.efficiencyMatched[g][st].count(photon.globalIndex())) + qaState.efficiencyNumerator[g][st]->fill(photon.pt(), weight); } } } @@ -1298,22 +1480,24 @@ struct PromptphotonQA { { if (trackLabels.size() != tracks.size() || trackLabels.offset() != tracks.offset()) throw std::runtime_error("JMcTrackLbs must be row-aligned with JTracks; enable the JE MC track-label producer"); - std::array,3> trackingMatched; + std::array, 3> trackingMatched; for (auto group = qaState.efficiencyMatched.begin(); group != qaState.efficiencyMatched.end(); ++group) - for (auto stage = group->begin(); stage != group->end(); ++stage) stage->clear(); + for (auto stage = group->begin(); stage != group->end(); ++stage) + stage->clear(); std::set reconstructedMC, acceptedMC; for (auto const& collision : collisions) { - if (!collision.has_mcCollision()) continue; + if (!collision.has_mcCollision()) + continue; const auto mcCollision = collision.mcCollision_as(); reconstructedMC.insert(mcCollision.globalIndex()); const double weight = mcCollision.weight(); - auto eventClusters = clusters.sliceBy(clustersPerCollision,collision.globalIndex()); + auto eventClusters = clusters.sliceBy(clustersPerCollision, collision.globalIndex()); auto rawEventClusters = rawClusters.sliceBy(clustersPerCollision, collision.globalIndex()); - auto eventTracks = tracks.sliceBy(tracksPerCollision,collision.globalIndex()); + auto eventTracks = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); processReco(collision, eventClusters, rawEventClusters, eventTracks, emctracks, 0, weight, - [&](auto const& c) { - return tagClusterForCollision(c, particles, mcCollision.globalIndex(), weight); - }); + [&](auto const& c) { + return tagClusterForCollision(c, particles, mcCollision.globalIndex(), weight); + }); if (qaState.recoEventAccepted) { acceptedMC.insert(mcCollision.globalIndex()); collectTrackingMatches(eventTracks, trackLabels, particles, mcCollision.globalIndex(), emctracks, weight, trackingMatched); @@ -1323,21 +1507,26 @@ struct PromptphotonQA { fillPhotonPerformance(particles, reconstructedMC, acceptedMC); } - template + template void fillGenCollision(Collision const& collision, ParticleTable const& particles) { const double weight = collision.weight(); - if(!std::isfinite(weight)) throw std::runtime_error("Non-finite generator weight"); - qaState.genEvents->fill(0.5,weight); + if (!std::isfinite(weight)) + throw std::runtime_error("Non-finite generator weight"); + qaState.genEvents->fill(0.5, weight); std::vector constituents; for (auto const& particle : particles) { - if (!particle.isPhysicalPrimary() || !std::isfinite(particle.eta()) || !std::isfinite(particle.phi())) continue; + if (!particle.isPhysicalPrimary() || !std::isfinite(particle.eta()) || !std::isfinite(particle.phi())) + continue; const auto phi = phi02pi(particle.phi()); qaState.allPrimaryMap->fill(particle.eta(), phi, weight); - if (particle.pdgCode() == kGamma) qaState.allPhotonMap->fill(particle.eta(), phi, weight); - if (!std::isfinite(particle.pt()) || particle.pt() < cfgtrkMinPt || std::abs(particle.eta()) > cfgtrkMaxEta) continue; + if (particle.pdgCode() == kGamma) + qaState.allPhotonMap->fill(particle.eta(), phi, weight); + if (!std::isfinite(particle.pt()) || particle.pt() < cfgtrkMinPt || std::abs(particle.eta()) > cfgtrkMaxEta) + continue; auto* species = pdg->GetParticle(particle.pdgCode()); - if (!species || species->Charge() == 0.) continue; + if (!species || species->Charge() == 0.) + continue; constituents.push_back({particle.pt(), particle.eta(), particle.phi()}); qaState.genParticleMap->fill(particle.eta(), phi02pi(particle.phi()), weight); } @@ -1345,9 +1534,10 @@ struct PromptphotonQA { for (auto const& photon : particles) { if (photon.pdgCode() != kGamma || !photon.isPhysicalPrimary() || photon.getGenStatusCode() < 0 || !std::isfinite(photon.pt()) || !std::isfinite(photon.e()) || !std::isfinite(photon.phi()) || - !inCalo(photon.eta(), photon.phi())) continue; + !inCalo(photon.eta(), photon.phi())) + continue; const auto phi = phi02pi(photon.phi()); - const auto result = classifyOriginal(photon,particles); + const auto result = classifyOriginal(photon, particles); if (isPrompt(result.origin)) { qaState.genPromptEta->fill(photon.eta(), weight); qaState.genPromptPhi->fill(phi, weight); @@ -1355,35 +1545,40 @@ struct PromptphotonQA { } const bool direct = !result.directHardIndices.empty(); const bool frag = !result.fragmentationHardIndices.empty(); - qaState.ancestryStatus->fill(direct ? (frag ? 3 : 1) : (frag ? 2 : (result.leptonBranch ? 4 : 0)),weight); + qaState.ancestryStatus->fill(direct ? (frag ? 3 : 1) : (frag ? 2 : (result.leptonBranch ? 4 : 0)), weight); // Groups overlap only for mixed ancestry. Each group counts the photon once. - for (int g=0;g<2;++g) { - if (!(g == 0 ? direct : frag)) continue; + for (int g = 0; g < 2; ++g) { + if (!(g == 0 ? direct : frag)) + continue; qaState.genPhotonMap[g]->fill(photon.eta(), phi, weight); - qaState.genPhi[g]->fill(phi,weight); qaState.genEta[g]->fill(photon.eta(),weight); qaState.genPt[g]->fill(photon.pt(),weight); + qaState.genPhi[g]->fill(phi, weight); + qaState.genEta[g]->fill(photon.eta(), weight); + qaState.genPt[g]->fill(photon.pt(), weight); } } } //****************************************** - //ANALYSIS HELPERS COMPLETE!!! - //NOW FOR THE MAIN PROCESS FUNCTIONS + // ANALYSIS HELPERS COMPLETE!!! + // NOW FOR THE MAIN PROCESS FUNCTIONS //****************************************** void processMCRec(MCCollisions const& collisions, FilteredMCClusters const& clusters, aod::JClusters const& rawClusters, RecoTracks const& tracks, TrackCandidates const&, aod::JMcParticles const& particles, BcCandidates const&, aod::JEMCTracks const& emctracks, aod::JMcCollisions const&, - aod::JMcTrackLbs const& trackLabels) { + aod::JMcTrackLbs const& trackLabels) + { fillMCReco(collisions, clusters, rawClusters, tracks, particles, emctracks, trackLabels); } - + PROCESS_SWITCH(PromptphotonQA, processMCRec, "Weighted reconstructed MC QA", false); void processData(DataCollisions::iterator const& collision, FilteredDataClusters const& clusters, aod::JClusters const& rawClusters, RecoTracks const& tracks, TrackCandidates const&, BcCandidates const&, - aod::JEMCTracks const& emctracks) { + aod::JEMCTracks const& emctracks) + { processReco(collision, clusters, rawClusters, tracks, emctracks, 1, 1., [](auto const&) { return TruthTag{}; }); } PROCESS_SWITCH(PromptphotonQA, processData, "Real-data reconstructed QA", false); @@ -1401,12 +1596,14 @@ struct PromptphotonQA { soa::SmallGroups const& recoCollisions, aod::JMcParticles const& particles) { - if (cfgGenReqRec && recoCollisions.size() == 0) return; + if (cfgGenReqRec && recoCollisions.size() == 0) + return; fillGenCollision(collision, particles); } PROCESS_SWITCH(PromptphotonQA, processMCGenJE, "Generator QA on JE truth with optional reconstructed counterpart", false); }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ return WorkflowSpec{adaptAnalysisTask(cfgc)}; }