Skip to content

Commit a390108

Browse files
omassenalibuildvkucera
authored
[PWGEM] New Dielectron analysis type: Bootstrapping for v2. (#15576)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch> Co-authored-by: Vít Kučera <26327373+vkucera@users.noreply.github.com>
1 parent cecb2fe commit a390108

File tree

3 files changed

+130
-11
lines changed

3 files changed

+130
-11
lines changed

‎PWGEM/Dilepton/Core/Dilepton.h‎

Lines changed: 82 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,7 @@
6464
#include <TH1.h>
6565
#include <TH2.h>
6666
#include <TList.h>
67+
#include <TRandom3.h>
6768
#include <TString.h>
6869

6970
#include <sys/types.h>
@@ -113,7 +114,7 @@ struct Dilepton {
113114
o2::framework::Configurable<std::string> spresoPath{"spresoPath", "Users/d/dsekihat/PWGEM/dilepton/Qvector/resolution/LHC23zzh/pass3/test", "Path to SP resolution file"};
114115
o2::framework::Configurable<std::string> spresoHistName{"spresoHistName", "h1_R2_FT0M_BPos_BNeg", "histogram name of SP resolution file"};
115116

116-
o2::framework::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"};
117+
o2::framework::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"};
117118
o2::framework::Configurable<int> cfgEP2Estimator_for_Mix{"cfgEP2Estimator_for_Mix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5, FV0A:6"};
118119
o2::framework::Configurable<int> cfgQvecEstimator{"cfgQvecEstimator", 2, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5, FV0A:6"};
119120
o2::framework::Configurable<int> cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"};
@@ -142,6 +143,8 @@ struct Dilepton {
142143
o2::framework::ConfigurableAxis ConfPolarizationPhiBins{"ConfPolarizationPhiBins", {1, -M_PI, M_PI}, "phi bins for polarization analysis"};
143144
o2::framework::ConfigurableAxis ConfPolarizationQuadMomBins{"ConfPolarizationQuadMomBins", {15, -0.5, 1}, "quadrupole moment bins for polarization analysis"}; // quardrupole moment <(3 x cos^2(theta) -1)/2>
144145

146+
o2::framework::Configurable<int> cfgNumBootstrapSamples{"cfgNumBootstrapSamples", 1, "Number of Bootstrap Samples"};
147+
145148
EMEventCut fEMEventCut;
146149
struct : o2::framework::ConfigurableGroup {
147150
std::string prefix = "eventcut_group";
@@ -673,6 +676,17 @@ struct Dilepton {
673676
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
674677
fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/");
675678
fRegistry.addClone("Pair/same/", "Pair/mix/");
679+
} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) {
680+
nmod = 2;
681+
const o2::framework::AxisSpec axis_sp{ConfSPBins, Form("#vec{u}_{%d,ll} #upoint #vec{Q}_{%d}^{%s}", nmod, nmod, qvec_det_names[cfgQvecEstimator].data())};
682+
const o2::framework::AxisSpec axis_bootstrap{cfgNumBootstrapSamples, 0.5, static_cast<double>(cfgNumBootstrapSamples) + 0.5, "sample"}; // for bootstrap samples
683+
fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_sp, axis_bootstrap}, true);
684+
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
685+
fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/");
686+
fRegistry.add("Pair/mix/uls/hs", "dilepton", o2::framework::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true);
687+
fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lspp/");
688+
fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lsmm/");
689+
o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistogramsBootstrap(&fRegistry, cfgNumBootstrapSamples);
676690
} else { // same as kQC to avoid seg. fault
677691
fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true);
678692
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
@@ -844,7 +858,7 @@ struct Dilepton {
844858
}
845859

846860
template <int ev_id, typename TCollision, typename TTrack1, typename TTrack2, typename TCut, typename TAllTracks>
847-
bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&)
861+
bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&, const std::vector<float> weightvector)
848862
{
849863
if constexpr (ev_id == 0) {
850864
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
@@ -1057,7 +1071,6 @@ struct Dilepton {
10571071
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
10581072
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);
10591073
}
1060-
10611074
} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kHFll)) {
10621075
float dphi = v1.Phi() - v2.Phi();
10631076
dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf);
@@ -1071,6 +1084,51 @@ struct Dilepton {
10711084
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), dphi, deta, weight);
10721085
}
10731086

