Skip to content

[PWGEM/Dilepton] add flexibility of binning #13582

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

Merged
dsekihat merged 1 commit into from
Oct 29, 2025
Merged

[PWGEM/Dilepton] add flexibility of binning #13582

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
77 changes: 40 additions & 37 deletions PWGEM/Dilepton/Tasks/createResolutionMap.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,13 +82,15 @@
ConfigurableAxis ConfEtaCBGenBins{"ConfEtaCBGenBins", {30, -1.5, +1.5}, "gen. eta bins at midrapidity for output histograms"};
ConfigurableAxis ConfEtaFWDGenBins{"ConfEtaFWDGenBins", {40, -5.5, -1.5}, "gen. eta bins at forward rapidity for output histograms"};
ConfigurableAxis ConfPhiGenBins{"ConfPhiGenBins", {36, 0, 2.f * M_PI}, "gen. phi bins at forward rapidity for output histograms"};
ConfigurableAxis ConfPhiPositionCBGenBins{"ConfPhiPositionCBGenBins", {VARIABLE_WIDTH, 2.3 - M_PI, 0.85, 2.3, 0.85 + M_PI, 2.3 + M_PI}, "gen. phi psotion bins at forward rapidity for output histograms"}; // default is adjusted at R = 0.50 m
ConfigurableAxis ConfPhiPositionFWDGenBins{"ConfPhiPositionFWDGenBins", {1, 0, 2 * M_PI}, "gen. phi psotion bins at forward rapidity for output histograms"};
// ConfigurableAxis ConfPhiPositionCBGenBins{"ConfPhiPositionCBGenBins", {VARIABLE_WIDTH, 2.3 - M_PI, 0.85, 2.3, 0.85 + M_PI, 2.3 + M_PI}, "gen. phi psotion bins at forward rapidity for output histograms"}; // default is adjusted at R = 0.50 m
// ConfigurableAxis ConfPhiPositionFWDGenBins{"ConfPhiPositionFWDGenBins", {1, 0, 2 * M_PI}, "gen. phi psotion bins at forward rapidity for output histograms"};
Configurable<float> cfgRefR{"cfgRefR", 0.50, "ref. radius (m) for calculating phi position"}; // 0.50 +/- 0.06 can be syst. unc.

ConfigurableAxis ConfRelDeltaPtBins{"ConfRelDeltaPtBins", {200, -1.f, +1.f}, "rel. dpt for output histograms"};
ConfigurableAxis ConfDeltaEtaCBBins{"ConfDeltaEtaCBBins", {200, -0.5f, +0.5f}, "deta bins for output histograms"};
ConfigurableAxis ConfDeltaEtaFWDBins{"ConfDeltaEtaFWDBins", {200, -0.5f, +0.5f}, "deta bins for output histograms"};
ConfigurableAxis ConfRelDeltaPtCBBins{"ConfRelDeltaPtCBBins", {200, -1.f, +1.f}, "rel. dpt for output histograms at midrapidity"};
ConfigurableAxis ConfRelDeltaPtFWDBins{"ConfRelDeltaPtFWDBins", {200, -1.f, +1.f}, "rel. dpt for output histograms at fwd rapidity"};

ConfigurableAxis ConfDeltaEtaCBBins{"ConfDeltaEtaCBBins", {200, -0.5f, +0.5f}, "deta bins for output histograms at midrapidity"};
ConfigurableAxis ConfDeltaEtaFWDBins{"ConfDeltaEtaFWDBins", {200, -0.5f, +0.5f}, "deta bins for output histograms at fwd rapidity"};
ConfigurableAxis ConfDeltaPhiBins{"ConfDeltaPhiBins", {200, -0.5f, +0.5f}, "dphi bins for output histograms"};

Configurable<bool> cfgFillTHnSparse{"cfgFillTHnSparse", true, "fill THnSparse for output"};
Expand Down Expand Up @@ -193,12 +195,12 @@
o2::dataformats::VertexBase mVtx;
const o2::dataformats::MeanVertexObject* mMeanVtx = nullptr;
o2::base::MatLayerCylSet* lut = nullptr;
std::vector<float> phiPosition_bin_edges;
// std::vector<float> phiPosition_bin_edges;

~CreateResolutionMap()
{
phiPosition_bin_edges.clear();
phiPosition_bin_edges.shrink_to_fit();
// phiPosition_bin_edges.clear();
// phiPosition_bin_edges.shrink_to_fit();
}

