Skip to content

Commit 05d0bb2

Browse files
omvazqueomvazque
authored
[PWGLF] Adds process function for Poisson bootstrap (#11929)
1 parent 8a5f6ae commit 05d0bb2

File tree

1 file changed

+267
-0
lines changed

1 file changed

+267
-0
lines changed

‎PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx‎

Lines changed: 267 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,6 +37,7 @@
3737

3838
#include "TPDGCode.h"
3939
#include <TRandom.h>
40+
#include <TRandom3.h>
4041

4142
#include <chrono>
4243
#include <cmath>
@@ -46,6 +47,7 @@
4647
#include <numeric>
4748
#include <string>
4849
#include <string_view>
50+
#include <typeinfo>
4951
#include <vector>
5052

5153
using namespace std;
@@ -65,6 +67,16 @@ using SimCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLa
6567
using SimTracks = soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::McTrackLabels>;
6668
} // namespace o2::aod
6769

70+
static constexpr int kSizeBootStrapEnsemble{8};
71+
72+
std::array<std::shared_ptr<TH1>, kSizeBootStrapEnsemble> hPoisson{};
73+
std::array<std::shared_ptr<TH1>, kSizeBootStrapEnsemble> hNch{};
74+
std::array<std::shared_ptr<TProfile>, kSizeBootStrapEnsemble> pNchVsOneParCorr{};
75+
std::array<std::shared_ptr<TProfile2D>, kSizeBootStrapEnsemble> pNchVsOneParCorrVsZN{};
76+
std::array<std::shared_ptr<TProfile2D>, kSizeBootStrapEnsemble> pNchVsTwoParCorrVsZN{};
77+
std::array<std::shared_ptr<TProfile2D>, kSizeBootStrapEnsemble> pNchVsThreeParCorrVsZN{};
78+
std::array<std::shared_ptr<TProfile2D>, kSizeBootStrapEnsemble> pNchVsFourParCorrVsZN{};
79+
6880
struct UccZdc {
6981

7082
static constexpr float kCollEnergy{2.68};
@@ -229,6 +241,18 @@ struct UccZdc {
229241
registry.add("NchVsFourParCorrVsZN", ";#it{N}_{ch} (|#eta| < 0.8, Corrected);ZNA+ZNC;#LT[#it{p}_{T}^{(4)}]#GT", kTProfile2D, {{{nBinsNch, minNch, maxNch}, {nBinsZN, -0.5, maxZN}}});
230242
}
231243

244+
if (doprocessEventSampling) {
245+
for (int i = 0; i < kSizeBootStrapEnsemble; i++) {
246+
hNch[i] = registry.add<TH1>(Form("Nch_Replica%d", i), ";#it{N}_{ch} (|#eta| < 0.8);Entries", kTH1F, {{nBinsNch, minNch, maxNch}});
247+
hPoisson[i] = registry.add<TH1>(Form("Poisson_Replica%d", i), ";#it{k};Entries", kTH1F, {{21, -0.5, 20.5}});
248+
pNchVsOneParCorrVsZN[i] = registry.add<TProfile2D>(Form("NchVsOneParCorrVsZN_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; ZNA+ZNC; #LT[#it{p}_{T}^{(1)}]#GT (GeV/#it{c})", kTProfile2D, {{{nBinsNch, minNch, maxNch}, {nBinsZN, -0.5, maxZN}}});
249+
pNchVsTwoParCorrVsZN[i] = registry.add<TProfile2D>(Form("NchVsTwoParCorrVsZN_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; ZNA+ZNC; #LT[#it{p}_{T}^{(1)}]#GT (GeV/#it{c})", kTProfile2D, {{{nBinsNch, minNch, maxNch}, {nBinsZN, -0.5, maxZN}}});
250+
pNchVsThreeParCorrVsZN[i] = registry.add<TProfile2D>(Form("NchVsThreeParCorrVsZN_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; ZNA+ZNC; #LT[#it{p}_{T}^{(1)}]#GT (GeV/#it{c})", kTProfile2D, {{{nBinsNch, minNch, maxNch}, {nBinsZN, -0.5, maxZN}}});
251+
pNchVsFourParCorrVsZN[i] = registry.add<TProfile2D>(Form("NchVsFourParCorrVsZN_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; ZNA+ZNC; #LT[#it{p}_{T}^{(1)}]#GT (GeV/#it{c})", kTProfile2D, {{{nBinsNch, minNch, maxNch}, {nBinsZN, -0.5, maxZN}}});
252+
pNchVsOneParCorr[i] = registry.add<TProfile>(Form("NchVsOneParCorr_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8;#LT[#it{p}_{T}^{(1)}]#GT (GeV/#it{c})", kTProfile, {{nBinsNch, minNch, maxNch}});
253+
}
254+
}
255+
232256
if (doprocessMCclosure) {
233257
registry.add("RandomNumber", "", kTH1F, {{50, 0., 1.}});
234258
registry.add("EvtsDivided", ";Event type;Entries;", kTH1F, {{2, -0.5, 1.5}});
@@ -826,6 +850,249 @@ struct UccZdc {
826850
}
827851
PROCESS_SWITCH(UccZdc, processZdcCollAss, "Process ZDC W/Coll Ass.", true);
828852

853+
void processEventSampling(o2::aod::ColEvSels::iterator const& collision, o2::aod::BCsRun3 const& /*bcs*/, aod::Zdcs const& /*zdcs*/, aod::FV0As const& /*fv0as*/, aod::FT0s const& /*ft0s*/, TheFilteredTracks const& tracks)
854+
{
855+
if (!isEventSelected(collision)) {
856+
return;
857+
}
858+
859+
const auto& foundBC = collision.foundBC_as<o2::aod::BCsRun3>();
860+
861+
// has ZDC?
862+
if (!foundBC.has_zdc()) {
863+
return;
864+
}
865+
registry.fill(HIST("hEventCounter"), EvCutLabel::Zdc);
866+
867+
float aT0A = 0., aT0C = 0.;
868+
if (foundBC.has_ft0()) {
869+
for (const auto& amplitude : foundBC.ft0().amplitudeA()) {
870+
aT0A += amplitude;
871+
}
872+
for (const auto& amplitude : foundBC.ft0().amplitudeC()) {
873+
aT0C += amplitude;
874+
}
875+
} else {
876+
return;
877+
}
878+
879+
registry.fill(HIST("hEventCounter"), EvCutLabel::TZero);
880+
881+
const double normT0M{(aT0A + aT0C) / 100.};
882+
float znA{foundBC.zdc().amplitudeZNA()};
883+
float znC{foundBC.zdc().amplitudeZNC()};
884+
float aZEM1{foundBC.zdc().amplitudeZEM1()};
885+
float aZEM2{foundBC.zdc().amplitudeZEM2()};
886+
float tZNA{foundBC.zdc().timeZNA()};
887+
float tZNC{foundBC.zdc().timeZNC()};
888+
float tZPA{foundBC.zdc().timeZPA()};
889+
float tZPC{foundBC.zdc().timeZPC()};
890+
float tZDCdif{tZNC + tZPC - tZNA - tZPA};
891+
float tZDCsum{tZNC + tZPC + tZNA + tZPA};
892+
znA /= kCollEnergy;
893+
znC /= kCollEnergy;
894+
float sumZNs{znA + znC};
895+
float sumZEMs{aZEM1 + aZEM2};
896+
897+
// TDC cut
898+
if (isTDCcut) {
899+
if (std::sqrt(std::pow(tZDCdif, 2.) + std::pow(tZDCsum, 2.)) > tdcCut) {
900+
return;
901+
}
902+
registry.fill(HIST("hEventCounter"), EvCutLabel::Tdc);
903+
}
904+
905+
// ZEM cut
906+
if (isZEMcut) {
907+
if (sumZEMs < zemCut) {
908+
return;
909+
}
910+
registry.fill(HIST("hEventCounter"), EvCutLabel::Zem);
911+
}
912+
913+
registry.fill(HIST("zPos"), collision.posZ());
914+
registry.fill(HIST("T0Ccent"), collision.centFT0C());
915+
916+
// Nch-based selection
917+
int glbTracks = 0;
918+
for (const auto& track : tracks) {
919+
// Track Selection
920+
// if (track.hasITS()) { itsTracks++; }
921+
if (!track.isGlobalTrack()) {
922+
continue;
923+
}
924+
if ((track.pt() < minPt) || (track.pt() > maxPt)) {
925+
continue;
926+
}
927+
registry.fill(HIST("ZposVsEta"), collision.posZ(), track.eta());
928+
registry.fill(HIST("EtaVsPhi"), track.eta(), track.phi());
929+
registry.fill(HIST("sigma1Pt"), track.pt(), track.sigma1Pt());
930+
registry.fill(HIST("dcaXYvspT"), track.dcaXY(), track.pt());
931+
glbTracks++;
932+
}
933+
934+
bool skipEvent{false};
935+
if (useMidRapNchSel) {
936+
auto hMeanNch = ccdb->getForTimeStamp<TH1F>(paTHmeanNch.value, foundBC.timestamp());
937+
auto hSigmaNch = ccdb->getForTimeStamp<TH1F>(paTHsigmaNch.value, foundBC.timestamp());
938+
if (!hMeanNch) {
939+
LOGF(info, "hMeanNch NOT LOADED!");
940+
return;
941+
}
942+
if (!hSigmaNch) {
943+
LOGF(info, "hSigmaNch NOT LOADED!");
944+
return;
945+
}
946+
947+
const int binT0M{hMeanNch->FindBin(normT0M)};
948+
const double meanNch{hMeanNch->GetBinContent(binT0M)};
949+
const double sigmaNch{hSigmaNch->GetBinContent(binT0M)};
950+
const double nSigmaSelection{nSigmaNchCut * sigmaNch};
951+
const double diffMeanNch{meanNch - glbTracks};
952+
953+
if (!(std::abs(diffMeanNch) < nSigmaSelection)) {
954+
registry.fill(HIST("ExcludedEvtVsFT0M"), normT0M);
955+
registry.fill(HIST("ExcludedEvtVsNch"), glbTracks);
956+
} else {
957+
skipEvent = true;
958+
}
959+
}
960+
961+
// Skip event based on number of Nch sigmas
962+
if (!skipEvent) {
963+
return;
964+
}
965+
966+
auto efficiency = ccdb->getForTimeStamp<TH2F>(paTHEff.value, foundBC.timestamp());
967+
if (!efficiency) {
968+
return;
969+
}
970+
971+
auto feedDown = ccdb->getForTimeStamp<TH2F>(paTHFD.value, foundBC.timestamp());
972+
if (!feedDown) {
973+
return;
974+
}
975+
976+
//---------------------------------------------------
977+
978+
uint64_t timeStamp{foundBC.timestamp()};
979+
980+
TRandom3 rndGen(timeStamp);
981+
std::vector<uint64_t> vPoisson;
982+
for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) {
983+
vPoisson.emplace_back(rndGen.Poisson(1.));
984+
}
985+
986+
for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) {
987+
988+
hPoisson[replica]->Fill(vPoisson.at(replica));
989+
990+
for (uint64_t evtRep = 0; evtRep < vPoisson.at(replica); ++evtRep) {
991+
992+
double nchMult{0.};
993+
std::vector<float> pTs;
994+
std::vector<float> vecFD;
995+
std::vector<float> vecEff;
996+
997+
// Calculates the Nch multiplicity
998+
for (const auto& track : tracks) {
999+
// Track Selection
1000+
if (!track.isGlobalTrack()) {
1001+
continue;
1002+
}
1003+
if ((track.pt() < minPt) || (track.pt() > maxPt)) {
1004+
continue;
1005+
}
1006+
1007+
float pt{track.pt()};
1008+
int foundNchBin{efficiency->GetXaxis()->FindBin(glbTracks)};
1009+
int foundPtBin{efficiency->GetYaxis()->FindBin(pt)};
1010+
float effValue{1.};
1011+
float fdValue{1.};
1012+
if (applyEff) {
1013+
effValue = efficiency->GetBinContent(foundNchBin, foundPtBin);
1014+
}
1015+
if (applyFD) {
1016+
fdValue = feedDown->GetBinContent(foundNchBin, foundPtBin);
1017+
}
1018+
if ((effValue > 0.) && (fdValue > 0.)) {
1019+
nchMult += (std::pow(effValue, -1.) * fdValue);
1020+
}
1021+
}
1022+
1023+
if (!applyEff)
1024+
nchMult = static_cast<double>(glbTracks);
1025+
if (applyEff && !correctNch)
1026+
nchMult = static_cast<double>(glbTracks);
1027+
if (nchMult < minNchSel) {
1028+
return;
1029+
}
1030+
1031+
// Fill vectors for [pT] measurement
1032+
pTs.clear();
1033+
vecFD.clear();
1034+
vecEff.clear();
1035+
for (const auto& track : tracks) {
1036+
// Track Selection
1037+
if (!track.isGlobalTrack()) {
1038+
continue;
1039+
}
1040+
if ((track.pt() < minPt) || (track.pt() > maxPtSpectra)) {
1041+
continue;
1042+
}
1043+
1044+
float pt{track.pt()};
1045+
int foundNchBin{efficiency->GetXaxis()->FindBin(glbTracks)};
1046+
int foundPtBin{efficiency->GetYaxis()->FindBin(pt)};
1047+
float effValue{1.};
1048+
float fdValue{1.};
1049+
if (applyEff) {
1050+
effValue = efficiency->GetBinContent(foundNchBin, foundPtBin);
1051+
fdValue = feedDown->GetBinContent(foundNchBin, foundPtBin);
1052+
}
1053+
if (applyEff && !applyFD) {
1054+
fdValue = 1.0;
1055+
}
1056+
if ((effValue > 0.) && (fdValue > 0.)) {
1057+
pTs.emplace_back(pt);
1058+
vecEff.emplace_back(effValue);
1059+
vecFD.emplace_back(fdValue);
1060+
}
1061+
}
1062+
1063+
double p1, p2, p3, p4, w1, w2, w3, w4;
1064+
p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0;
1065+
getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4);
1066+
1067+
// EbE one-particle pT correlation
1068+
double oneParCorr{p1 / w1};
1069+
1070+
// EbE two-particle pT correlation
1071+
double denTwoParCorr{std::pow(w1, 2.) - w2};
1072+
double numTwoParCorr{std::pow(p1, 2.) - p2};
1073+
double twoParCorr{numTwoParCorr / denTwoParCorr};
1074+
1075+
// EbE three-particle pT correlation
1076+
double denThreeParCorr{std::pow(w1, 3.) - 3. * w2 * w1 + 2. * w3};
1077+
double numThreeParCorr{std::pow(p1, 3.) - 3. * p2 * p1 + 2. * p3};
1078+
double threeParCorr{numThreeParCorr / denThreeParCorr};
1079+
1080+
// EbE four-particle pT correlation
1081+
double denFourParCorr{std::pow(w1, 4.) - 6. * w2 * std::pow(w1, 2.) + 3. * std::pow(w2, 2.) + 8 * w3 * w1 - 6. * w4};
1082+
double numFourParCorr{std::pow(p1, 4.) - 6. * p2 * std::pow(p1, 2.) + 3. * std::pow(p2, 2.) + 8 * p3 * p1 - 6. * p4};
1083+
double fourParCorr{numFourParCorr / denFourParCorr};
1084+
1085+
hNch[replica]->Fill(nchMult);
1086+
pNchVsOneParCorr[replica]->Fill(nchMult, oneParCorr, w1);
1087+
pNchVsOneParCorrVsZN[replica]->Fill(nchMult, sumZNs, oneParCorr, w1);
1088+
pNchVsTwoParCorrVsZN[replica]->Fill(nchMult, sumZNs, twoParCorr, denTwoParCorr);
1089+
pNchVsThreeParCorrVsZN[replica]->Fill(nchMult, sumZNs, threeParCorr, denThreeParCorr);
1090+
pNchVsFourParCorrVsZN[replica]->Fill(nchMult, sumZNs, fourParCorr, denFourParCorr);
1091+
}
1092+
}
1093+
}
1094+
PROCESS_SWITCH(UccZdc, processEventSampling, "Process Event Sampling 4 Bootstrap", true);
1095+
8291096
// Preslice<aod::McParticles> perMCCollision = aod::mcparticle::mcCollisionId;
8301097
Preslice<TheFilteredSimTracks> perCollision = aod::track::collisionId;
8311098
Service<o2::framework::O2DatabasePDG> pdg;

0 commit comments

Comments
 (0)