1087+
} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) {
1088+
std::array<float, 2> q2ft0m = {collision.q2xft0m(), collision.q2yft0m()};
1089+
std::array<float, 2> q2ft0a = {collision.q2xft0a(), collision.q2yft0a()};
1090+
std::array<float, 2> q2ft0c = {collision.q2xft0c(), collision.q2yft0c()};
1091+
std::array<float, 2> q2btot = {collision.q2xbtot(), collision.q2ybtot()};
1092+
std::array<float, 2> q2bpos = {collision.q2xbpos(), collision.q2ybpos()};
1093+
std::array<float, 2> q2bneg = {collision.q2xbneg(), collision.q2ybneg()};
1094+
std::array<float, 2> q2fv0a = {collision.q2xfv0a(), collision.q2yfv0a()};
1095+
std::array<float, 2> q3ft0m = {collision.q3xft0m(), collision.q3yft0m()};
1096+
std::array<float, 2> q3ft0a = {collision.q3xft0a(), collision.q3yft0a()};
1097+
std::array<float, 2> q3ft0c = {collision.q3xft0c(), collision.q3yft0c()};
1098+
std::array<float, 2> q3btot = {collision.q3xbtot(), collision.q3ybtot()};
1099+
std::array<float, 2> q3bpos = {collision.q3xbpos(), collision.q3ybpos()};
1100+
std::array<float, 2> q3bneg = {collision.q3xbneg(), collision.q3ybneg()};
1101+
std::array<float, 2> q3fv0a = {collision.q3xfv0a(), collision.q3yfv0a()};
1102+
1103+
std::vector<std::vector<std::array<float, 2>>> qvectors = {
1104+
{{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
1105+
{{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
1106+
{q2ft0m, q2ft0a, q2ft0c, q2btot, q2bpos, q2bneg, q2fv0a}, // 2nd harmonics
1107+
{q3ft0m, q3ft0a, q3ft0c, q3btot, q3bpos, q3bneg, q3fv0a}, // 3rd harmonics
1108+
};
1109+
1110+
if constexpr (ev_id == 0) {
1111+
// LOGF(info, "collision.centFT0C() = %f, collision.trackOccupancyInTimeRange() = %d, getSPresolution = %f", collision.centFT0C(), collision.trackOccupancyInTimeRange(), getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange()));
1112+
1113+
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());
1114+
for (int i = 0; i < cfgNumBootstrapSamples; i++) {
1115+
if (t1.sign() * t2.sign() < 0) { // ULS
1116+
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));
1117+
} else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
1118+
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));
1119+
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
1120+
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));
1121+
}
1122+
}
1123+
} else if constexpr (ev_id == 1) {
1124+
if (t1.sign() * t2.sign() < 0) { // ULS
1125+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
1126+
} else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
1127+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
1128+
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
1129+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
1130+
}
1131+
}
10741132
} else { // same as kQC to avoid seg. fault
10751133
if (t1.sign() * t2.sign() < 0) { // ULS
10761134
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
@@ -1266,6 +1324,20 @@ struct Dilepton {
12661324
}
12671325
}
12681326

1327+
std::vector<float> bootstrapweights = {};
1328+
if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) { // bootstrapping for accepted events
1329+
int randomSeed = static_cast<int>(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
1330+
TRandom3 randomNumber(randomSeed);
1331+
for (int i = 0; i < cfgNumBootstrapSamples; i++) {
1332+
float poissonweight = 0.;
1333+
poissonweight = static_cast<float>(randomNumber.PoissonD(1.0));
1334+
bootstrapweights.push_back(poissonweight);
1335+
o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfoBootstrap(&fRegistry, collision, i, poissonweight);
1336+
}
1337+
} else {
1338+
bootstrapweights.push_back(1.0); // to pass as non-empyt dummy to use
1339+
}
1340+
12691341
if (nmod == 2) {
12701342
o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, 2>(&fRegistry, collision);
12711343
} else if (nmod == 3) {
@@ -1286,19 +1358,19 @@ struct Dilepton {
12861358
used_trackIds_per_col.reserve(posTracks_per_coll.size() + negTracks_per_coll.size());
12871359
int nuls = 0, nlspp = 0, nlsmm = 0;
12881360
for (const auto& [pos, neg] : combinations(o2::soa::CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS
1289-
bool is_pair_ok = fillPairInfo<0>(collision, pos, neg, cut, tracks);
1361+
bool is_pair_ok = fillPairInfo<0>(collision, pos, neg, cut, tracks, bootstrapweights);
12901362
if (is_pair_ok) {
12911363
nuls++;
12921364
}
12931365
}
12941366
for (const auto& [pos1, pos2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++
1295-
bool is_pair_ok = fillPairInfo<0>(collision, pos1, pos2, cut, tracks);
1367+
bool is_pair_ok = fillPairInfo<0>(collision, pos1, pos2, cut, tracks, bootstrapweights);
12961368
if (is_pair_ok) {
12971369
nlspp++;
12981370
}
12991371
}
13001372
for (const auto& [neg1, neg2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS--
1301-
bool is_pair_ok = fillPairInfo<0>(collision, neg1, neg2, cut, tracks);
1373+
bool is_pair_ok = fillPairInfo<0>(collision, neg1, neg2, cut, tracks, bootstrapweights);
13021374
if (is_pair_ok) {
13031375
nlsmm++;
13041376
}
@@ -1380,25 +1452,25 @@ struct Dilepton {
13801452

13811453
for (const auto& pos : selected_posTracks_in_this_event) { // ULS mix
13821454
for (const auto& neg : negTracks_from_event_pool) {
1383-
fillPairInfo<1>(collision, pos, neg, cut, nullptr);
1455+
fillPairInfo<1>(collision, pos, neg, cut, nullptr, bootstrapweights);
13841456
}
13851457
}
13861458

13871459
for (const auto& neg : selected_negTracks_in_this_event) { // ULS mix
13881460
for (const auto& pos : posTracks_from_event_pool) {
1389-
fillPairInfo<1>(collision, neg, pos, cut, nullptr);
1461+
fillPairInfo<1>(collision, neg, pos, cut, nullptr, bootstrapweights);
13901462
}
13911463
}
13921464

13931465
for (const auto& pos1 : selected_posTracks_in_this_event) { // LS++ mix
13941466
for (const auto& pos2 : posTracks_from_event_pool) {
1395-
fillPairInfo<1>(collision, pos1, pos2, cut, nullptr);
1467+
fillPairInfo<1>(collision, pos1, pos2, cut, nullptr, bootstrapweights);
13961468
}
13971469
}
13981470

13991471
for (const auto& neg1 : selected_negTracks_in_this_event) { // LS-- mix
14001472
for (const auto& neg2 : negTracks_from_event_pool) {
1401-
fillPairInfo<1>(collision, neg1, neg2, cut, nullptr);
1473+
fillPairInfo<1>(collision, neg1, neg2, cut, nullptr, bootstrapweights);
14021474
}
14031475
}
14041476
} // end of loop over mixed event pool

0 commit comments

Comments
 (0)