Skip to content

Commit 1d75a04

Browse files
omassenomassen
committed
Introducing Bootstrapping to Dilepton analysis framework
1 parent 192165d commit 1d75a04

File tree

1 file changed

+84
-10
lines changed

1 file changed

+84
-10
lines changed

‎PWGEM/Dilepton/Core/Dilepton.h‎

Lines changed: 84 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,7 @@
4747

4848
#include "Math/Vector4D.h"
4949
#include "TString.h"
50+
#include "TRandom3.h"
5051

5152
#include <algorithm>
5253
#include <array>
@@ -97,7 +98,7 @@ struct Dilepton {
9798
Configurable<std::string> spresoPath{"spresoPath", "Users/d/dsekihat/PWGEM/dilepton/Qvector/resolution/LHC23zzh/pass3/test", "Path to SP resolution file"};
9899
Configurable<std::string> spresoHistName{"spresoHistName", "h1_R2_FT0M_BPos_BNeg", "histogram name of SP resolution file"};
99100

100-
Configurable<int> cfgAnalysisType{"cfgAnalysisType", static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kQC), "kQC:0, kUPC:1, kFlowV2:2, kFlowV3:3, kPolarization:4, kHFll:5"};
101+
Configurable<int> cfgAnalysisType{"cfgAnalysisType", static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kQC), "kQC:0, kUPC:1, kFlowV2:2, kFlowV3:3, kPolarization:4, kHFll:5, kBootstrapv2:6"};
101102
Configurable<int> cfgEP2Estimator_for_Mix{"cfgEP2Estimator_for_Mix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5, FV0A:6"};
102103
Configurable<int> cfgQvecEstimator{"cfgQvecEstimator", 2, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5, FV0A:6"};
103104
Configurable<int> cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"};
@@ -126,6 +127,8 @@ struct Dilepton {
126127
ConfigurableAxis ConfPolarizationPhiBins{"ConfPolarizationPhiBins", {1, -M_PI, M_PI}, "phi bins for polarization analysis"};
127128
ConfigurableAxis ConfPolarizationQuadMomBins{"ConfPolarizationQuadMomBins", {15, -0.5, 1}, "quadrupole moment bins for polarization analysis"}; // quardrupole moment <(3 x cos^2(theta) -1)/2>
128129

130+
Configurable<int> cfgNumBootstrapSamples{"cfgNumBootstrapSamples", 1, "Number of Bootstrap Samples"};
131+
129132
EMEventCut fEMEventCut;
130133
struct : ConfigurableGroup {
131134
std::string prefix = "eventcut_group";
@@ -657,6 +660,17 @@ struct Dilepton {
657660
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
658661
fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/");
659662
fRegistry.addClone("Pair/same/", "Pair/mix/");
663+
} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) {
664+
nmod = 2;
665+
const AxisSpec axis_sp{ConfSPBins, Form("#vec{u}_{%d,ll} #upoint #vec{Q}_{%d}^{%s}", nmod, nmod, qvec_det_names[cfgQvecEstimator].data())};
666+
const AxisSpec axis_bootstrap{cfgNumBootstrapSamples, 0.5, static_cast<double>(cfgNumBootstrapSamples) + 0.5, "sample"}; // for bootstrap samples
667+
fRegistry.add("Pair/same/uls/hs", "dilepton", kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_sp, axis_bootstrap}, true);
668+
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
669+
fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/");
670+
fRegistry.add("Pair/mix/uls/hs", "dilepton", kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true);
671+
fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lspp/");
672+
fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lsmm/");
673+
o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistogramsBootstrap(&fRegistry, cfgNumBootstrapSamples);
660674
} else { // same as kQC to avoid seg. fault
661675
fRegistry.add("Pair/same/uls/hs", "dilepton", kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true);
662676
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
@@ -672,6 +686,7 @@ struct Dilepton {
672686
} else {
673687
o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<-1>(&fRegistry);
674688
}
689+
675690
fRegistry.add("Event/before/hEP2_CentFT0C_forMix", Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].data()), kTH2F, {{110, 0, 110}, {180, -M_PI_2, +M_PI_2}}, false);
676691
fRegistry.add("Event/after/hEP2_CentFT0C_forMix", Form("2nd harmonics event plane for mix;centrality FT0C (%%);#Psi_{2}^{%s} (rad.)", qvec_det_names[cfgEP2Estimator_for_Mix].data()), kTH2F, {{110, 0, 110}, {180, -M_PI_2, +M_PI_2}}, false);
677692
}
@@ -828,7 +843,7 @@ struct Dilepton {
828843
}
829844