void init(o2::framework::InitContext&)
Expand All @@ -224,35 +226,36 @@
d_bz = 0;
mBzMFT = 0;

if (ConfPhiPositionCBGenBins.value[0] == VARIABLE_WIDTH) {
phiPosition_bin_edges = std::vector<float>(ConfPhiPositionCBGenBins.value.begin(), ConfPhiPositionCBGenBins.value.end());
phiPosition_bin_edges.erase(phiPosition_bin_edges.begin());
// for (const auto& edge : phiPosition_bin_edges) {
// LOGF(info, "VARIABLE_WIDTH: phiPosition_bin_edges = %f", edge);
// }
} else { // FIXED bin width
int nbins = static_cast<int>(ConfPhiPositionCBGenBins.value[0]);
float xmin = static_cast<float>(ConfPhiPositionCBGenBins.value[1]);
float xmax = static_cast<float>(ConfPhiPositionCBGenBins.value[2]);
phiPosition_bin_edges.resize(nbins + 1);
for (int i = 0; i < nbins + 1; i++) {
phiPosition_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin;
// LOGF(info, "FIXED_WIDTH: phiPosition_bin_edges[%d] = %f", i, phiPosition_bin_edges[i]);
}
}
// if (ConfPhiPositionCBGenBins.value[0] == VARIABLE_WIDTH) {
// phiPosition_bin_edges = std::vector<float>(ConfPhiPositionCBGenBins.value.begin(), ConfPhiPositionCBGenBins.value.end());
// phiPosition_bin_edges.erase(phiPosition_bin_edges.begin());
// // for (const auto& edge : phiPosition_bin_edges) {
// // LOGF(info, "VARIABLE_WIDTH: phiPosition_bin_edges = %f", edge);
// // }
// } else { // FIXED bin width
// int nbins = static_cast<int>(ConfPhiPositionCBGenBins.value[0]);
// float xmin = static_cast<float>(ConfPhiPositionCBGenBins.value[1]);
// float xmax = static_cast<float>(ConfPhiPositionCBGenBins.value[2]);
// phiPosition_bin_edges.resize(nbins + 1);
// for (int i = 0; i < nbins + 1; i++) {
// phiPosition_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin;
// // LOGF(info, "FIXED_WIDTH: phiPosition_bin_edges[%d] = %f", i, phiPosition_bin_edges[i]);
// }
// }

const AxisSpec axis_cent{ConfCentBins, "centrality (%)"};
const AxisSpec axis_pt_gen{ConfPtGenBins, "p_{T,l}^{gen} (GeV/c)"};
const AxisSpec axis_eta_cb_gen{ConfEtaCBGenBins, "#eta_{l}^{gen}"};
const AxisSpec axis_eta_fwd_gen{ConfEtaFWDGenBins, "#eta_{l}^{gen}"};
const AxisSpec axis_phi_gen{ConfPhiGenBins, "#varphi_{l}^{gen} (rad.)"};
const AxisSpec axis_dpt{ConfRelDeltaPtBins, "(p_{T,l}^{gen} - p_{T,l}^{rec})/p_{T,l}^{gen}"};
const AxisSpec axis_dpt_cb{ConfRelDeltaPtCBBins, "(p_{T,l}^{gen} - p_{T,l}^{rec})/p_{T,l}^{gen}"};
const AxisSpec axis_dpt_fwd{ConfRelDeltaPtFWDBins, "(p_{T,l}^{gen} - p_{T,l}^{rec})/p_{T,l}^{gen}"};
const AxisSpec axis_deta_cb{ConfDeltaEtaCBBins, "#eta_{l}^{gen} - #eta_{l}^{rec}"};
const AxisSpec axis_deta_fwd{ConfDeltaEtaFWDBins, "#eta_{l}^{gen} - #eta_{l}^{rec}"};
const AxisSpec axis_dphi{ConfDeltaPhiBins, "#varphi_{l}^{gen} - #varphi_{l}^{rec} (rad.)"};
const AxisSpec axis_charge_gen{3, -1.5, +1.5, "true sign"};
const AxisSpec axis_phiPositionCB_gen{ConfPhiPositionCBGenBins, Form("#varphi^{*, gen} (rad.) at r_{xy} = %3.2f m", cfgRefR.value)};
const AxisSpec axis_phiPositionFWD_gen{ConfPhiPositionFWDGenBins, "#varphi^{*, gen} (rad.)"};
// const AxisSpec axis_phiPositionCB_gen{ConfPhiPositionCBGenBins, Form("#varphi^{*, gen} (rad.) at r_{xy} = %3.2f m", cfgRefR.value)};
// const AxisSpec axis_phiPositionFWD_gen{ConfPhiPositionFWDGenBins, "#varphi^{*, gen} (rad.)"};

