Skip to content

[PWGEM] New Dielectron analysis type: Bootstrapping for v2. #15576

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
omassen wants to merge 5 commits into
base: master
Choose a base branch
from
+130 −10
Draft

[PWGEM] New Dielectron analysis type: Bootstrapping for v2. #15576

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9f46bfe
changes to EventHistograms to inlude Bootstrapping Histograms
omassen Mar 30, 2026
192165d
Adding bootstrapping as new analysis type
omassen Mar 30, 2026
1d75a04
Introducing Bootstrapping to Dilepton analysis framework
omassen Mar 30, 2026
9741bb9
Please consider the following formatting changes
alibuild Mar 30, 2026
3a77b20
Merge pull request #10 from alibuild/alibot-cleanup-15576
omassen Mar 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
93 changes: 83 additions & 10 deletions PWGEM/Dilepton/Core/Dilepton.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
#include "MathUtils/Utils.h"

#include "Math/Vector4D.h"
#include "TRandom3.h"
#include "TString.h"

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

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"};
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"};
Configurable<int> cfgEP2Estimator_for_Mix{"cfgEP2Estimator_for_Mix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5, FV0A:6"};
Configurable<int> cfgQvecEstimator{"cfgQvecEstimator", 2, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5, FV0A:6"};
Configurable<int> cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"};
Expand Down Expand Up @@ -126,6 +127,8 @@ struct Dilepton {
ConfigurableAxis ConfPolarizationPhiBins{"ConfPolarizationPhiBins", {1, -M_PI, M_PI}, "phi bins for polarization analysis"};
ConfigurableAxis ConfPolarizationQuadMomBins{"ConfPolarizationQuadMomBins", {15, -0.5, 1}, "quadrupole moment bins for polarization analysis"}; // quardrupole moment <(3 x cos^2(theta) -1)/2>

Configurable<int> cfgNumBootstrapSamples{"cfgNumBootstrapSamples", 1, "Number of Bootstrap Samples"};

EMEventCut fEMEventCut;
struct : ConfigurableGroup {
std::string prefix = "eventcut_group";
Expand Down Expand Up @@ -657,6 +660,17 @@ struct Dilepton {
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/");
fRegistry.addClone("Pair/same/", "Pair/mix/");
} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) {
nmod = 2;
const AxisSpec axis_sp{ConfSPBins, Form("#vec{u}_{%d,ll} #upoint #vec{Q}_{%d}^{%s}", nmod, nmod, qvec_det_names[cfgQvecEstimator].data())};
const AxisSpec axis_bootstrap{cfgNumBootstrapSamples, 0.5, static_cast<double>(cfgNumBootstrapSamples) + 0.5, "sample"}; // for bootstrap samples
fRegistry.add("Pair/same/uls/hs", "dilepton", kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_sp, axis_bootstrap}, true);
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/");
fRegistry.add("Pair/mix/uls/hs", "dilepton", kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true);
fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lspp/");
fRegistry.addClone("Pair/mix/uls/", "Pair/mix/lsmm/");
o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistogramsBootstrap(&fRegistry, cfgNumBootstrapSamples);
} else { // same as kQC to avoid seg. fault
fRegistry.add("Pair/same/uls/hs", "dilepton", kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y}, true);
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
Expand All @@ -672,6 +686,7 @@ struct Dilepton {
} else {
o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<-1>(&fRegistry);
}

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);
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);
}
Expand Down Expand Up @@ -828,7 +843,7 @@ struct Dilepton {
}