830845
template <int ev_id, typename TCollision, typename TTrack1, typename TTrack2, typename TCut, typename TAllTracks>
831-
bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&)
846+
bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&, const std::vector<float> weightvector)
832847
{
833848
if constexpr (ev_id == 0) {
834849
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
@@ -1041,7 +1056,6 @@ struct Dilepton {
10411056
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
10421057
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), cos_thetaPol, phiPol, quadmom, weight);
10431058
}
1044-
10451059
} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kHFll)) {
10461060
float dphi = v1.Phi() - v2.Phi();
10471061
dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf);
@@ -1055,6 +1069,51 @@ struct Dilepton {
10551069
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), dphi, deta, weight);
10561070
}
10571071

1072+
} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) {
1073+
std::array<float, 2> q2ft0m = {collision.q2xft0m(), collision.q2yft0m()};
1074+
std::array<float, 2> q2ft0a = {collision.q2xft0a(), collision.q2yft0a()};
1075+
std::array<float, 2> q2ft0c = {collision.q2xft0c(), collision.q2yft0c()};
1076+
std::array<float, 2> q2btot = {collision.q2xbtot(), collision.q2ybtot()};
1077+
std::array<float, 2> q2bpos = {collision.q2xbpos(), collision.q2ybpos()};
1078+
std::array<float, 2> q2bneg = {collision.q2xbneg(), collision.q2ybneg()};
1079+
std::array<float, 2> q2fv0a = {collision.q2xfv0a(), collision.q2yfv0a()};
1080+
std::array<float, 2> q3ft0m = {collision.q3xft0m(), collision.q3yft0m()};
1081+
std::array<float, 2> q3ft0a = {collision.q3xft0a(), collision.q3yft0a()};
1082+
std::array<float, 2> q3ft0c = {collision.q3xft0c(), collision.q3yft0c()};
1083+
std::array<float, 2> q3btot = {collision.q3xbtot(), collision.q3ybtot()};
1084+
std::array<float, 2> q3bpos = {collision.q3xbpos(), collision.q3ybpos()};
1085+
std::array<float, 2> q3bneg = {collision.q3xbneg(), collision.q3ybneg()};
1086+
std::array<float, 2> q3fv0a = {collision.q3xfv0a(), collision.q3yfv0a()};
1087+
1088+
std::vector<std::vector<std::array<float, 2>>> qvectors = {
1089+
{{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 0th harmonics
1090+
{{999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}, {999.f, 999.f}}, // 1st harmonics
1091+
{q2ft0m, q2ft0a, q2ft0c, q2btot, q2bpos, q2bneg, q2fv0a}, // 2nd harmonics
1092+
{q3ft0m, q3ft0a, q3ft0c, q3btot, q3bpos, q3bneg, q3fv0a}, // 3rd harmonics
1093+
};
1094+
1095+
if constexpr (ev_id == 0) {
1096+
// LOGF(info, "collision.centFT0C() = %f, collision.trackOccupancyInTimeRange() = %d, getSPresolution = %f", collision.centFT0C(), collision.trackOccupancyInTimeRange(), getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange()));
1097+
1098+
float sp = RecoDecay::dotProd(std::array<float, 2>{static_cast<float>(std::cos(nmod * v12.Phi())), static_cast<float>(std::sin(nmod * v12.Phi()))}, qvectors[nmod][cfgQvecEstimator]) / getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange());
1099+
for(int i = 0; i<cfgNumBootstrapSamples; i++){
1100+
if (t1.sign() * t2.sign() < 0) { // ULS
1101+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, i + 0.5, weightvector.at(i));
1102+
} else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
1103+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, i + 0.5, weightvector.at(i));
1104+
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
1105+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, i + 0.5, weightvector.at(i));
1106+
}
1107+
}
1108+
} else if constexpr (ev_id == 1) {
1109+
if (t1.sign() * t2.sign() < 0) { // ULS
1110+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
1111+
} else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
1112+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
1113+
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
1114+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
1115+
}
1116+
}
10581117
} else { // same as kQC to avoid seg. fault
10591118
if (t1.sign() * t2.sign() < 0) { // ULS
10601119
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
@@ -1250,6 +1309,21 @@ struct Dilepton {
12501309
}
12511310
}
12521311

1312+
std::vector<float> bootstrapweights = {};
1313+
int randomSeed = static_cast<int>(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
1314+
TRandom3 randomNumber(randomSeed);
1315+
if(cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)){ // bootstrapping for accepted events
1316+
for(int i = 0; i<cfgNumBootstrapSamples; i++){
1317+
float poissonweight = 0.;
1318+
poissonweight = static_cast<float>(randomNumber.PoissonD(1.0));
1319+
bootstrapweights.push_back(poissonweight);
1320+
o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfoBootstrap(&fRegistry, collision, i, poissonweight);
1321+
}
1322+
}
1323+
else{
1324+
bootstrapweights.push_back(1.0); // to pass as non-empyt dummy to use
1325+
}
1326+
12531327
if (nmod == 2) {
12541328
o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, 2>(&fRegistry, collision);
12551329
} else if (nmod == 3) {
@@ -1270,19 +1344,19 @@ struct Dilepton {
12701344
used_trackIds_per_col.reserve(posTracks_per_coll.size() + negTracks_per_coll.size());
12711345
int nuls = 0, nlspp = 0, nlsmm = 0;
12721346
for (const auto& [pos, neg] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS
1273-
bool is_pair_ok = fillPairInfo<0>(collision, pos, neg, cut, tracks);
1347+
bool is_pair_ok = fillPairInfo<0>(collision, pos, neg, cut, tracks, bootstrapweights);
12741348
if (is_pair_ok) {
12751349
nuls++;
12761350
}
12771351
}
12781352
for (const auto& [pos1, pos2] : combinations(CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++
1279-
bool is_pair_ok = fillPairInfo<0>(collision, pos1, pos2, cut, tracks);
1353+
bool is_pair_ok = fillPairInfo<0>(collision, pos1, pos2, cut, tracks, bootstrapweights);
12801354
if (is_pair_ok) {
12811355
nlspp++;
12821356
}
12831357
}
12841358
for (const auto& [neg1, neg2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS--
1285-
bool is_pair_ok = fillPairInfo<0>(collision, neg1, neg2, cut, tracks);
1359+
bool is_pair_ok = fillPairInfo<0>(collision, neg1, neg2, cut, tracks, bootstrapweights);
12861360
if (is_pair_ok) {
12871361
nlsmm++;
12881362
}
@@ -1364,25 +1438,25 @@ struct Dilepton {
13641438

13651439
for (const auto& pos : selected_posTracks_in_this_event) { // ULS mix
13661440
for (const auto& neg : negTracks_from_event_pool) {
1367-
fillPairInfo<1>(collision, pos, neg, cut, nullptr);
1441+
fillPairInfo<1>(collision, pos, neg, cut, nullptr, bootstrapweights);
13681442
}
13691443
}
13701444

13711445
for (const auto& neg : selected_negTracks_in_this_event) { // ULS mix
13721446
for (const auto& pos : posTracks_from_event_pool) {
1373-
fillPairInfo<1>(collision, neg, pos, cut, nullptr);
1447+
fillPairInfo<1>(collision, neg, pos, cut, nullptr, bootstrapweights);
13741448
}
13751449
}
13761450

13771451
for (const auto& pos1 : selected_posTracks_in_this_event) { // LS++ mix
13781452
for (const auto& pos2 : posTracks_from_event_pool) {
1379-
fillPairInfo<1>(collision, pos1, pos2, cut, nullptr);
1453+
fillPairInfo<1>(collision, pos1, pos2, cut, nullptr, bootstrapweights);
13801454
}
13811455
}
13821456

13831457
for (const auto& neg1 : selected_negTracks_in_this_event) { // LS-- mix
13841458
for (const auto& neg2 : negTracks_from_event_pool) {
1385-
fillPairInfo<1>(collision, neg1, neg2, cut, nullptr);
1459+
fillPairInfo<1>(collision, neg1, neg2, cut, nullptr, bootstrapweights);
13861460
}
13871461
}
13881462
} // end of loop over mixed event pool

0 commit comments

Comments
 (0)