Skip to content

Commit f7471ea

Browse files
omassenomassen
committed
Implement NUA-weights for flow analysis
1 parent 7e64819 commit f7471ea

1 file changed

Lines changed: 44 additions & 0 deletions

File tree

‎PWGEM/Dilepton/Core/Dilepton.h‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -149,6 +149,7 @@ struct Dilepton {
149149
o2::framework::ConfigurableAxis ConfPolarizationQuadMomBins{"ConfPolarizationQuadMomBins", {15, -0.5, 1}, "quadrupole moment bins for polarization analysis"}; // quardrupole moment <(3 x cos^2(theta) -1)/2>
150150

151151
o2::framework::Configurable<int> cfgNumBootstrapSamples{"cfgNumBootstrapSamples", 1, "Number of Bootstrap Samples"};
152+
152153

153154
EMEventCut fEMEventCut;
154155
struct : o2::framework::ConfigurableGroup {
@@ -323,6 +324,15 @@ struct Dilepton {
323324
o2::framework::Configurable<uint64_t> bcMarginForSoftwareTrigger{"bcMarginForSoftwareTrigger", 100, "Number of BCs of margin for software triggers"};
324325
} zorroGroup;
325326

327+
struct : o2::framework::ConfigurableGroup {
328+
std::string prefix = "flowcorrection_group";
329+
o2::framework::Configurable<bool> cfgApplyWeightNUA{"cfgApplyWeightNUA", false, "flag to apply q-vector Non-uniform acceptance weighting"};
330+
o2::framework::Configurable<std::string> nuaPath{"nuaPath", "Users/o/omassen/Dielectron/NUAWeights/LHC23_PbPb_pass5", "Path to NUA-weights file"};
331+
o2::framework::ConfigurableAxis ConfNUAPhiBins{"ConfNUAPhiBins", {60, -M_PI, M_PI}, "NUA histogram bins - phi-angle"};
332+
o2::framework::ConfigurableAxis ConfNUAEtaBins{"ConfNUAEtaBins", {40, -1., +1.}, "NUA histogram bins - eta"};
333+
o2::framework::ConfigurableAxis ConfNUAZVtxBins{"ConfNUAZVtxBins", {40, -10., 10.}, "NUA histogram bins - z-vtx"};
334+
} flowcorrectionGroup;
335+
326336
Zorro zorro;
327337
int mToIidx = 0;
328338
int mTOICounter = 0;
@@ -355,6 +365,7 @@ struct Dilepton {
355365
float beamP1 = 0.f; // beam momentum
356366
float beamP2 = 0.f; // beam momentum
357367
TH2D* h2sp_resolution = nullptr;
368+
std::vector<TH3D*> h3nua_weights;
358369

359370
void init(o2::framework::InitContext& /*context*/)
360371
{
@@ -558,6 +569,16 @@ struct Dilepton {
558569
h2sp_resolution = reinterpret_cast<TH2D*>(list->FindObject(spresoHistName.value.data()));
559570
LOGF(info, "h2sp_resolution.GetBinContent(40, 1) = %f", h2sp_resolution->GetBinContent(40, 1));
560571
}
572+
573+
if(flowcorrectionGroup.cfgApplyWeightNUA){
574+
auto list_nua = ccdb->getForTimeStamp<TList>(flowcorrectionGroup.nuaPath, collision.timestamp());
575+
h3nua_weights.push_back(reinterpret_cast<TH3D*>(list_nua->FindObject("weights_uls_NUA")));
576+
LOGF(info, "h3nua_weights.at(0)->GetBinContent(5, 5, 5) = %f", h3nua_weights.at(0)->GetBinContent(5, 5, 5));
577+
h3nua_weights.push_back(reinterpret_cast<TH3D*>(list_nua->FindObject("weights_lspp_NUA")));
578+
LOGF(info, "h3nua_weights.at(1)->GetBinContent(5, 5, 5) = %f", h3nua_weights.at(1)->GetBinContent(5, 5, 5));
579+
h3nua_weights.push_back(reinterpret_cast<TH3D*>(list_nua->FindObject("weights_lsmm_NUA")));
580+
LOGF(info, "h3nua_weights.at(2)->GetBinContent(5, 5, 5) = %f", h3nua_weights.at(2)->GetBinContent(5, 5, 5));
581+
}
561582
}
562583

563584
~Dilepton()
@@ -648,6 +669,7 @@ struct Dilepton {
648669
const o2::framework::AxisSpec axis_sp{ConfSPBins, Form("#vec{u}_{%d,ll} #upoint #vec{Q}_{%d}^{%s}", nmod, nmod, qvec_det_names[cfgQvecEstimator].data())};
649670

650671
fRegistry.add("Pair/same/uls/hs", "dilepton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_sp}, true);
672+
fRegistry.add("Pair/same/uls/hNUA", "NUA Histogram;#phi (rad.);#eta;VtxZ;", o2::framework::HistType::kTH3D, {flowcorrectionGroup.ConfNUAPhiBins, flowcorrectionGroup.ConfNUAEtaBins, flowcorrectionGroup.ConfNUAZVtxBins}, true);
651673
fRegistry.addClone("Pair/same/uls/", "Pair/same/lspp/");
652674
fRegistry.addClone("Pair/same/uls/", "Pair/same/lsmm/");
653675

@@ -880,6 +902,22 @@ struct Dilepton {
880902
}
881903
}
882904

905+
float getNUAweight(const int type_int, const float phi, const float eta, const float zVtx)
906+
{
907+
if (h3nua_weights.at(type_int) == nullptr){
908+
return 1.f;
909+
}
910+
int binId_phi = h3nua_weights.at(type_int)->GetXaxis()->FindBin(phi);
911+
int binId_eta = h3nua_weights.at(type_int)->GetYaxis()->FindBin(eta);
912+
int binId_zVtx = h3nua_weights.at(type_int)->GetZaxis()->FindBin(zVtx);
913+
float nuaWeight = h3nua_weights.at(type_int)->GetBinContent(binId_phi, binId_eta, binId_zVtx);
914+
915+
if (nuaWeight == 0 || std::isnan(nuaWeight) || std::isinf(nuaWeight)) {
916+
nuaWeight = 1.f;
917+
}
918+
return nuaWeight;
919+
}
920+
883921
template <int ev_id, typename TCollision, typename TTrack1, typename TTrack2, typename TCut, typename TAllTracks>
884922
bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const&, const std::vector<float> weightvector)
885923
{
@@ -1054,14 +1092,20 @@ struct Dilepton {
10541092

10551093
if constexpr (ev_id == 0) {
10561094
// LOGF(info, "collision.centFT0C() = %f, collision.trackOccupancyInTimeRange() = %d, getSPresolution = %f", collision.centFT0C(), collision.trackOccupancyInTimeRange(), getSPresolution(collision.centFT0C(), collision.trackOccupancyInTimeRange()));
1095+
if(flowcorrectionGroup.cfgApplyWeightNUA){
1096+
weight *= 1./getNUAweight(t1.sign() * t2.sign() < 0 ? 0 : (t1.sign() > 0 && t2.sign() > 0 ? 1 : 2), v12.Phi(), v12.Eta(), collision.posZ());
1097+
}
10571098

10581099
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());
10591100
if (t1.sign() * t2.sign() < 0) { // ULS
10601101
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, weight);
1102+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hNUA"), v12.Phi(), v12.Eta(), collision.posZ(), weight);
10611103
} else if (t1.sign() > 0 && t2.sign() > 0) { // LS++
10621104
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, weight);
1105+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lspp/hNUA"), v12.Phi(), v12.Eta(), collision.posZ(), weight);
10631106
} else if (t1.sign() < 0 && t2.sign() < 0) { // LS--
10641107
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), sp, weight);
1108+
fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hNUA"), v12.Phi(), v12.Eta(), collision.posZ(), weight);
10651109
}
10661110
} else if constexpr (ev_id == 1) {
10671111
if (t1.sign() * t2.sign() < 0) { // ULS

0 commit comments

Comments
 (0)