Skip to content

Commit bc5cea5

Browse files
peressounkoperessounko
peressounko
and
peressounko
authored
[Common,PWGEM] Fix data nonlinearity; add pi0 single histos (#9846)
Co-authored-by: peressounko <Dmitri.Peresunko@cern.ch>
1 parent 613cf98 commit bc5cea5

File tree

2 files changed

+207
-38
lines changed

2 files changed

+207
-38
lines changed

‎Common/TableProducer/caloClusterProducer.cxx‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1728,14 +1728,14 @@ struct CaloClusterProducer {
17281728
case 0:
17291729
return en;
17301730
case 1: { // Data Run3
1731-
const double a = 9.3494e-01;
1732-
const double b = 1.00526e-02;
1733-
const double c = 8.45164e-02;
1734-
const double d = -1.03364e-02;
1735-
const double f = 5.4803e-03;
1736-
const double g = 0.779983;
1737-
const double h = 0.622282;
1738-
const double k = 8.0182e-05;
1731+
const double a = 0.892787;
1732+
const double b = 0.004053;
1733+
const double c = 0.074652;
1734+
const double d = -0.016306;
1735+
const double f = 7.616314;
1736+
const double g = -104.409;
1737+
const double h = 1837.17;
1738+
const double k = 0.000091;
17391739
double eMin = std::max(static_cast<float>(0.1), en); // Parameterization valid down to 100 MeV
17401740
return en * (a + b * eMin + c / eMin + d / (eMin * eMin) + f / ((eMin - g) * (eMin - g) + h * h) + k / std::pow(eMin, 4));
17411741
}

‎PWGEM/Tasks/phosPi0.cxx‎

Lines changed: 199 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,10 @@ struct PhosPi0 {
6363
Configurable<int> nMixedEvents{"nMixedEvents", 10, "number of events to mix"};
6464
Configurable<bool> fillQC{"fillQC", true, "Fill QC histos"};
6565
Configurable<float> minOccE{"minOccE", 0.5, "Min. cluster energy of occupancy plots"};
66+
Configurable<float> nonlinA{"nonlinA", 1., "nonlinsrity param A (scale)"};
67+
Configurable<float> nonlinB{"nonlinB", 0., "nonlinsrity param B (a+b*exp(-e/c))"};
68+
Configurable<float> nonlinC{"nonlinC", 1., "nonlinsrity param C (a+b*exp(-e/c))"};
69+
Configurable<int> tofEffParam{"tofEffParam", 0, "parameterization of TOF cut efficiency"};
6670

6771
using SelCollisions = soa::Join<aod::Collisions, aod::EvSels>;
6872
using SelCollisionsMC = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels>;
@@ -81,20 +85,22 @@ struct PhosPi0 {
8185
{
8286
public:
8387
Photon() = default;
84-
Photon(double x, double y, double z, double ee, int m, bool isDispOK, bool isCPVOK, int mcLabel) : px(x), py(y), pz(z), e(ee), mod(m), mPID(isDispOK << 1 | isCPVOK << 2), label(mcLabel) {}
88+
Photon(double x, double y, double z, double ee, double t, int m, bool isDispOK, bool isCPVOK, int mcLabel) : px(x), py(y), pz(z), e(ee), time(t), mod(m), mPID(isDispOK << 1 | isCPVOK << 2), label(mcLabel) {}
8589
~Photon() = default;
8690

8791
bool isCPVOK() const { return (mPID >> 2) & 1; }
8892
bool isDispOK() const { return (mPID >> 1) & 1; }
93+
double pt() const { return std::sqrt(px * px + py * py); }
8994

9095
public:
91-
double px = 0.; // px
92-
double py = 0.; // py
93-
double pz = 0.; // pz
94-
double e = 0.; // energy
95-
int mod = 0; // module
96-
int mPID = 0; // store PID bits
97-
int label = -1; // label of MC particle
96+
double px = 0.; // px
97+
double py = 0.; // py
98+
double pz = 0.; // pz
99+
double e = 0.; // energy
100+
double time = 0.; // time
101+
int mod = 0; // module
102+
int mPID = 0; // store PID bits
103+
int label = -1; // label of MC particle
98104
};
99105

100106
int mRunNumber = 0; // Current run number
@@ -110,11 +116,13 @@ struct PhosPi0 {
110116
TH3 *hReMod, *hMiMod;
111117
TH2 *hReAll, *hReDisp, *hReCPV, *hReBoth, *hSignalAll, *hPi0SignalAll, *hPi0SignalCPV, *hPi0SignalDisp,
112118
*hPi0SignalBoth, *hMiAll, *hMiDisp, *hMiCPV, *hMiBoth;
119+
TH2 *hReOneAll, *hReOneDisp, *hReOneCPV, *hReOneBoth, *hMiOneAll, *hMiOneDisp, *hMiOneCPV, *hMiOneBoth;
120+
TH2 *hReTime12, *hReTime30, *hReTime50, *hReTime100;
113121

114-
std::vector<double> pt = {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2,
115-
1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0,
116-
6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 10., 11., 12., 13., 14., 15., 16., 18., 20., 22., 24., 26., 28.,
117-
30., 34., 38., 42., 46., 50., 55., 60., 70., 80., 90., 100., 110., 120., 150.};
122+
std::vector<double> pt = {0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1,
123+
1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.5, 4.6, 4.8, 5.0,
124+
5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 22., 24., 26., 28.,
125+
30., 34., 38., 42., 46., 50., 55., 60., 70., 75., 80., 85., 90., 95., 100., 110., 120., 130., 140., 150., 160., 180., 200.};
118126

119127
/// \brief Create output histograms
120128
void init(InitContext const&)
@@ -194,6 +202,31 @@ struct PhosPi0 {
194202
hReBoth = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReBoth", "inv mass for centrality",
195203
HistType::kTH2F, {mggAxis, ptAxis}))
196204
.get();
205+
hReOneAll = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReOneAll", "inv mass for centrality",
206+
HistType::kTH2F, {mggAxis, ptAxis}))
207+
.get();
208+
hReOneCPV = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReOneCPV", "inv mass for centrality",
209+
HistType::kTH2F, {mggAxis, ptAxis}))
210+
.get();
211+
hReOneDisp = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReOneDisp", "inv mass for centrality",
212+
HistType::kTH2F, {mggAxis, ptAxis}))
213+
.get();
214+
hReOneBoth = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReOneBoth", "inv mass for centrality",
215+
HistType::kTH2F, {mggAxis, ptAxis}))
216+
.get();
217+
218+
hReTime12 = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReTime12", "inv mass for centrality",
219+
HistType::kTH2F, {mggAxis, ptAxis}))
220+
.get();
221+
hReTime30 = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReTime30", "inv mass for centrality",
222+
HistType::kTH2F, {mggAxis, ptAxis}))
223+
.get();
224+
hReTime50 = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReTime50", "inv mass for centrality",
225+
HistType::kTH2F, {mggAxis, ptAxis}))
226+
.get();
227+
hReTime100 = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggReTime100", "inv mass for centrality",
228+
HistType::kTH2F, {mggAxis, ptAxis}))
229+
.get();
197230

