Skip to content

Commit 0d0736a

Browse files
SteffimroSteffimro
committed
add adjustable QA settings + qinv gate for mixed pairs + fix headers
1 parent f8565ec commit 0d0736a

1 file changed

Lines changed: 58 additions & 14 deletions

File tree

‎PWGEM/PhotonMeson/Tasks/photonhbt.cxx‎

Lines changed: 58 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@
2929
#include <CCDB/BasicCCDBManager.h>
3030
#include <CommonConstants/MathConstants.h>
3131
#include <DataFormatsParameters/GRPMagField.h>
32+
#include <Framework/ASoA.h>
3233
#include <Framework/ASoAHelpers.h>
3334
#include <Framework/AnalysisDataModel.h>
3435
#include <Framework/AnalysisHelpers.h>
@@ -47,6 +48,7 @@
4748
#include <Math/Vector3Dfwd.h>
4849
#include <Math/Vector4D.h> // IWYU pragma: keep
4950
#include <Math/Vector4Dfwd.h>
51+
#include <TH1.h>
5052
#include <TPDGCode.h>
5153
#include <TString.h>
5254

@@ -222,6 +224,9 @@ struct Photonhbt {
222224
Configurable<float> cfgMaxQinvForQA{"cfgMaxQinvForQA", 0.1f, "fill per-step pair QA histograms only when q_inv < this value"};
223225
Configurable<float> cfgMaxQinvForFullRange{"cfgMaxQinvForFullRange", 0.3f, "fill full-range histograms only when q_inv < this value"};
224226
Configurable<float> cfgMaxQinvForMCQA{"cfgMaxQinvForMCQA", 0.3f, "fill MC truth 1D histograms only when q_inv < this value"};
227+
Configurable<int> cfgQaLevel{"cfgQaLevel", 2, "QA: 0 no QA, 1 standard, 2 full diagnostics with leg information etc."};
228+
Configurable<bool> cfgFillDRDZSparse{"cfgFillDRDZSparse", true, "book/fill the FullRange |R1-R2|-Deltaz-qinv sparse (large; needs cfgQaLevel>=2)"};
229+
Configurable<float> cfgMaxQinvForProcessing{"cfgMaxQinvForProcessing", 0.5, "skip mixed pairs with q_inv above this before building observables"};
225230
} qaflags;
226231