// registry.add("Event/Electron/hImpPar_Centrality", "true imapact parameter vs. estimated centrality;impact parameter (fm);centrality (%)", kTH2F, {{200, 0, 20}, {110, 0, 110}}, true);
// registry.add("Event/Electron/hImpPar_Centrality", "true imapact parameter vs. estimated centrality;impact parameter (fm);centrality (%)", kTH2F, {{200, 0, 20}, {110, 0, 110}}, true);
Expand All @@ -262,24 +265,24 @@
if (cfgFillTH2) {
registry.add("Electron/hPt", "rec. p_{T,e};p_{T,e} (GeV/c)", kTH1F, {{1000, 0, 10}}, false);
registry.add("Electron/hEtaPhi", "rec. #eta vs. #varphi;#varphi_{e} (rad.);#eta_{e}", kTH2F, {{90, 0, 2 * M_PI}, {100, -5, +5}}, false);
registry.add("Electron/Ptgen_RelDeltaPt", "resolution", kTH2F, {{axis_pt_gen}, {axis_dpt}}, true);
registry.add("Electron/Ptgen_RelDeltaPt", "resolution", kTH2F, {{axis_pt_gen}, {axis_dpt_cb}}, true);
registry.add("Electron/Ptgen_DeltaEta", "resolution", kTH2F, {{axis_pt_gen}, {axis_deta_cb}}, true);
registry.add("Electron/Ptgen_DeltaPhi_Pos", "resolution", kTH2F, {{axis_pt_gen}, {axis_dphi}}, true);
registry.add("Electron/Ptgen_DeltaPhi_Neg", "resolution", kTH2F, {{axis_pt_gen}, {axis_dphi}}, true);

registry.add("StandaloneMuon/hPt", "rec. p_{T,#mu};p_{T,#mu} (GeV/c)", kTH1F, {{1000, 0, 10}}, false);
registry.add("StandaloneMuon/hEtaPhi", "rec. #eta vs. #varphi;#varphi_{#mu} (rad.);#eta_{#mu}", kTH2F, {{90, 0, 2 * M_PI}, {100, -5, +5}}, false);
registry.add("StandaloneMuon/Ptgen_RelDeltaPt", "resolution", kTH2F, {{axis_pt_gen}, {axis_dpt}}, true);
registry.add("StandaloneMuon/Ptgen_RelDeltaPt", "resolution", kTH2F, {{axis_pt_gen}, {axis_dpt_fwd}}, true);
registry.add("StandaloneMuon/Ptgen_DeltaEta", "resolution", kTH2F, {{axis_pt_gen}, {axis_deta_fwd}}, true);
registry.add("StandaloneMuon/Ptgen_DeltaPhi_Pos", "resolution", kTH2F, {{axis_pt_gen}, {axis_dphi}}, true);
registry.add("StandaloneMuon/Ptgen_DeltaPhi_Neg", "resolution", kTH2F, {{axis_pt_gen}, {axis_dphi}}, true);
registry.addClone("StandaloneMuon/", "GlobalMuon/");
}