template <int ev_id, typename TCollision, typename TTrack1, typename TTrack2, typename TCut, typename TAllTracks>
bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&)
bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&, const std::vector<float> weightvector)
{
if constexpr (ev_id == 0) {
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
Expand Down Expand Up @@ -1041,7 +1056,6 @@ struct Dilepton {
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
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);
}

} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kHFll)) {
float dphi = v1.Phi() - v2.Phi();
dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf);
Expand All @@ -1055,6 +1069,51 @@ struct Dilepton {
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), dphi, deta, weight);
}

} else if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) {
std::array<float, 2> q2ft0m = {collision.q2xft0m(), collision.q2yft0m()};
std::array<float, 2> q2ft0a = {collision.q2xft0a(), collision.q2yft0a()};
std::array<float, 2> q2ft0c = {collision.q2xft0c(), collision.q2yft0c()};
std::array<float, 2> q2btot = {collision.q2xbtot(), collision.q2ybtot()};
std::array<float, 2> q2bpos = {collision.q2xbpos(), collision.q2ybpos()};
std::array<float, 2> q2bneg = {collision.q2xbneg(), collision.q2ybneg()};
std::array<float, 2> q2fv0a = {collision.q2xfv0a(), collision.q2yfv0a()};
std::array<float, 2> q3ft0m = {collision.q3xft0m(), collision.q3yft0m()};
std::array<float, 2> q3ft0a = {collision.q3xft0a(), collision.q3yft0a()};
std::array<float, 2> q3ft0c = {collision.q3xft0c(), collision.q3yft0c()};
std::array<float, 2> q3btot = {collision.q3xbtot(), collision.q3ybtot()};
std::array<float, 2> q3bpos = {collision.q3xbpos(), collision.q3ybpos()};
std::array<float, 2> q3bneg = {collision.q3xbneg(), collision.q3ybneg()};
std::array<float, 2> q3fv0a = {collision.q3xfv0a(), collision.q3yfv0a()};

std::vector<std::vector<std::array<float, 2>>> qvectors = {
{{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
{{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
{q2ft0m, q2ft0a, q2ft0c, q2btot, q2bpos, q2bneg, q2fv0a}, // 2nd harmonics
{q3ft0m, q3ft0a, q3ft0c, q3btot, q3bpos, q3bneg, q3fv0a}, // 3rd harmonics
};

if constexpr (ev_id == 0) {
// LOGF(info, "collision.centFT0C() = %f, collision.trackOccupancyInTimeRange() = %d, getSPresolution = %f", collision.centFT0C(), collision.trackOccupancyInTimeRange(), getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange()));

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());
for (int i = 0; i < cfgNumBootstrapSamples; i++) {
if (t1.sign() * t2.sign() < 0) { // ULS
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));
} else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
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));
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
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));
}
}
} else if constexpr (ev_id == 1) {
if (t1.sign() * t2.sign() < 0) { // ULS
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
} else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
}
}
} else { // same as kQC to avoid seg. fault
if (t1.sign() * t2.sign() < 0) { // ULS
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), weight);
Expand Down Expand Up @@ -1250,6 +1309,20 @@ struct Dilepton {
}
}

std::vector<float> bootstrapweights = {};
int randomSeed = static_cast<int>(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
TRandom3 randomNumber(randomSeed);
if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kBootstrapv2)) { // bootstrapping for accepted events
for (int i = 0; i < cfgNumBootstrapSamples; i++) {
float poissonweight = 0.;
poissonweight = static_cast<float>(randomNumber.PoissonD(1.0));
bootstrapweights.push_back(poissonweight);
o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfoBootstrap(&fRegistry, collision, i, poissonweight);
}
} else {
bootstrapweights.push_back(1.0); // to pass as non-empyt dummy to use
}

if (nmod == 2) {
o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, 2>(&fRegistry, collision);
} else if (nmod == 3) {
Expand All @@ -1270,19 +1343,19 @@ struct Dilepton {
used_trackIds_per_col.reserve(posTracks_per_coll.size() + negTracks_per_coll.size());
int nuls = 0, nlspp = 0, nlsmm = 0;
for (const auto& [pos, neg] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS
bool is_pair_ok = fillPairInfo<0>(collision, pos, neg, cut, tracks);
bool is_pair_ok = fillPairInfo<0>(collision, pos, neg, cut, tracks, bootstrapweights);
if (is_pair_ok) {
nuls++;
}
}
for (const auto& [pos1, pos2] : combinations(CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++
bool is_pair_ok = fillPairInfo<0>(collision, pos1, pos2, cut, tracks);
bool is_pair_ok = fillPairInfo<0>(collision, pos1, pos2, cut, tracks, bootstrapweights);
if (is_pair_ok) {
nlspp++;
}
}
for (const auto& [neg1, neg2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS--
bool is_pair_ok = fillPairInfo<0>(collision, neg1, neg2, cut, tracks);
bool is_pair_ok = fillPairInfo<0>(collision, neg1, neg2, cut, tracks, bootstrapweights);
if (is_pair_ok) {
nlsmm++;
}
Expand Down Expand Up @@ -1364,25 +1437,25 @@ struct Dilepton {

for (const auto& pos : selected_posTracks_in_this_event) { // ULS mix
for (const auto& neg : negTracks_from_event_pool) {
fillPairInfo<1>(collision, pos, neg, cut, nullptr);
fillPairInfo<1>(collision, pos, neg, cut, nullptr, bootstrapweights);
}
}

for (const auto& neg : selected_negTracks_in_this_event) { // ULS mix
for (const auto& pos : posTracks_from_event_pool) {
fillPairInfo<1>(collision, neg, pos, cut, nullptr);
fillPairInfo<1>(collision, neg, pos, cut, nullptr, bootstrapweights);
}
}

for (const auto& pos1 : selected_posTracks_in_this_event) { // LS++ mix
for (const auto& pos2 : posTracks_from_event_pool) {
fillPairInfo<1>(collision, pos1, pos2, cut, nullptr);
fillPairInfo<1>(collision, pos1, pos2, cut, nullptr, bootstrapweights);
}
}

for (const auto& neg1 : selected_negTracks_in_this_event) { // LS-- mix
for (const auto& neg2 : negTracks_from_event_pool) {
fillPairInfo<1>(collision, neg1, neg2, cut, nullptr);
fillPairInfo<1>(collision, neg1, neg2, cut, nullptr, bootstrapweights);
}
}
} // end of loop over mixed event pool
Expand Down
Loading
Loading