227232
// ─── HBT analysis mode ───────────────────────────────────────────────────────────
@@ -507,9 +512,21 @@ struct Photonhbt {
507512
// INITS
508513
/*************************************************/
509514

515+
516+
bool mDoPairQa{true}, mDoSinglePhotonQa{true}, mDoLegPairQA{true};
517+
bool mDoPairSepQA{true}, mFillDRDZSparse{true};
518+
510519
void init(InitContext& context)
511520
{
512521
isMC = context.mOptions.get<bool>("processMC");
522+
const int qaLevel = qaflags.cfgQaLevel.value;
523+
mDoPairQa = qaflags.doPairQa.value && qaLevel >= 1;
524+
mDoSinglePhotonQa = qaflags.doSinglePhotonQa.value && qaLevel >= 1;
525+
mDoLegPairQA = qaflags.doLegPairQA.value && qaLevel >= 2;
526+
mDoPairSepQA = pairsep.cfgDoPairSepQA.value && qaLevel >= 2;
527+
mFillDRDZSparse = qaflags.cfgFillDRDZSparse.value && qaLevel >= 2;
528+
LOGF(info, "photonhbt QA level %d -> pairQA %d, singlePhotonQA %d, legPairQA %d, pairSep %d, dRdZ sparse %d",
529+
qaLevel, mDoPairQa, mDoSinglePhotonQa, mDoLegPairQA, mDoPairSepQA, mFillDRDZSparse);
513530
mRunNumber = 0;
514531
parseBins(mixing.confVtxBins, ztxBinEdges);
515532
parseBins(mixing.confCentBins, centBinEdges);
@@ -676,7 +693,7 @@ struct Photonhbt {
676693

677694
if (isMC) {
678695
addPairMCHistograms();
679-
if (qaflags.doLegPairQA) {
696+
if (mDoLegPairQA) {
680697
addLegPairMCHistograms();
681698
}
682699
addTruthMCHistograms();
@@ -763,8 +780,10 @@ struct Photonhbt {
763780
}
764781

765782
fRegistryPairQA.addClone("Pair/same/QA/", "Pair/mix/QA/");
766-
addLegPairQAForStep("Pair/same/QA/Before/");
767-
addLegPairQAForStep("Pair/same/QA/AfterPairCuts/");
783+
if (mDoLegPairQA) {
784+
addLegPairQAForStep("Pair/same/QA/Before/");
785+
addLegPairQAForStep("Pair/same/QA/AfterPairCuts/");
786+
}
768787
}
769788

770789
void addPairMCHistograms()
@@ -1005,7 +1024,9 @@ struct Photonhbt {
10051024
fRegistryCF.add((path + "hDeltaR3DVsQinv").c_str(), "#Delta r_{3D} vs q_{inv};q_{inv} (GeV/c);#Delta r_{3D} (cm)", kTH2D, {axisQinv, axisDeltaR3D}, true);
10061025
fRegistryCF.add((path + "hQinvVsCent").c_str(), "q_{inv} vs centrality;centrality (%);q_{inv} (GeV/c)", kTH2D, {axisCentQA, axisQinv}, true);
10071026
fRegistryCF.add((path + "hQinvVsOccupancy").c_str(), "q_{inv} vs occupancy;occupancy;q_{inv} (GeV/c)", kTH2D, {axisOccupancy, axisQinv}, true);
1008-
fRegistryCF.add((path + "hSparseDeltaRDeltaZQinv").c_str(), "|R_{1}-R_{2}|,#Delta z,q_{inv}", kTHnSparseD, {axisDeltaR, axisDeltaZ, axisQinv}, true);
1027+
if (mFillDRDZSparse) {
1028+
fRegistryCF.add((path + "hSparseDeltaRDeltaZQinv").c_str(), "|R_{1}-R_{2}|,#Delta z,q_{inv}", kTHnSparseD, {axisDeltaR, axisDeltaZ, axisQinv}, true);
1029+
}
10091030
fRegistryCF.add((path + "hDeltaRCosOAVsQinv").c_str(), "#Delta r/cos(#theta_{op}/2) vs q_{inv};q_{inv} (GeV/c);#Delta r/cos(#theta_{op}/2) (cm)", kTH2D, {axisQinv, {100, 0, 100}}, true);
10101031
}
10111032

@@ -1266,6 +1287,19 @@ struct Photonhbt {
12661287
return s;
12671288
}
12681289

1290+
template <typename TG1, typename TG2>
1291+
[[nodiscard]] inline bool passFastQinvGate(TG1 const& g1, TG2 const& g2) const
1292+
{
1293+
const float qMax = qaflags.cfgMaxQinvForProcessing.value;
1294+
if (qMax > 1e9f) { // o2-linter: disable=magic-number (skip if non-sensical value is chosen)
1295+
return true;
1296+
}
1297+
const float dEta = g1.eta() - g2.eta();
1298+
const float dPhi = RecoDecay::constrainAngle(g1.phi() - g2.phi(), -o2::constants::math::PI);
1299+
const float q2 = 2.f * g1.pt() * g2.pt() * (std::cosh(dEta) - std::cos(dPhi));
1300+
return q2 <= qMax * qMax;
1301+
}
1302+
12691303
[[nodiscard]] inline bool passLegSepCut(PairSep const& s) const
12701304
{
12711305
if (!pairsep.cfgDoLegSepCut.value) {
@@ -1386,13 +1420,15 @@ struct Photonhbt {
13861420
fRegistryCF.fill(HIST(base) + HIST("hDeltaR3DVsQinv"), obs.qinv, obs.deltaR3D);
13871421
fRegistryCF.fill(HIST(base) + HIST("hQinvVsCent"), cent, obs.qinv);
13881422
fRegistryCF.fill(HIST(base) + HIST("hQinvVsOccupancy"), occupancy, obs.qinv);
1389-
fRegistryCF.fill(HIST(base) + HIST("hSparseDeltaRDeltaZQinv"), obs.deltaR, obs.deltaZ, obs.qinv);
1423+
if (mFillDRDZSparse) {
1424+
fRegistryCF.fill(HIST(base) + HIST("hSparseDeltaRDeltaZQinv"), obs.deltaR, obs.deltaZ, obs.qinv);
1425+
}
13901426
}
13911427

13921428
template <int ev_id, bool after_cut = false>
13931429
inline void fillPairSep(PairSep const& s, PairQAObservables const& obs)
13941430
{
1395-
if (!pairsep.cfgDoPairSepQA.value) {
1431+
if (!mDoPairSepQA) {
13961432
return;
13971433
}
13981434
const float limit = pairsep.cfgPairSepMaxQinv.value;
@@ -1422,7 +1458,7 @@ struct Photonhbt {
14221458
template <int step_id, typename TPhoton>
14231459
inline void fillSinglePhotonQAStep(TPhoton const& g)
14241460
{
1425-
if (!qaflags.doSinglePhotonQa) {
1461+
if (!mDoSinglePhotonQa) {
14261462
return;
14271463
}
14281464
constexpr auto base = singlePhotonQAPrefix<step_id>();
@@ -1519,7 +1555,7 @@ struct Photonhbt {
15191555
template <int ev_id, int step_id>
15201556
inline void fillPairQAStep(PairQAObservables const& o, float /*cent*/, float /*occupancy*/)
15211557
{
1522-
if (!qaflags.doPairQa) {
1558+
if (!mDoPairQa) {
15231559
return;
15241560
}
15251561
constexpr auto base = qaPrefix<ev_id, step_id>();
@@ -1559,7 +1595,7 @@ struct Photonhbt {
15591595
template <int step_id>
15601596
inline void fillLegPairQAStep(LegPairObservables const& lo, float kt)
15611597
{
1562-
if (!qaflags.doPairQa) {
1598+
if (!mDoLegPairQA) {
15631599
return;
15641600
}
15651601
constexpr auto base = qaPrefix<0, step_id>();
@@ -1948,7 +1984,7 @@ struct Photonhbt {
19481984
auto keyDFCollision = std::make_pair(ndf, collision.globalIndex());
19491985
auto photons1Coll = photons1.sliceBy(perCollision1, collision.globalIndex());
19501986
auto photons2Coll = photons2.sliceBy(perCollision2, collision.globalIndex());
1951-
if (qaflags.doSinglePhotonQa) {
1987+
if (mDoSinglePhotonQa) {
19521988
for (const auto& g : photons1Coll) {
19531989
if (cut1.template IsSelected<decltype(g), TSubInfos1>(g)) {
19541990
fillSinglePhotonQAStep<0>(g);
@@ -2039,7 +2075,7 @@ struct Photonhbt {
20392075
addToPool(g1, pwl1);
20402076
addToPool(g2, pwl2);
20412077
}
2042-
if (qaflags.doSinglePhotonQa) {
2078+
if (mDoSinglePhotonQa) {
20432079
for (const auto& g : photons1Coll) {
20442080
if (cut1.template IsSelected<decltype(g), TSubInfos1>(g)) {
20452081
if (idsAfterPairCuts.contains(g.globalIndex())) {
@@ -2070,6 +2106,9 @@ struct Photonhbt {
20702106
if (!passAsymmetryCut(g1.pt(), g2.pt())) {
20712107
continue;
20722108
}
2109+
if (!passFastQinvGate(g1, g2)) {
2110+
continue;
2111+
}
20732112
auto obs = buildPairQAObservables(g1, g2);
20742113
if (!obs.valid) {
20752114
continue;
@@ -2170,7 +2209,7 @@ struct Photonhbt {
21702209
auto keyBin = std::make_tuple(zbin, centbin, epbin, occbin);
21712210
auto keyDFCollision = std::make_pair(ndf, collision.globalIndex());
21722211
auto photonsColl = photons.sliceBy(perCollision, collision.globalIndex());
2173-
if (qaflags.doSinglePhotonQa) {
2212+
if (mDoSinglePhotonQa) {
21742213
for (const auto& g : photonsColl) {
21752214
if (cut.template IsSelected<decltype(g), TLegs>(g)) {
21762215
fillSinglePhotonQAStep<0>(g);
@@ -2265,7 +2304,9 @@ struct Photonhbt {
22652304
} else {
22662305
const bool doMCQA = passQinvMCQAGate(obs.qinv);
22672306
fillMCPairQA<false>(truthType, obs, doQA, doMCQA);
2268-
fillLegPairMC(truthType, legObs, obs.kt);
2307+
if (mDoLegPairQA) {
2308+
fillLegPairMC(truthType, legObs, obs.kt);
2309+
}
22692310
if (doFR) {
22702311
fillMCPairQAFullRange<false>(truthType, obs);
22712312
}
@@ -2318,7 +2359,7 @@ struct Photonhbt {
23182359
addToPool(g1, pwl1);
23192360
addToPool(g2, pwl2);
23202361
}
2321-
if (qaflags.doSinglePhotonQa) {
2362+
if (mDoSinglePhotonQa) {
23222363
for (const auto& g : photonsColl) {
23232364
if (cut.template IsSelected<decltype(g), TLegs>(g)) {
23242365
if (idsAfterPairCuts.contains(g.globalIndex())) {
@@ -2349,6 +2390,9 @@ struct Photonhbt {
23492390
if (!passAsymmetryCut(g1.pt(), g2.pt())) {
23502391
continue;
23512392
}
2393+
if (!passFastQinvGate(g1, g2)) {
2394+
continue;
2395+
}
23522396
auto obs = buildPairQAObservables(g1, g2);
23532397
if (!obs.valid) {
23542398
continue;

0 commit comments

Comments
 (0)