Skip to content

Commit 48a02b0

Browse files
mhemmer-cernmhemmer-cern
authored
[PWGEM,PWGEM-36] Pi0 Flow Task: Enable running Mixed and rotation tog… (#10050)
1 parent 91a41c9 commit 48a02b0

File tree

1 file changed

+133
-27
lines changed

1 file changed

+133
-27
lines changed

‎PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx‎

Lines changed: 133 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -95,14 +95,16 @@ struct TaskPi0FlowEMC {
9595
Configurable<int> cfgEMCalMapLevelSameEvent{"cfgEMCalMapLevelSameEvent", 1, "Different levels of correction for the same event, the smaller number includes the level of the higher number (4: none, 3: only inside EMCal, 2: exclude bad channels, 1: remove edges)"};
9696
Configurable<float> cfgRotAngle{"cfgRotAngle", std::move(const_cast<float&>(o2::constants::math::PIHalf)), "Angle used for the rotation method"};
9797
Configurable<int> cfgDistanceToEdge{"cfgDistanceToEdge", 1, "Distance to edge in cells required for rotated cluster to be accepted"};
98+
Configurable<bool> cfgDoM02{"cfgDoM02", false, "Flag to enable flow vs M02 for single photons"};
9899

99100
// configurable axis
100-
ConfigurableAxis thnConfigAxisInvMass{"thnConfigAxisInvMass", {200, 0.0, 0.4}, ""};
101+
ConfigurableAxis thnConfigAxisInvMass{"thnConfigAxisInvMass", {400, 0.0, 0.8}, ""};
101102
ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {100, 0., 20.}, ""};
102103
ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {20, 0., 100.}, ""};
103104
ConfigurableAxis thnConfigAxisCosNPhi{"thnConfigAxisCosNPhi", {100, -1., 1.}, ""};
104105
ConfigurableAxis thnConfigAxisCosDeltaPhi{"thnConfigAxisCosDeltaPhi", {100, -1., 1.}, ""};
105106
ConfigurableAxis thnConfigAxisScalarProd{"thnConfigAxisScalarProd", {100, -5., 5.}, ""};
107+
ConfigurableAxis thnConfigAxisM02{"thnConfigAxisM02", {200, 0., 5.}, ""};
106108