198231
if (isMC) {
199232
hSignalAll = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggSignal", "inv mass for correlated pairs",
@@ -228,6 +261,18 @@ struct PhosPi0 {
228261
hMiBoth = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggMiBoth", "inv mass for centrality",
229262
HistType::kTH2F, {mggAxis, ptAxis}))
230263
.get();
264+
hMiOneAll = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggMiAll", "inv mass for centrality",
265+
HistType::kTH2F, {mggAxis, ptAxis}))
266+
.get();
267+
hMiOneCPV = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggMiCPV", "inv mass for centrality",
268+
HistType::kTH2F, {mggAxis, ptAxis}))
269+
.get();
270+
hMiOneDisp = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggMiDisp", "inv mass for centrality",
271+
HistType::kTH2F, {mggAxis, ptAxis}))
272+
.get();
273+
hMiOneBoth = std::get<std::shared_ptr<TH2>>(mHistManager.add("mggMiBoth", "inv mass for centrality",
274+
HistType::kTH2F, {mggAxis, ptAxis}))
275+
.get();
231276
if (isMC) {
232277
mHistManager.add("hMCPi0SpAll", "pi0 spectrum inclusive", HistType::kTH1F, {ptAxis});
233278
mHistManager.add("hMCPi0SpPrim", "pi0 spectrum Primary", HistType::kTH1F, {ptAxis});
@@ -395,7 +440,7 @@ struct PhosPi0 {
395440
if (mcPart->begin().mcCollisionId() != mPrevMCColId) {
396441
mPrevMCColId = mcPart->begin().mcCollisionId(); // to avoid scanning full MC table each BC
397442
for (const auto& part : *mcPart) {
398-
if (part.mcCollision().bcId() != cluMcBCId) {
443+
if (part.mcCollision().bcId() != col.bcId()) {
399444
continue;
400445
}
401446
if (part.pdgCode() == 111) {
@@ -462,7 +507,11 @@ struct PhosPi0 {
462507
mcLabel = mcList[0];
463508
}
464509
}
465-
Photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod(), testLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > cpvCut, mcLabel);
510+
double enCorr = 1;
511+
if constexpr (isMC) { // correct MC energy
512+
enCorr = nonlinearity(clu.e());
513+
}
514+
Photon ph1(clu.px() * enCorr, clu.py() * enCorr, clu.pz() * enCorr, clu.e() * enCorr, clu.time(), clu.mod(), testLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > cpvCut, mcLabel);
466515
// Mix with other photons added to stack
467516
for (const auto& ph2 : mCurEvent) {
468517
double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) -
@@ -473,35 +522,89 @@ struct PhosPi0 {
473522
double pt = std::sqrt(std::pow(ph1.px + ph2.px, 2) +
474523
std::pow(ph1.py + ph2.py, 2));
475524
int modComb = moduleCombination(ph1.mod, ph2.mod);
476-
hReMod->Fill(m, pt, modComb);
477-
hReAll->Fill(m, pt);
525+
double w = 1.;
526+
if constexpr (isMC) { // correct MC energy
527+
w = tofCutEff(ph1.e) * tofCutEff(ph2.e);
528+
}
529+
hReMod->Fill(m, pt, modComb, w);
530+
hReAll->Fill(m, pt, w);
531+
hReOneAll->Fill(m, ph1.pt(), w);
532+
hReOneAll->Fill(m, ph2.pt(), w);
533+
if (ph1.isCPVOK()) {
534+
hReOneCPV->Fill(m, ph1.pt(), w);
535+
}
536+
if (ph2.isCPVOK()) {
537+
hReOneCPV->Fill(m, ph2.pt(), w);
538+
}
539+
if (ph1.isDispOK()) {
540+
hReOneDisp->Fill(m, ph1.pt(), w);
541+
if (ph1.isCPVOK()) {
542+
hReOneBoth->Fill(m, ph1.pt(), w);
543+
}
544+
}
545+
if (ph2.isDispOK()) {
546+
hReOneDisp->Fill(m, ph2.pt(), w);
547+
if (ph2.isCPVOK()) {
548+
hReOneBoth->Fill(m, ph2.pt(), w);
549+
}
550+
}
551+
// Test time eff
552+
if (std::abs(ph1.time) < 12.5e-9) { // strict cut on first photon
553+
if (std::abs(ph2.time) < 100.e-9) {
554+
hReTime100->Fill(m, ph2.pt());
555+
if (std::abs(ph2.time) < 50.e-9) {
556+
hReTime50->Fill(m, ph2.pt());
557+
if (std::abs(ph2.time) < 30.e-9) {
558+
hReTime30->Fill(m, ph2.pt());
559+
if (std::abs(ph2.time) < 12.5e-9) {
560+
hReTime12->Fill(m, ph2.pt());
561+
}
562+
}
563+
}
564+
}
565+
}
566+
if (std::abs(ph2.time) < 12.5e-9) { // strict cut on first photon
567+
if (std::abs(ph1.time) < 100.e-9) {
568+
hReTime100->Fill(m, ph1.pt());
569+
if (std::abs(ph1.time) < 50.e-9) {
570+
hReTime50->Fill(m, ph1.pt());
571+
if (std::abs(ph1.time) < 30.e-9) {
572+
hReTime30->Fill(m, ph1.pt());
573+
if (std::abs(ph1.time) < 12.5e-9) {
574+
hReTime12->Fill(m, ph1.pt());
575+
}
576+
}
577+
}
578+
}
579+
}
580+
478581
bool isPi0 = false;
479582
if constexpr (isMC) { // test parent
480583
int cp = commonParentPDG(ph1.label, ph2.label, mcPart);
481584
if (cp != 0) {
482-
hSignalAll->Fill(m, pt);
585+
hSignalAll->Fill(m, pt, w);
483586
if (cp == 111) {
484587
isPi0 = true;
485-
hPi0SignalAll->Fill(m, pt);
588+
hPi0SignalAll->Fill(m, pt, w);
486589
}
487590
}
488591
}
489592

490593
if (ph1.isCPVOK() && ph2.isCPVOK()) {
491-
hReCPV->Fill(m, pt);
594+
hReCPV->Fill(m, pt, w);
492595
if (isPi0) {
493-
hPi0SignalCPV->Fill(m, pt);
596+
hPi0SignalCPV->Fill(m, pt, w);
494597
}
495598
}
496599
if (ph1.isDispOK() && ph2.isDispOK()) {
497-
hReDisp->Fill(m, pt);
600+
hReDisp->Fill(m, pt, w);
498601
if (isPi0) {
499-
hPi0SignalDisp->Fill(m, pt);
602+
hPi0SignalDisp->Fill(m, pt, w);
500603
}
501604
if (ph1.isCPVOK() && ph2.isCPVOK()) {
502-
hReBoth->Fill(m, pt);
605+
hReBoth->Fill(m, pt, w);
503606
if (isPi0) {
504-
hPi0SignalBoth->Fill(m, pt);
607+
hPi0SignalBoth->Fill(m, pt, w);
505608
}
506609
}
507610
}
@@ -523,15 +626,39 @@ struct PhosPi0 {
523626
double pt = std::sqrt(std::pow(ph1.px + ph2.px, 2) +
524627
std::pow(ph1.py + ph2.py, 2));
525628
int modComb = moduleCombination(ph1.mod, ph2.mod);
526-
hMiMod->Fill(m, pt, modComb);
527-
hMiAll->Fill(m, pt);
629+
double w = 1.;
630+
if constexpr (isMC) { // correct MC energy
631+
w = tofCutEff(ph1.e) * tofCutEff(ph2.e);
632+
}
633+
hMiMod->Fill(m, pt, modComb, w);
634+
hMiAll->Fill(m, pt, w);
635+
hMiOneAll->Fill(m, ph1.pt(), w);
636+
hMiOneAll->Fill(m, ph2.pt(), w);
637+
if (ph1.isCPVOK()) {
638+
hMiOneCPV->Fill(m, ph1.pt(), w);
639+
}
640+
if (ph2.isCPVOK()) {
641+
hMiOneCPV->Fill(m, ph2.pt(), w);
642+
}
643+
if (ph1.isDispOK()) {
644+
hMiOneDisp->Fill(m, ph1.pt(), w);
645+
if (ph1.isCPVOK()) {
646+
hMiOneBoth->Fill(m, ph1.pt(), w);
647+
}
648+
}
649+
if (ph2.isDispOK()) {
650+
hMiOneDisp->Fill(m, ph2.pt(), w);
651+
if (ph2.isCPVOK()) {
652+
hMiOneBoth->Fill(m, ph2.pt(), w);
653+
}
654+
}
528655
if (ph1.isCPVOK() && ph2.isCPVOK()) {
529-
hMiCPV->Fill(m, pt);
656+
hMiCPV->Fill(m, pt, w);
530657
}
531658
if (ph1.isDispOK() && ph2.isDispOK()) {
532-
hMiDisp->Fill(m, pt);
659+
hMiDisp->Fill(m, pt, w);
533660
if (ph1.isCPVOK() && ph2.isCPVOK()) {
534-
hMiBoth->Fill(m, pt);
661+
hMiBoth->Fill(m, pt, w);
535662
}
536663
}
537664
}
@@ -597,7 +724,12 @@ struct PhosPi0 {
597724
}
598725

599726
int mcLabel = -1;
600-
Photon ph1(clu.px(), clu.py(), clu.pz(), clu.e(), clu.mod(), testLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > cpvCut, mcLabel);
727+
double enCorr = 1;
728+
if (isMC) { // correct MC energy
729+
enCorr = nonlinearity(clu.e());
730+
}
731+
Photon ph1(clu.px() * enCorr, clu.py() * enCorr, clu.pz() * enCorr, clu.e() * enCorr, clu.time(), clu.mod(), testLambda(clu.e(), clu.m02(), clu.m20()), clu.trackdist() > cpvCut, mcLabel);
732+
601733
// Mix with other photons added to stack
602734
for (const auto& ph2 : mCurEvent) {
603735
double m = std::pow(ph1.e + ph2.e, 2) - std::pow(ph1.px + ph2.px, 2) -
@@ -737,6 +869,43 @@ struct PhosPi0 {
737869
}
738870
return 0; // nothing found
739871
}
872+
double nonlinearity(double e)
873+
{
874+
return nonlinA + nonlinB * std::exp(-e / nonlinC);
875+
}
876+
double tofCutEff(double en)
877+
{
878+
if (tofEffParam == 0) {
879+
return 1.;
880+
}
881+
if (tofEffParam == 1) { // Run2 100 ns
882+
// parameterization 01.08.2020
883+
if (en > 1.1)
884+
en = 1.1;
885+
if (en < 0.11)
886+
en = 0.11;
887+
return std::exp((-1.15295e+05 + 2.26754e+05 * en - 1.26063e+05 * en * en + en * en * en) /
888+
(1. - 3.16443e+05 * en + 3.68044e+06 * en * en + en * en * en));
889+
}
890+
if (tofEffParam == 2) { // Run2 30 ns
891+
if (en > 1.6)
892+
en = 1.6;
893+
return 1. / (1. + std::exp((4.83230e+01 - 8.89758e+01 * en + 1.10897e+03 * en * en - 5.73755e+03 * en * en * en -
894+
1.43777e+03 * en * en * en * en) /
895+
(1. - 1.23667e+02 * en + 1.07255e+03 * en * en + 5.87221e+02 * en * en * en)));
896+
}
897+
if (tofEffParam == 2) { // Run2 12.5 ns
898+
if (en < 4.6) {
899+
return std::exp(3.64952e-03 *
900+
(-5.80032e+01 - 1.53442e+02 * en + 1.30994e+02 * en * en + -3.53094e+01 * en * en * en + en * en * en * en) /
901+
(-7.75638e-02 + 8.64761e-01 * en + 1.22320e+00 * en * en - 1.00177e+00 * en * en * en + en * en * en * en));
902+
} else {
903+
return 0.63922783 * (1. - 1.63273e-01 * std::tanh((en - 7.94528e+00) / 1.28997e+00)) *
904+
(-4.39257e+00 * en + 2.25503e+00 * en * en + en * en * en) / (2.37160e+00 * en - 6.93786e-01 * en * en + en * en * en);
905+
}
906+
}
907+
return 1.;
908+
}
740909
};
741910

742911
o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc)

0 commit comments

Comments
 (0)