if (cfgFillTHnSparse) {
registry.add("Electron/hs_reso", "9D resolution", kTHnSparseF, {axis_cent, axis_pt_gen, axis_eta_cb_gen, axis_phi_gen, axis_phiPositionCB_gen, axis_charge_gen, axis_dpt, axis_deta_cb, axis_dphi}, true);
registry.add("StandaloneMuon/hs_reso", "9D resolution", kTHnSparseF, {axis_cent, axis_pt_gen, axis_eta_fwd_gen, axis_phi_gen, axis_phiPositionFWD_gen, axis_charge_gen, axis_dpt, axis_deta_fwd, axis_dphi}, true);
registry.add("GlobalMuon/hs_reso", "9D resolution", kTHnSparseF, {axis_cent, axis_pt_gen, axis_eta_fwd_gen, axis_phi_gen, axis_phiPositionFWD_gen, axis_charge_gen, axis_dpt, axis_deta_fwd, axis_dphi}, true);
registry.add("Electron/hs_reso", "8D resolution", kTHnSparseF, {axis_cent, axis_pt_gen, axis_eta_cb_gen, axis_phi_gen, axis_charge_gen, axis_dpt_cb, axis_deta_cb, axis_dphi}, true);
registry.add("StandaloneMuon/hs_reso", "8D resolution", kTHnSparseF, {axis_cent, axis_pt_gen, axis_eta_fwd_gen, axis_phi_gen, axis_charge_gen, axis_dpt_fwd, axis_deta_fwd, axis_dphi}, true);
registry.add("GlobalMuon/hs_reso", "8D resolution", kTHnSparseF, {axis_cent, axis_pt_gen, axis_eta_fwd_gen, axis_phi_gen, axis_charge_gen, axis_dpt_fwd, axis_deta_fwd, axis_dphi}, true);
}
}

Expand All @@ -298,10 +301,10 @@
}

// In case override, don't proceed, please - no CCDB access required
if (d_bz_input > -990) {

Check failure on line 304 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
d_bz = d_bz_input;
o2::parameters::GRPMagField grpmag;
if (std::fabs(d_bz) > 1e-5) {

Check failure on line 307 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
grpmag.setL3Current(30000.f / (d_bz / 5.0f));
}
o2::base::Propagator::initFieldFromGRP(&grpmag);
Expand Down Expand Up @@ -535,7 +538,7 @@
void fillMuon(TCollision const& collision, TMuon const& muon, TMFTTracksCov const& mftCovs, const float centrality)
{
auto mcparticle = muon.template mcParticle_as<aod::McParticles>();
if (std::abs(mcparticle.pdgCode()) != 13 || !(mcparticle.isPhysicalPrimary() || mcparticle.producedByGenerator())) {

Check failure on line 541 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
return;
}
if (cfg_require_true_mc_collision_association && mcparticle.mcCollisionId() != collision.mcCollisionId()) {
Expand Down Expand Up @@ -697,7 +700,7 @@
return;
}
if (cfgFillTHnSparse) {
registry.fill(HIST("StandaloneMuon/hs_reso"), centrality, mcparticle.pt(), mcparticle.eta(), mcparticle.phi(), 1.f, -mcparticle.pdgCode() / 13, (mcparticle.pt() - pt) / mcparticle.pt(), mcparticle.eta() - eta, mcparticle.phi() - phi);
registry.fill(HIST("StandaloneMuon/hs_reso"), centrality, mcparticle.pt(), mcparticle.eta(), mcparticle.phi(), -mcparticle.pdgCode() / 13, (mcparticle.pt() - pt) / mcparticle.pt(), mcparticle.eta() - eta, mcparticle.phi() - phi);
}

if (cfgFillTH2) {
Expand All @@ -705,9 +708,9 @@
registry.fill(HIST("StandaloneMuon/hEtaPhi"), phi, eta);
registry.fill(HIST("StandaloneMuon/Ptgen_RelDeltaPt"), mcparticle.pt(), (mcparticle.pt() - pt) / mcparticle.pt());
registry.fill(HIST("StandaloneMuon/Ptgen_DeltaEta"), mcparticle.pt(), mcparticle.eta() - eta);
if (mcparticle.pdgCode() == -13) { // positive muon

Check failure on line 711 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registry.fill(HIST("StandaloneMuon/Ptgen_DeltaPhi_Pos"), mcparticle.pt(), mcparticle.phi() - phi);
} else if (mcparticle.pdgCode() == 13) { // negative muon

Check failure on line 713 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registry.fill(HIST("StandaloneMuon/Ptgen_DeltaPhi_Neg"), mcparticle.pt(), mcparticle.phi() - phi);
}
}
Expand All @@ -716,16 +719,16 @@
return;
}
if (cfgFillTHnSparse) {
registry.fill(HIST("GlobalMuon/hs_reso"), centrality, mcparticle.pt(), mcparticle.eta(), mcparticle.phi(), 1.f, -mcparticle.pdgCode() / 13, (mcparticle.pt() - pt) / mcparticle.pt(), mcparticle.eta() - eta, mcparticle.phi() - phi);
registry.fill(HIST("GlobalMuon/hs_reso"), centrality, mcparticle.pt(), mcparticle.eta(), mcparticle.phi(), -mcparticle.pdgCode() / 13, (mcparticle.pt() - pt) / mcparticle.pt(), mcparticle.eta() - eta, mcparticle.phi() - phi);
}
if (cfgFillTH2) {
registry.fill(HIST("GlobalMuon/hPt"), pt);
registry.fill(HIST("GlobalMuon/hEtaPhi"), phi, eta);
registry.fill(HIST("GlobalMuon/Ptgen_RelDeltaPt"), mcparticle.pt(), (mcparticle.pt() - pt) / mcparticle.pt());
registry.fill(HIST("GlobalMuon/Ptgen_DeltaEta"), mcparticle.pt(), mcparticle.eta() - eta);
if (mcparticle.pdgCode() == -13) { // positive muon

Check failure on line 729 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registry.fill(HIST("GlobalMuon/Ptgen_DeltaPhi_Pos"), mcparticle.pt(), mcparticle.phi() - phi);
} else if (mcparticle.pdgCode() == 13) { // negative muon

Check failure on line 731 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registry.fill(HIST("GlobalMuon/Ptgen_DeltaPhi_Neg"), mcparticle.pt(), mcparticle.phi() - phi);
}
}
Expand Down Expand Up @@ -811,7 +814,7 @@
}
const auto& mcparticle = track.template mcParticle_as<aod::McParticles>();

