Skip to content

Commit 401d3ba

Browse files
omvazqueomvazque
authored
[PWGLF] Updated the MC correction procedure (#11590)
1 parent 517903c commit 401d3ba

File tree

1 file changed

+124
-67
lines changed

1 file changed

+124
-67
lines changed

‎PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx‎

Lines changed: 124 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -15,25 +15,14 @@
1515
/// \author Omar Vazquez (omar.vazquez.rueda@cern.ch)
1616
/// \since January 29, 2025
1717

18-
#include <CCDB/BasicCCDBManager.h>
19-
20-
#include <chrono>
21-
#include <cmath>
22-
#include <cstddef>
23-
#include <cstdint>
24-
#include <cstdlib>
25-
#include <string>
26-
#include <string_view>
27-
#include <vector>
28-
#include <TRandom.h>
29-
3018
#include "Common/CCDB/EventSelectionParams.h"
3119
#include "Common/CCDB/TriggerAliases.h"
3220
#include "Common/Core/TrackSelection.h"
3321
#include "Common/DataModel/Centrality.h"
3422
#include "Common/DataModel/EventSelection.h"
3523
#include "Common/DataModel/Multiplicity.h"
3624
#include "Common/DataModel/TrackSelectionTables.h"
25+
3726
#include "CommonConstants/MathConstants.h"
3827
#include "CommonConstants/ZDCConstants.h"
3928
#include "Framework/ASoAHelpers.h" // required for Filter op.
@@ -44,7 +33,20 @@
4433
#include "Framework/runDataProcessing.h"
4534
#include "ReconstructionDataFormats/GlobalTrackID.h"
4635
#include "ReconstructionDataFormats/Track.h"
36+
#include <CCDB/BasicCCDBManager.h>
37+
4738
#include "TPDGCode.h"
39+
#include <TRandom.h>
40+
41+
#include <chrono>
42+
#include <cmath>
43+
#include <cstddef>
44+
#include <cstdint>
45+
#include <cstdlib>
46+
#include <numeric>
47+
#include <string>
48+
#include <string_view>
49+
#include <vector>
4850

4951
using namespace std;
5052
using namespace o2;
@@ -81,6 +83,7 @@ struct UccZdc {
8183
Configurable<bool> isZEMcut{"isZEMcut", true, "Use ZEM cut"};
8284
Configurable<bool> useMidRapNchSel{"useMidRapNchSel", true, "Use mid-rapidit Nch selection"};
8385
Configurable<bool> applyEff{"applyEff", true, "Apply track-by-track efficiency correction"};
86+
Configurable<bool> correctNch{"correctNch", true, "Correct also Nch"};
8487

8588
// Event selection
8689
Configurable<float> posZcut{"posZcut", +10.0, "z-vertex position cut"};
@@ -126,7 +129,8 @@ struct UccZdc {
126129
ConfigurableAxis binsCent{"binsCent", {VARIABLE_WIDTH, 0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.}, "T0C binning"};
127130

128131
// CCDB paths
129-
Configurable<std::string> paTH{"paTH", "Users/o/omvazque/TrackingEfficiency", "base path to the ccdb object"};
132+
Configurable<std::string> paTHEff{"paTHEff", "Users/o/omvazque/MCcorrection/perTimeStamp/TrackingEff", "base path to the ccdb object"};
133+
Configurable<std::string> paTHFD{"paTHFD", "Users/o/omvazque/MCcorrection/perTimeStamp/FeedDown", "base path to the ccdb object"};
130134
Configurable<std::string> paTHmeanNch{"paTHmeanNch", "Users/o/omvazque/FitMeanNch_9May2025", "base path to the ccdb object"};
131135
Configurable<std::string> paTHsigmaNch{"paTHsigmaNch", "Users/o/omvazque/FitSigmaNch_9May2025", "base path to the ccdb object"};
132136
Configurable<int64_t> ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"};
@@ -238,6 +242,9 @@ struct UccZdc {
238242
registry.add("NchVsThreeParCorr", "MC closure;#it{N}_{ch} (|#eta| < 0.8, Corrected);#LT[#it{p}_{T}^{(3)}]#GT", kTProfile, {{nBinsNch, minNch, maxNch}});
239243
registry.add("NchVsFourParCorr", "MC closure;#it{N}_{ch} (|#eta| < 0.8, Corrected);#LT[#it{p}_{T}^{(4)}]#GT", kTProfile, {{nBinsNch, minNch, maxNch}});
240244
// Corrections
245+
registry.add("NchRec", "Corrections;#it{N}_{ch} (|#eta| < 0.8);Entries;", kTH1F, {{nBinsNch, minNch, maxNch}});
246+
registry.add("NchTrue", "Corrections;#it{N}_{ch} (|#eta| < 0.8);Entries;", kTH1F, {{nBinsNch, minNch, maxNch}});
247+
241248
registry.add("zPosMC", "Filled at MC closure + Corrections;;Entries;", kTH1F, {axisZpos});
242249
registry.add("hEventCounterMC", "Event counter", kTH1F, {axisEvent});
243250
registry.add("nRecColvsCent", "", kTH2F, {{6, -0.5, 5.5}, {{axisCent}}});
@@ -292,8 +299,11 @@ struct UccZdc {
292299

293300
LOG(info) << "\tccdbNoLaterThan=" << ccdbNoLaterThan.value;
294301
LOG(info) << "\tapplyEff=" << applyEff.value;
295-
LOG(info) << "\tpaTH=" << paTH.value;
302+
LOG(info) << "\tcorrectNch=" << correctNch.value;
303+
LOG(info) << "\tpaTHEff=" << paTHEff.value;
304+
LOG(info) << "\tpaTHFD=" << paTHFD.value;
296305
LOG(info) << "\tuseMidRapNchSel=" << useMidRapNchSel.value;
306+
LOG(info) << "\tnSigmaNchCut=" << nSigmaNchCut.value;
297307
LOG(info) << "\tpaTHmeanNch=" << paTHmeanNch.value;
298308
LOG(info) << "\tpaTHsigmaNch=" << paTHsigmaNch.value;
299309
LOG(info) << "\tminPt=" << minPt.value;
@@ -681,48 +691,55 @@ struct UccZdc {
681691
}
682692
}
683693

694+
// Skip event based on number of Nch sigmas
684695
if (!skipEvent) {
685696
return;
686697
}
687698

688-
auto efficiency = ccdb->getForTimeStamp<TH1F>(paTH.value, foundBC.timestamp());
689-
// auto efficiency = ccdb->getForRun<TH1F>(paTH.value, foundBC.runNumber());
690-
if (!efficiency) {
699+
auto efficiency = ccdb->getForTimeStamp<TH1F>(paTHEff.value, foundBC.timestamp());
700+
auto fd = ccdb->getForTimeStamp<TH1F>(paTHFD.value, foundBC.timestamp());
701+
if (!efficiency || !fd) {
691702
return;
692703
}
693704

694705
std::vector<float> pTs;
695-
std::vector<float> wIs;
696-
// Calculates the event weight, W_k
706+
std::vector<float> vecFD;
707+
std::vector<float> vecOneOverEff;
708+
709+
// Calculates the Nch multiplicity
697710
for (const auto& track : tracks) {
698711
// Track Selection
699712
if (!track.isGlobalTrack()) {
700713
continue;
701714
}
702-
if ((track.pt() < minPt) || (track.pt() > maxPtSpectra)) {
715+
if ((track.pt() < minPt) || (track.pt() > maxPt)) {
703716
continue;
704717
}
705718

706719
float pt{track.pt()};
707-
double weight{1.};
720+
float effValue{1.0};
708721
if (applyEff) {
709-
weight = efficiency->GetBinContent(efficiency->FindBin(pt));
722+
effValue = efficiency->GetBinContent(efficiency->FindBin(pt));
710723
}
711-
if (weight > 0.) {
712-
pTs.emplace_back(pt);
713-
wIs.emplace_back(weight);
724+
if (effValue > 0.) {
725+
vecOneOverEff.emplace_back(1. / effValue);
714726
}
715727
}
716728

717-
double p1, p2, p3, p4, w1, w2, w3, w4;
718-
p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0;
719-
getPTpowers(pTs, wIs, p1, w1, p2, w2, p3, w3, p4, w4);
720-
const double nch{static_cast<double>(pTs.size())};
721-
if (nch < minNchSel) {
729+
double nchMult{0.};
730+
nchMult = std::accumulate(vecOneOverEff.begin(), vecOneOverEff.end(), 0);
731+
if (!applyEff)
732+
nchMult = static_cast<double>(glbTracks);
733+
if (applyEff && !correctNch)
734+
nchMult = static_cast<double>(glbTracks);
735+
if (nchMult < minNchSel) {
722736
return;
723737
}
724738

725-
// To calculate event-averaged <pt>
739+
// Fill vectors for [pT] measurement
740+
pTs.clear();
741+
vecFD.clear();
742+
vecOneOverEff.clear();
726743
for (const auto& track : tracks) {
727744
// Track Selection
728745
if (!track.isGlobalTrack()) {
@@ -731,9 +748,27 @@ struct UccZdc {
731748
if ((track.pt() < minPt) || (track.pt() > maxPtSpectra)) {
732749
continue;
733750
}
734-
registry.fill(HIST("NchVsZNVsPt"), w1, sumZNs, track.pt());
751+
752+
float pt{track.pt()};
753+
float effValue{1.};
754+
float fdValue{1.};
755+
if (applyEff) {
756+
effValue = efficiency->GetBinContent(efficiency->FindBin(pt));
757+
fdValue = fd->GetBinContent(fd->FindBin(pt));
758+
}
759+
if ((effValue > 0.) && (fdValue > 0.)) {
760+
pTs.emplace_back(pt);
761+
vecOneOverEff.emplace_back(1. / effValue);
762+
vecFD.emplace_back(fdValue);
763+
}
764+
// To calculate event-averaged <pt>
765+
registry.fill(HIST("NchVsZNVsPt"), nchMult, sumZNs, track.pt());
735766
}
736767

768+
double p1, p2, p3, p4, w1, w2, w3, w4;
769+
p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0;
770+
getPTpowers(pTs, vecOneOverEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4);
771+
737772
// EbE one-particle pT correlation
738773
double oneParCorr{p1 / w1};
739774

@@ -752,20 +787,20 @@ struct UccZdc {
752787
double numFourParCorr{std::pow(p1, 4.) - 6. * p2 * std::pow(p1, 2.) + 3. * std::pow(p2, 2.) + 8 * p3 * p1 - 6. * p4};
753788
double fourParCorr{numFourParCorr / denFourParCorr};
754789

755-
registry.fill(HIST("Nch"), w1);
790+
registry.fill(HIST("Nch"), nchMult);
756791
registry.fill(HIST("ZNamp"), sumZNs);
757-
registry.fill(HIST("NchVsZN"), w1, sumZNs);
758-
registry.fill(HIST("NchVsZP"), w1, sumZPs);
792+
registry.fill(HIST("NchVsZN"), nchMult, sumZNs);
793+
registry.fill(HIST("NchVsZP"), nchMult, sumZPs);
759794
registry.fill(HIST("NITSTacksVsZN"), itsTracks, sumZNs);
760795
registry.fill(HIST("NITSTacksVsZP"), itsTracks, sumZPs);
761796
registry.fill(HIST("T0MVsZN"), normT0M, sumZNs);
762797
registry.fill(HIST("T0MVsZP"), normT0M, sumZPs);
763798
registry.fill(HIST("NchUncorrected"), glbTracks);
764-
registry.fill(HIST("NchVsOneParCorr"), w1, oneParCorr, w1);
765-
registry.fill(HIST("NchVsOneParCorrVsZN"), w1, sumZNs, oneParCorr, w1);
766-
registry.fill(HIST("NchVsTwoParCorrVsZN"), w1, sumZNs, twoParCorr, denTwoParCorr);
767-
registry.fill(HIST("NchVsThreeParCorrVsZN"), w1, sumZNs, threeParCorr, denThreeParCorr);
768-
registry.fill(HIST("NchVsFourParCorrVsZN"), w1, sumZNs, fourParCorr, denFourParCorr);
799+
registry.fill(HIST("NchVsOneParCorr"), nchMult, oneParCorr, w1);
800+
registry.fill(HIST("NchVsOneParCorrVsZN"), nchMult, sumZNs, oneParCorr, w1);
801+
registry.fill(HIST("NchVsTwoParCorrVsZN"), nchMult, sumZNs, twoParCorr, denTwoParCorr);
802+
registry.fill(HIST("NchVsThreeParCorrVsZN"), nchMult, sumZNs, threeParCorr, denThreeParCorr);
803+
registry.fill(HIST("NchVsFourParCorrVsZN"), nchMult, sumZNs, fourParCorr, denFourParCorr);
769804
}
770805
PROCESS_SWITCH(UccZdc, processZdcCollAss, "Process ZDC W/Coll Ass.", true);
771806

@@ -774,7 +809,6 @@ struct UccZdc {
774809
TRandom* randPointer = new TRandom();
775810
void processMCclosure(aod::McCollisions::iterator const& mccollision, soa::SmallGroups<o2::aod::SimCollisions> const& collisions, o2::aod::BCsRun3 const& /*bcs*/, aod::McParticles const& mcParticles, TheFilteredSimTracks const& simTracks)
776811
{
777-
778812
float rndNum = randPointer->Uniform(0.0, 1.0);
779813
registry.fill(HIST("RandomNumber"), rndNum);
780814

@@ -809,14 +843,17 @@ struct UccZdc {
809843

810844
// To use run-by-run efficiency
811845
const auto& foundBC = collision.foundBC_as<o2::aod::BCsRun3>();
812-
// auto efficiency = ccdb->getForTimeStamp<TH1F>(paTH.value, foundBC.timestamp());
813-
auto efficiency = ccdb->getForRun<TH1F>(paTH.value, foundBC.runNumber());
814-
if (!efficiency) {
815-
continue;
846+
auto efficiency = ccdb->getForTimeStamp<TH1F>(paTHEff.value, foundBC.timestamp());
847+
auto fd = ccdb->getForTimeStamp<TH1F>(paTHFD.value, foundBC.timestamp());
848+
if (!efficiency || !fd) {
849+
return;
816850
}
817851

852+
int nchRaw{0};
818853
std::vector<float> pTs;
819-
std::vector<float> wIs;
854+
std::vector<float> vecFD;
855+
std::vector<float> vecOneOverEff;
856+
// std::vector<float> wIs;
820857
const auto& groupedTracks{simTracks.sliceBy(perCollision, collision.globalIndex())};
821858
// Calculates the event weight, W_k
822859
for (const auto& track : groupedTracks) {
@@ -826,22 +863,29 @@ struct UccZdc {
826863
}
827864

828865
float pt{track.pt()};
829-
double weight{efficiency->GetBinContent(efficiency->FindBin(pt))};
830-
if (!(weight > 0.)) {
831-
continue;
866+
float effValue{1.};
867+
float fdValue{1.};
868+
nchRaw++;
869+
if (applyEff) {
870+
effValue = efficiency->GetBinContent(efficiency->FindBin(pt));
871+
fdValue = fd->GetBinContent(fd->FindBin(pt));
872+
}
873+
if ((effValue > 0.) && (fdValue > 0.)) {
874+
pTs.emplace_back(pt);
875+
vecOneOverEff.emplace_back(1. / effValue);
876+
vecFD.emplace_back(fdValue);
832877
}
833-
pTs.emplace_back(pt);
834-
wIs.emplace_back(weight);
835878
}
836879

837-
const double nch{static_cast<double>(pTs.size())};
838-
if (nch < minNchSel) {
839-
continue;
880+
double nchMult{0.};
881+
nchMult = std::accumulate(vecOneOverEff.begin(), vecOneOverEff.end(), 0);
882+
if (nchMult < minNchSel) {
883+
return;
840884
}
841885

842886
double p1, p2, p3, p4, w1, w2, w3, w4;
843887
p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0;
844-
getPTpowers(pTs, wIs, p1, w1, p2, w2, p3, w3, p4, w4);
888+
getPTpowers(pTs, vecOneOverEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4);
845889

846890
const double denTwoParCorr{std::pow(w1, 2.) - w2};
847891
const double numTwoParCorr{std::pow(p1, 2.) - p2};
@@ -855,16 +899,18 @@ struct UccZdc {
855899
const double threeParCorr{numThreeParCorr / denThreeParCorr};
856900
const double fourParCorr{numFourParCorr / denFourParCorr};
857901

858-
registry.fill(HIST("Nch"), w1);
859-
registry.fill(HIST("NchUncorrected"), nch);
860-
registry.fill(HIST("NchVsOneParCorr"), w1, oneParCorr, w1);
861-
registry.fill(HIST("NchVsTwoParCorr"), w1, twoParCorr, denTwoParCorr);
862-
registry.fill(HIST("NchVsThreeParCorr"), w1, threeParCorr, denThreeParCorr);
863-
registry.fill(HIST("NchVsFourParCorr"), w1, fourParCorr, denFourParCorr);
902+
registry.fill(HIST("Nch"), nchMult);
903+
registry.fill(HIST("NchUncorrected"), nchRaw);
904+
registry.fill(HIST("NchVsOneParCorr"), nchMult, oneParCorr, w1);
905+
registry.fill(HIST("NchVsTwoParCorr"), nchMult, twoParCorr, denTwoParCorr);
906+
registry.fill(HIST("NchVsThreeParCorr"), nchMult, threeParCorr, denThreeParCorr);
907+
registry.fill(HIST("NchVsFourParCorr"), nchMult, fourParCorr, denFourParCorr);
864908

865909
//--------------------------- Generated MC ---------------------------
866910
std::vector<float> pTsMC;
867-
std::vector<float> wIsMC;
911+
std::vector<float> vecFullEff;
912+
std::vector<float> vecFDEqualOne;
913+
868914
// Calculates the event weight, W_k
869915
for (const auto& particle : mcParticles) {
870916
if (particle.eta() < minEta || particle.eta() > maxEta) {
@@ -879,17 +925,19 @@ struct UccZdc {
879925

880926
float pt{particle.pt()};
881927
pTsMC.emplace_back(pt);
882-
wIsMC.emplace_back(1.);
928+
vecFullEff.emplace_back(1.);
929+
vecFDEqualOne.emplace_back(1.);
883930
}
884931

885-
const double nchMC{static_cast<double>(pTsMC.size())};
932+
double nchMC{0};
933+
nchMult = std::accumulate(vecFullEff.begin(), vecFullEff.end(), 0);
886934
if (nchMC < minNchSel) {
887935
continue;
888936
}
889937

890938
double p1MC, p2MC, p3MC, p4MC, w1MC, w2MC, w3MC, w4MC;
891939
p1MC = p2MC = p3MC = p4MC = w1MC = w2MC = w3MC = w4MC = 0.0;
892-
getPTpowers(pTsMC, wIsMC, p1MC, w1MC, p2MC, w2MC, p3MC, w3MC, p4MC, w4MC);
940+
getPTpowers(pTsMC, vecFullEff, vecFDEqualOne, p1MC, w1MC, p2MC, w2MC, p3MC, w3MC, p4MC, w4MC);
893941

894942
const double denTwoParCorrMC{std::pow(w1MC, 2.) - w2MC};
895943
const double numTwoParCorrMC{std::pow(p1MC, 2.) - p2MC};
@@ -911,6 +959,8 @@ struct UccZdc {
911959
} else { // Correction with the remaining half of the sample
912960
registry.fill(HIST("EvtsDivided"), 1);
913961
//----- MC reconstructed -----//
962+
int nchTrue{0};
963+
int nchRec{0};
914964
const auto& groupedTracks{simTracks.sliceBy(perCollision, collision.globalIndex())};
915965
for (const auto& track : groupedTracks) {
916966
// Track Selection
@@ -932,6 +982,7 @@ struct UccZdc {
932982
continue;
933983
}
934984

985+
nchRec++;
935986
registry.fill(HIST("Pt_ch"), cent, track.pt());
936987
if (particle.pdgCode() == PDG_t::kPiPlus || particle.pdgCode() == PDG_t::kPiMinus) {
937988
registry.fill(HIST("Pt_pi"), cent, track.pt());
@@ -960,6 +1011,7 @@ struct UccZdc {
9601011
continue;
9611012
}
9621013

1014+
nchTrue++;
9631015
registry.fill(HIST("PtMC_ch"), cent, particle.pt());
9641016
if (particle.pdgCode() == PDG_t::kPiPlus || particle.pdgCode() == PDG_t::kPiMinus) { // pion
9651017
registry.fill(HIST("PtMC_pi"), cent, particle.pt());
@@ -975,18 +1027,23 @@ struct UccZdc {
9751027
registry.fill(HIST("PtMC_re"), cent, particle.pt());
9761028
}
9771029
}
1030+
1031+
registry.fill(HIST("NchRec"), nchRec);
1032+
registry.fill(HIST("NchTrue"), nchTrue);
9781033
} // Half of statistics for corrections
9791034
} // Collisions
9801035
}
9811036
PROCESS_SWITCH(UccZdc, processMCclosure, "Process MC closure", false);
9821037

9831038
template <typename T, typename U>
984-
void getPTpowers(const T& pTs, const T& wIs, U& pOne, U& wOne, U& pTwo, U& wTwo, U& pThree, U& wThree, U& pFour, U& wFour)
1039+
void getPTpowers(const T& pTs, const T& vecOneOverEff, const T& vecFD, U& pOne, U& wOne, U& pTwo, U& wTwo, U& pThree, U& wThree, U& pFour, U& wFour)
9851040
{
9861041
pOne = wOne = pTwo = wTwo = pThree = wThree = pFour = wFour = 0.;
9871042
for (std::size_t i = 0; i < pTs.size(); ++i) {
9881043
const float pTi{pTs.at(i)};
989-
const float wEighti{wIs.at(i)};
1044+
const float eFFi{vecOneOverEff.at(i)};
1045+
const float fDi{vecFD.at(i)};
1046+
const float wEighti{eFFi * fDi};
9901047
pOne += wEighti * pTi;
9911048
wOne += wEighti;
9921049
pTwo += std::pow(wEighti * pTi, 2.);

0 commit comments

Comments
 (0)