107109
EMPhotonEventCut fEMEventCut;
108110
struct : ConfigurableGroup {
@@ -276,6 +278,7 @@ struct TaskPi0FlowEMC {
276278
const AxisSpec thnAxisCosNPhi{thnConfigAxisCosNPhi, Form("cos(%d#varphi)", harmonic.value)};
277279
const AxisSpec thnAxisCosDeltaPhi{thnConfigAxisCosDeltaPhi, Form("cos(%d(#varphi - #Psi_{sub}))", harmonic.value)};
278280
const AxisSpec thnAxisScalarProd{thnConfigAxisScalarProd, "SP"};
281+
const AxisSpec thnAxisM02{thnConfigAxisM02, "M_{02}"};
279282
const AxisSpec thAxisTanThetaPhi{mesonConfig.thConfigAxisTanThetaPhi, "atan(#Delta#theta/#Delta#varphi)"};
280283
const AxisSpec thAxisClusterEnergy{thnConfigAxisPt, "#it{E} (GeV)"};
281284
const AxisSpec thAxisAlpha{100, -1., +1, "#alpha"};
@@ -291,8 +294,14 @@ struct TaskPi0FlowEMC {
291294
const AxisSpec thAxisSN{8, 0.5, 8.5, "#it{s}_{n}"};
292295
const AxisSpec thAxisCPUTime{1000, 0, 10000, "#it{t} (#mus)"};
293296

297+
const AxisSpec thnAxisMixingVtx{mixingConfig.cfgVtxBins, "#it{z} (cm)"};
298+
const AxisSpec thnAxisMixingCent{mixingConfig.cfgCentBins, "Centrality (%)"};
299+
const AxisSpec thnAxisMixingEP{mixingConfig.cfgEPBins, Form("cos(%d#varphi)", harmonic.value)};
300+
294301
registry.add("hSparsePi0Flow", "THn for SP", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisCent, thnAxisScalarProd});
295-
registry.add("hSparseBkgFlow", "THn for SP", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisCent, thnAxisScalarProd});
302+
registry.add("hSparseBkgRotFlow", "THn for SP", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisCent, thnAxisScalarProd});
303+
registry.add("hSparseBkgMixFlow", "THn for SP", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisCent, thnAxisScalarProd});
304+
registry.add("h3DMixingCount", "THn Event Mixing QA", HistType::kTH3D, {thnAxisMixingVtx, thnAxisMixingCent, thnAxisMixingEP});
296305
auto hClusterCuts = registry.add<TH1>("hClusterCuts", "hClusterCuts;;Counts", kTH1D, {{6, 0.5, 6.5}}, false);
297306
hClusterCuts->GetXaxis()->SetBinLabel(1, "in");
298307
hClusterCuts->GetXaxis()->SetBinLabel(2, "opening angle");
@@ -301,6 +310,14 @@ struct TaskPi0FlowEMC {
301310
hClusterCuts->GetXaxis()->SetBinLabel(5, "conversion cut");
302311
hClusterCuts->GetXaxis()->SetBinLabel(6, "out");
303312

313+
auto hClusterCutsMixed = registry.add<TH1>("hClusterCutsMixed", "hClusterCutsMixed;;Counts", kTH1D, {{6, 0.5, 6.5}}, false);
314+
hClusterCutsMixed->GetXaxis()->SetBinLabel(1, "in");
315+
hClusterCutsMixed->GetXaxis()->SetBinLabel(2, "opening angle");
316+
hClusterCutsMixed->GetXaxis()->SetBinLabel(3, "#it{M}_{#gamma#gamma}");
317+
hClusterCutsMixed->GetXaxis()->SetBinLabel(4, "#it{p}_{T}");
318+
hClusterCutsMixed->GetXaxis()->SetBinLabel(5, "conversion cut");
319+
hClusterCutsMixed->GetXaxis()->SetBinLabel(6, "out");
320+
304321
if (saveSPResoHist) {
305322
registry.add("spReso/hSpResoFT0cFT0a", "hSpResoFT0cFT0a; centrality; Q_{FT0c} #bullet Q_{FT0a}", HistType::kTH2D, {thnAxisCent, thnConfigAxisScalarProd});
306323
registry.add("spReso/hSpResoFT0cTPCpos", "hSpResoFT0cTPCpos; centrality; Q_{FT0c} #bullet Q_{TPCpos}", HistType::kTH2D, {thnAxisCent, thnConfigAxisScalarProd});
@@ -352,10 +369,6 @@ struct TaskPi0FlowEMC {
352369
hCollisionEMCCheck->GetXaxis()->SetBinLabel(5, "EMC MB Readout but no clusters");
353370
hCollisionEMCCheck->GetXaxis()->SetBinLabel(6, "No EMC MB Readout but has clusters");
354371
hCollisionEMCCheck->GetXaxis()->SetBinLabel(7, "No EMC MB Readout and no clusters");
355-
registry.add("LED/hMult", "multiplicity in LED events", HistType::kTH1D, {thAxisMult});
356-
registry.add("LED/hClusterEtaPhi", "hClusterEtaPhi", HistType::kTH2D, {thAxisPhi, thAxisEta});
357-
registry.add("LED/clusterTimeVsE", "Cluster time vs energy", HistType::kTH2D, {thAxisTime, thAxisEnergy});
358-
registry.add("LED/hNCell", "hNCell", HistType::kTH1D, {thAxisNCell});
359372
}
360373

361374
if (emccuts.cfgEnableQA) {
@@ -369,6 +382,11 @@ struct TaskPi0FlowEMC {
369382
registry.add("hAlphaPt", "Histo of meson asymmetry vs pT", HistType::kTH2D, {thAxisAlpha, thnAxisPt});
370383
registry.add("mesonQA/hClusterEtaPhiBefore", "hClusterEtaPhiBefore", HistType::kTH2D, {thAxisPhi, thAxisEta});
371384
registry.add("mesonQA/hClusterEtaPhiAfter", "hClusterEtaPhiAfter", HistType::kTH2D, {thAxisPhi, thAxisEta});
385+
registry.add("hInvMassPtMixed", "Histo for inv pair mass vs pt for mixed event", HistType::kTH2D, {thnAxisInvMass, thnAxisPt});
386+
registry.add("hTanThetaPhiMixed", "Histo for identification of conversion cluster for mixed event", HistType::kTH2D, {thnAxisInvMass, thAxisTanThetaPhi});
387+
registry.add("hAlphaPtMixed", "Histo of meson asymmetry vs pT for mixed event", HistType::kTH2D, {thAxisAlpha, thnAxisPt});
388+
registry.add("mesonQA/hClusterEtaPhiBeforeMixed", "hClusterEtaPhiBefore for mixed event", HistType::kTH2D, {thAxisPhi, thAxisEta});
389+
registry.add("mesonQA/hClusterEtaPhiAfterMixed", "hClusterEtaPhiAfter for mixed event", HistType::kTH2D, {thAxisPhi, thAxisEta});
372390
if (cfgDoRotation) {
373391
registry.add("mesonQA/hClusterBackEtaPhiBefore", "hClusterBackEtaPhiBefore", HistType::kTH2D, {thAxisPhi, thAxisEta});
374392
registry.add("mesonQA/hClusterBackEtaPhiAfter", "hClusterBackEtaPhiAfter", HistType::kTH2D, {thAxisPhi, thAxisEta});
@@ -380,6 +398,10 @@ struct TaskPi0FlowEMC {
380398
registry.add("hSparseCalibBack", "THn for Calib background", HistType::kTHnSparseF, {thnAxisInvMass, thAxisEnergyCalib, thnAxisCent});
381399
}
382400

401+
if (cfgDoM02.value) {
402+
registry.add("hSparseFlow", "THn for SP", HistType::kTHnSparseF, {thnAxisM02, thnAxisPt, thnAxisCent, thnAxisScalarProd});
403+
}
404+
383405
ccdb->setURL(ccdbUrl);
384406
ccdb->setCaching(true);
385407
ccdb->setLocalObjectValidityChecking();
@@ -417,7 +439,7 @@ struct TaskPi0FlowEMC {
417439
float& cent,
418440
float& sp)
419441
{
420-
static constexpr std::string_view HistTypes[2] = {"hSparsePi0Flow", "hSparseBkgFlow"};
442+
static constexpr std::string_view HistTypes[3] = {"hSparsePi0Flow", "hSparseBkgRotFlow", "hSparseBkgMixFlow"};
421443
registry.fill(HIST(HistTypes[histType]), mass, pt, cent, sp);
422444
}
423445

@@ -693,10 +715,10 @@ struct TaskPi0FlowEMC {
693715
float dTheta = photon1.Theta() - photon3.Theta();
694716
float dPhi = photon1.Phi() - photon3.Phi();
695717
if (mesonConfig.enableTanThetadPhi && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) {
696-
registry.fill(HIST("hSparseBkgFlow"), mother1.M(), mother1.Pt(), cent, scalprodCand1);
718+
registry.fill(HIST("hSparseBkgRotFlow"), mother1.M(), mother1.Pt(), cent, scalprodCand1);
697719
}
698720
} else {
699-
registry.fill(HIST("hSparseBkgFlow"), mother1.M(), mother1.Pt(), cent, scalprodCand1);
721+
registry.fill(HIST("hSparseBkgRotFlow"), mother1.M(), mother1.Pt(), cent, scalprodCand1);
700722
}
701723
}
702724
}
@@ -716,10 +738,10 @@ struct TaskPi0FlowEMC {
716738
float dTheta = photon2.Theta() - photon3.Theta();
717739
float dPhi = photon2.Phi() - photon3.Phi();
718740
if (mesonConfig.enableTanThetadPhi && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) {
719-
registry.fill(HIST("hSparseBkgFlow"), mother2.M(), mother2.Pt(), cent, scalprodCand2);
741+
registry.fill(HIST("hSparseBkgRotFlow"), mother2.M(), mother2.Pt(), cent, scalprodCand2);
720742
}
721743
} else {
722-
registry.fill(HIST("hSparseBkgFlow"), mother2.M(), mother2.Pt(), cent, scalprodCand2);
744+
registry.fill(HIST("hSparseBkgRotFlow"), mother2.M(), mother2.Pt(), cent, scalprodCand2);
723745
}
724746
}
725747
}
@@ -853,12 +875,6 @@ struct TaskPi0FlowEMC {
853875
if (photonsPerCollision.size() > 0) {
854876
registry.fill(HIST("hCollisionEMCCheck"), 3.); // has EMC cluster
855877
registry.fill(HIST("hCollisionEMCCheck"), 6.); // has no EMC read out and clusters
856-
registry.fill(HIST("LED/hMult"), collision.multFT0C());
857-
for (const auto& photon : photonsPerCollision) {
858-
registry.fill(HIST("LED/hClusterEtaPhi"), photon.phi(), photon.eta());
859-
registry.fill(HIST("LED/clusterTimeVsE"), photon.time(), photon.e());
860-
registry.fill(HIST("LED/hNCell"), photon.nCells());
861-
}
862878
} else {
863879
registry.fill(HIST("hCollisionEMCCheck"), 7.); // has no EMC read out and no clusters
864880
}
@@ -1004,6 +1020,7 @@ struct TaskPi0FlowEMC {
10041020
initCCDB(c1);
10051021
runBefore = runNow;
10061022
}
1023+
registry.fill(HIST("h3DMixingCount"), c1.posZ(), getCentrality(c1), c1.ep2ft0m());
10071024
for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(clusters1, clusters2))) {
10081025
if (!(fEMCCut.IsSelected<EMCalPhotons::iterator>(g1)) || !(fEMCCut.IsSelected<EMCalPhotons::iterator>(g2))) {
10091026
continue;
@@ -1025,30 +1042,30 @@ struct TaskPi0FlowEMC {
10251042
float dPhi = v1.Phi() - v2.Phi();
10261043
float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P()));
10271044

1028-
registry.fill(HIST("hClusterCuts"), 1);
1045+
registry.fill(HIST("hClusterCutsMixed"), 1);
10291046
if (openingAngle <= mesonConfig.minOpenAngle) {
1030-
registry.fill(HIST("hClusterCuts"), 2);
1047+
registry.fill(HIST("hClusterCutsMixed"), 2);
10311048
continue;
10321049
}
10331050
if (thnConfigAxisInvMass.value[1] > vMeson.M() || thnConfigAxisInvMass.value.back() < vMeson.M()) {
1034-
registry.fill(HIST("hClusterCuts"), 3);
1051+
registry.fill(HIST("hClusterCutsMixed"), 3);
10351052
continue;
10361053
}
10371054
if (thnConfigAxisPt.value[1] > vMeson.Pt() || thnConfigAxisPt.value.back() < vMeson.Pt()) {
1038-
registry.fill(HIST("hClusterCuts"), 4);
1055+
registry.fill(HIST("hClusterCutsMixed"), 4);
10391056
continue;
10401057
}
10411058
if (mesonConfig.cfgEnableQA) {
1042-
registry.fill(HIST("hInvMassPt"), vMeson.M(), vMeson.Pt());
1043-
registry.fill(HIST("hTanThetaPhi"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi)));
1044-
registry.fill(HIST("hAlphaPt"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt());
1059+
registry.fill(HIST("hInvMassPtMixed"), vMeson.M(), vMeson.Pt());
1060+
registry.fill(HIST("hTanThetaPhiMixed"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi)));
1061+
registry.fill(HIST("hAlphaPtMixed"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt());
10451062
}
10461063
if (mesonConfig.enableTanThetadPhi && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) {
1047-
registry.fill(HIST("hClusterCuts"), 5);
1064+
registry.fill(HIST("hClusterCutsMixed"), 5);
10481065
continue;
10491066
}
1050-
registry.fill(HIST("hClusterCuts"), 6);
1051-
runFlowAnalysis<1>(c1, vMeson);
1067+
registry.fill(HIST("hClusterCutsMixed"), 6);
1068+
runFlowAnalysis<2>(c1, vMeson);
10521069
}
10531070
}
10541071
}
@@ -1283,6 +1300,95 @@ struct TaskPi0FlowEMC {
12831300
}
12841301
PROCESS_SWITCH(TaskPi0FlowEMC, processEMCalCalib, "Process EMCal calibration", false);
12851302