if (std::abs(mcparticle.pdgCode()) != 11 || !(mcparticle.isPhysicalPrimary() || mcparticle.producedByGenerator())) {

Check failure on line 817 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
return;
}
if (cfg_reject_fake_match_its_tpc && o2::aod::pwgem::dilepton::utils::mcutil::hasFakeMatchITSTPC(track)) {
Expand Down Expand Up @@ -840,24 +843,24 @@
float phi = trackParCov.getPhi();
o2::math_utils::bringTo02Pi(phi);

float phiPosition = mcparticle.phi() + std::asin(-0.30282 * (-mcparticle.pdgCode() / 11) * (d_bz * 0.1) * cfgRefR / (2.f * mcparticle.pt()));
phiPosition = RecoDecay::constrainAngle(phiPosition, phiPosition_bin_edges[0], 1U);
// float phiPosition = mcparticle.phi() + std::asin(-0.30282 * (-mcparticle.pdgCode() / 11) * (d_bz * 0.1) * cfgRefR / (2.f * mcparticle.pt()));
// phiPosition = RecoDecay::constrainAngle(phiPosition, phiPosition_bin_edges[0], 1U);

if (!isSelectedTrackWithKine(track, pt, eta, trackParCov.getTgl(), dcaXY, dcaZ)) {
return;
}

if (cfgFillTHnSparse) {
registry.fill(HIST("Electron/hs_reso"), centrality, mcparticle.pt(), mcparticle.eta(), mcparticle.phi(), phiPosition, -mcparticle.pdgCode() / 11, (mcparticle.pt() - pt) / mcparticle.pt(), mcparticle.eta() - eta, mcparticle.phi() - phi);
registry.fill(HIST("Electron/hs_reso"), centrality, mcparticle.pt(), mcparticle.eta(), mcparticle.phi(), -mcparticle.pdgCode() / 11, (mcparticle.pt() - pt) / mcparticle.pt(), mcparticle.eta() - eta, mcparticle.phi() - phi);
}
if (cfgFillTH2) {
registry.fill(HIST("Electron/hPt"), pt);
registry.fill(HIST("Electron/hEtaPhi"), phi, eta);
registry.fill(HIST("Electron/Ptgen_RelDeltaPt"), mcparticle.pt(), (mcparticle.pt() - pt) / mcparticle.pt());
registry.fill(HIST("Electron/Ptgen_DeltaEta"), mcparticle.pt(), mcparticle.eta() - eta);
if (mcparticle.pdgCode() == -11) { // positron

Check failure on line 861 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registry.fill(HIST("Electron/Ptgen_DeltaPhi_Pos"), mcparticle.pt(), mcparticle.phi() - phi);
} else if (mcparticle.pdgCode() == 11) { // electron

Check failure on line 863 in PWGEM/Dilepton/Tasks/createResolutionMap.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
registry.fill(HIST("Electron/Ptgen_DeltaPhi_Neg"), mcparticle.pt(), mcparticle.phi() - phi);
}
}
Expand Down
Loading