1303+
// Pi0 from EMCal
1304+
void processM02(CollsWithQvecs const& collisions, EMCalPhotons const& clusters)
1305+
{
1306+
for (const auto& collision : collisions) {
1307+
auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex());
1308+
1309+
if (eventcuts.cfgEnableQA) {
1310+
registry.fill(HIST("hCollisionEMCCheck"), 1.); // all
1311+
if (collision.alias_bit(kTVXinEMC) == true) {
1312+
registry.fill(HIST("hCollisionEMCCheck"), 2.); // has EMC read out
1313+
if (photonsPerCollision.size() > 0) {
1314+
registry.fill(HIST("hCollisionEMCCheck"), 3.); // has EMC cluster
1315+
registry.fill(HIST("hCollisionEMCCheck"), 4.); // has EMC read out and clusters
1316+
} else {
1317+
registry.fill(HIST("hCollisionEMCCheck"), 5.); // has EMC read out but no clusters
1318+
}
1319+
} else {
1320+
if (photonsPerCollision.size() > 0) {
1321+
registry.fill(HIST("hCollisionEMCCheck"), 3.); // has EMC cluster
1322+
registry.fill(HIST("hCollisionEMCCheck"), 6.); // has no EMC read out and clusters
1323+
} else {
1324+
registry.fill(HIST("hCollisionEMCCheck"), 7.); // has no EMC read out and no clusters
1325+
}
1326+
}
1327+
}
1328+
o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(&registry, collision);
1329+
if (!(fEMEventCut.IsSelected(collision))) {
1330+
// general event selection
1331+
continue;
1332+
}
1333+
if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) {
1334+
// occupancy selection
1335+
continue;
1336+
}
1337+
float cent = getCentrality(collision);
1338+
if (cent < eventcuts.cfgMinCent || cent > eventcuts.cfgMaxCent) {
1339+
// event selection
1340+
continue;
1341+
}
1342+
if (!isQvecGood(getAllQvec(collision))) {
1343+
// selection based on QVector
1344+
continue;
1345+
}
1346+
runNow = collision.runNumber();
1347+
if (runNow != runBefore) {
1348+
initCCDB(collision);
1349+
runBefore = runNow;
1350+
}
1351+
o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(&registry, collision);
1352+
registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted
1353+
registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted
1354+
1355+
for (const auto& photon : photonsPerCollision) {
1356+
if (mesonConfig.cfgEnableQA) {
1357+
registry.fill(HIST("hEClusterBefore"), photon.e()); // before cuts
1358+
registry.fill(HIST("mesonQA/hClusterEtaPhiBefore"), photon.phi(), photon.eta()); // before cuts
1359+
}
1360+
if (!(fEMCCut.IsSelected<EMCalPhotons::iterator>(photon))) {
1361+
continue;
1362+
}
1363+
if (cfgDistanceToEdge.value && (checkEtaPhi1D(photon.eta(), RecoDecay::constrainAngle(photon.phi())) >= cfgEMCalMapLevelSameEvent.value)) {
1364+
continue;
1365+
}
1366+
if (mesonConfig.cfgEnableQA) {
1367+
registry.fill(HIST("hEClusterAfter"), photon.e()); // accepted after cuts
1368+
registry.fill(HIST("mesonQA/hClusterEtaPhiAfter"), photon.phi(), photon.eta()); // before cuts
1369+
}
1370+
1371+
auto [xQVec, yQVec] = getQvec(collision, qvecDetector);
1372+
float cent = getCentrality(collision);
1373+
1374+
float phiCand = photon.phi();
1375+
1376+
float cosNPhi = std::cos(harmonic * phiCand);
1377+
float sinNPhi = std::sin(harmonic * phiCand);
1378+
float scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec;
1379+
1380+
if (correctionConfig.cfgApplySPresolution.value) {
1381+
scalprodCand = scalprodCand / h1SPResolution->GetBinContent(h1SPResolution->FindBin(cent + epsilon));
1382+
}
1383+
if (cfgDoM02.value) {
1384+
registry.fill(HIST("hSparseFlow"), photon.m02(), photon.pt(), cent, scalprodCand);
1385+
}
1386+
return;
1387+
} // end of loop over single cluster
1388+
} // end of loop over collisions
1389+
} // processM02
1390+
PROCESS_SWITCH(TaskPi0FlowEMC, processM02, "Process single EMCal clusters as function of M02", false);
1391+
12861392
}; // End struct TaskPi0FlowEMC
12871393

12881394
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)