Skip to content

[PWGEM/Dilepton,Common] update global muon propagation #13881

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
Nov 19, 2025
Merged

[PWGEM/Dilepton,Common] update global muon propagation #13881

Show file tree
Hide file tree
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
30 changes: 15 additions & 15 deletions Common/Core/fwdtrackUtilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,17 +18,17 @@
#ifndef COMMON_CORE_FWDTRACKUTILITIES_H_
#define COMMON_CORE_FWDTRACKUTILITIES_H_

#include "DetectorsBase/GeometryManager.h"
#include "Field/MagneticField.h"
#include "Framework/AnalysisDataModel.h"
#include <DetectorsBase/GeometryManager.h>
#include <Field/MagneticField.h>
#include <GlobalTracking/MatchGlobalFwd.h>
#include <MCHTracking/TrackExtrap.h>
#include <ReconstructionDataFormats/GlobalFwdTrack.h>
#include <ReconstructionDataFormats/TrackFwd.h>
#include "GlobalTracking/MatchGlobalFwd.h"
#include "MCHTracking/TrackExtrap.h"
#include "ReconstructionDataFormats/GlobalFwdTrack.h"
#include "ReconstructionDataFormats/TrackFwd.h"

#include <Math/MatrixRepresentationsStatic.h>
#include <Math/SMatrix.h>
#include <TGeoGlobalMagField.h>
#include "Math/MatrixRepresentationsStatic.h"
#include "Math/SMatrix.h"
#include "TGeoGlobalMagField.h"

#include <type_traits>
#include <vector>
Expand All @@ -51,8 +51,7 @@
template <typename TFwdTrack, typename TFwdTrackCov>
o2::track::TrackParCovFwd getTrackParCovFwd(TFwdTrack const& track, TFwdTrackCov const& cov)
{
// This function works for both (saMuon, saMuon) and (MFTTrack, MFTTrackCov).
// Don't use covariant matrix of global muons stored in AO2D.root.
// This function works for (glMuon, glMuon), (saMuon, saMuon) and (MFTTrack, MFTTrackCov).

double chi2 = track.chi2();
if constexpr (std::is_same_v<std::decay_t<TFwdTrackCov>, aod::MFTTracksCov::iterator>) {
Expand Down Expand Up @@ -102,7 +101,7 @@
o2::dataformats::GlobalFwdTrack propmuon;
o2::globaltracking::MatchGlobalFwd mMatching;

if (trackType > 2) { // MCH-MID or MCH standalone

Check failure on line 104 in Common/Core/fwdtrackUtilities.h

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.
o2::dataformats::GlobalFwdTrack track;
track.setParameters(fwdtrack.getParameters());
track.setZ(fwdtrack.getZ());
Expand All @@ -123,16 +122,17 @@
propmuon.setParameters(proptrack.getParameters());
propmuon.setZ(proptrack.getZ());
propmuon.setCovariances(proptrack.getCovariances());
} else if (trackType < 2) { // MFT-MCH-MID

Check failure on line 125 in Common/Core/fwdtrackUtilities.h

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.
// const double centerMFT[3] = {0, 0, -61.4};
// o2::field::MagneticField* field = static_cast<o2::field::MagneticField*>(TGeoGlobalMagField::Instance()->GetField());
// auto Bz = field->getBz(centerMFT); // Get field at centre of MFT in kG.

auto geoMan = o2::base::GeometryManager::meanMaterialBudget(fwdtrack.getX(), fwdtrack.getY(), fwdtrack.getZ(), collision.posX(), collision.posY(), collision.posZ());
auto x2x0 = static_cast<float>(geoMan.meanX2X0);

if (endPoint == propagationPoint::kToVertex) {
fwdtrack.propagateToVtxhelixWithMCS(collision.posZ(), {collision.posX(), collision.posY()}, {collision.covXX(), collision.covYY()}, bzkG, x2x0);
// auto geoMan = o2::base::GeometryManager::meanMaterialBudget(fwdtrack.getX(), fwdtrack.getY(), fwdtrack.getZ(), collision.posX(), collision.posY(), collision.posZ());
// auto x2x0 = static_cast<float>(geoMan.meanX2X0);
// fwdtrack.propagateToVtxhelixWithMCS(collision.posZ(), {collision.posX(), collision.posY()}, {collision.covXX(), collision.covYY()}, bzkG, x2x0);
std::array<double, 3> dcaInfOrig{999.f, 999.f, 999.f};
fwdtrack.propagateToDCAhelix(bzkG, {collision.posX(), collision.posY(), collision.posZ()}, dcaInfOrig);
} else if (endPoint == propagationPoint::kToDCA) {
fwdtrack.propagateToZhelix(collision.posZ(), bzkG);
} else if (endPoint == propagationPoint::kToMatchingPlane) {
Expand Down
205 changes: 69 additions & 136 deletions PWGEM/Dilepton/Core/Dilepton.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,14 +62,14 @@
#include <utility>
#include <vector>

using namespace o2;

Check failure on line 65 in PWGEM/Dilepton/Core/Dilepton.h

View workflow job for this annotation

GitHub Actions / O2 linter

[using-directive] 

Do not put using directives at global scope in headers.
using namespace o2::aod;

Check failure on line 66 in PWGEM/Dilepton/Core/Dilepton.h

View workflow job for this annotation

GitHub Actions / O2 linter

[using-directive] 

Do not put using directives at global scope in headers.
using namespace o2::framework;

Check failure on line 67 in PWGEM/Dilepton/Core/Dilepton.h

View workflow job for this annotation

GitHub Actions / O2 linter

[using-directive] 

Do not put using directives at global scope in headers.
using namespace o2::framework::expressions;

Check failure on line 68 in PWGEM/Dilepton/Core/Dilepton.h

View workflow job for this annotation

GitHub Actions / O2 linter

[using-directive] 

Do not put using directives at global scope in headers.
using namespace o2::soa;

Check failure on line 69 in PWGEM/Dilepton/Core/Dilepton.h

View workflow job for this annotation

GitHub Actions / O2 linter

[using-directive] 

Do not put using directives at global scope in headers.
using namespace o2::aod::pwgem::dilepton::utils;

Check failure on line 70 in PWGEM/Dilepton/Core/Dilepton.h

View workflow job for this annotation

GitHub Actions / O2 linter

[using-directive] 

Do not put using directives at global scope in headers.
using namespace o2::aod::pwgem::dilepton::utils::emtrackutil;

Check failure on line 71 in PWGEM/Dilepton/Core/Dilepton.h

View workflow job for this annotation

GitHub Actions / O2 linter

[using-directive] 

Do not put using directives at global scope in headers.
using namespace o2::aod::pwgem::dilepton::utils::pairutil;

Check failure on line 72 in PWGEM/Dilepton/Core/Dilepton.h

View workflow job for this annotation

GitHub Actions / O2 linter

[using-directive] 

Do not put using directives at global scope in headers.

using MyCollisions = soa::Join<aod::EMEvents, aod::EMEventsMult, aod::EMEventsCent, aod::EMEventsQvec>;
using MyCollision = MyCollisions::iterator;
Expand Down Expand Up @@ -418,111 +418,84 @@

emh_pos = new TEMH(ndepth);
emh_neg = new TEMH(ndepth);
emh_pair_uls = new MyEMH_pair(ndepth);
emh_pair_lspp = new MyEMH_pair(ndepth);
emh_pair_lsmm = new MyEMH_pair(ndepth);

if (accBins.ConfMllAccBins.value[0] == VARIABLE_WIDTH) {
mll_bin_edges = std::vector<float>(accBins.ConfMllAccBins.value.begin(), accBins.ConfMllAccBins.value.end());
mll_bin_edges.erase(mll_bin_edges.begin());
for (const auto& edge : mll_bin_edges) {
LOGF(info, "VARIABLE_WIDTH: mll_bin_edges = %f", edge);
}
} else {
int nbins = static_cast<int>(accBins.ConfMllAccBins.value[0]);
float xmin = static_cast<float>(accBins.ConfMllAccBins.value[1]);
float xmax = static_cast<float>(accBins.ConfMllAccBins.value[2]);
mll_bin_edges.resize(nbins + 1);
for (int i = 0; i < nbins + 1; i++) {
mll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin;
LOGF(info, "FIXED_WIDTH: mll_bin_edges[%d] = %f", i, mll_bin_edges[i]);
}
}
if (cfgAnalysisType == static_cast<int>(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kPolarization)) { // only for polarization

if (accBins.ConfPtllAccBins.value[0] == VARIABLE_WIDTH) {
ptll_bin_edges = std::vector<float>(accBins.ConfPtllAccBins.value.begin(), accBins.ConfPtllAccBins.value.end());
ptll_bin_edges.erase(ptll_bin_edges.begin());
for (const auto& edge : ptll_bin_edges) {
LOGF(info, "VARIABLE_WIDTH: ptll_bin_edges = %f", edge);
}
} else {
int nbins = static_cast<int>(accBins.ConfPtllAccBins.value[0]);
float xmin = static_cast<float>(accBins.ConfPtllAccBins.value[1]);
float xmax = static_cast<float>(accBins.ConfPtllAccBins.value[2]);
ptll_bin_edges.resize(nbins + 1);
for (int i = 0; i < nbins + 1; i++) {
ptll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin;
LOGF(info, "FIXED_WIDTH: ptll_bin_edges[%d] = %f", i, ptll_bin_edges[i]);
}
}
emh_pair_uls = new MyEMH_pair(ndepth);
emh_pair_lspp = new MyEMH_pair(ndepth);
emh_pair_lsmm = new MyEMH_pair(ndepth);

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

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

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

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

int nM = mll_bin_edges.size();
int nPt = ptll_bin_edges.size();
int nEta = etall_bin_edges.size();
int nPhi = phill_bin_edges.size();

// emhs_pair_uls.resize(nM);
// emhs_pair_lspp.resize(nM);
// emhs_pair_lsmm.resize(nM);
// for (int im = 0;im<nM;im++) {
// emhs_pair_uls[im].resize(nPt);
// emhs_pair_lspp[im].resize(nPt);
// emhs_pair_lsmm[im].resize(nPt);
// for (int ipt = 0;ipt<nPt;ipt++) {
// emhs_pair_uls[im][ipt].resize(nEta);
// emhs_pair_lspp[im][ipt].resize(nEta);
// emhs_pair_lsmm[im][ipt].resize(nEta);
// for (int ieta = 0;ieta<nEta;ieta++) {
// emhs_pair_uls[im][ipt][ieta].resize(nPhi);
// emhs_pair_lspp[im][ipt][ieta].resize(nPhi);
// emhs_pair_lsmm[im][ipt][ieta].resize(nPhi);
// for (int iphi = 0;iphi<nPhi;iphi++) {
// emhs_pair_uls[im][ipt][ieta][iphi] = new MyEMH_pair(ndepth);
// emhs_pair_lspp[im][ipt][ieta][iphi] = new MyEMH_pair(ndepth);
// emhs_pair_lsmm[im][ipt][ieta][iphi] = new MyEMH_pair(ndepth);
// }
// }
// }
// }

LOGF(info, "nM = %d, nPt = %d, nEta = %d, nPhi = %d", nM, nPt, nEta, nPhi);

std::random_device seed_gen;
engine = std::mt19937(seed_gen());
std::random_device seed_gen;
engine = std::mt19937(seed_gen());
}

DefineEMEventCut();
addhistograms();
Expand Down Expand Up @@ -637,30 +610,6 @@
delete emh_pair_lsmm;
emh_pair_lsmm = 0x0;

// for (int im = 0;im<emhs_pair_uls.size();im++) {
// for (int ipt = 0;ipt<emhs_pair_uls[im].size();ipt++) {
// for (int ieta = 0;ieta<emhs_pair_uls[im][ipt].size();ieta++) {
// for (int iphi = 0;iphi<emhs_pair_uls[im][ipt][ieta].size();iphi++) {
// delete emhs_pair_uls[im][ipt][ieta][iphi];
// emhs_pair_uls[im][ipt][ieta][iphi] = 0x0;

// delete emhs_pair_lspp[im][ipt][ieta][iphi];
// emhs_pair_lspp[im][ipt][ieta][iphi] = 0x0;

// delete emhs_pair_lsmm[im][ipt][ieta][iphi];
// emhs_pair_lsmm[im][ipt][ieta][iphi] = 0x0;
// }
// }
// }
// }

// emhs_pair_uls.clear();
// emhs_pair_uls.shrink_to_fit();
// emhs_pair_lspp.clear();
// emhs_pair_lspp.shrink_to_fit();
// emhs_pair_lsmm.clear();
// emhs_pair_lsmm.shrink_to_fit();

used_trackIds_per_col.clear();
used_trackIds_per_col.shrink_to_fit();
map_mixed_eventId_to_globalBC.clear();
Expand Down Expand Up @@ -1348,9 +1297,6 @@
MyEMH_pair* emh_pair_lspp = nullptr;
MyEMH_pair* emh_pair_lsmm = nullptr;

// std::vector<std::vector<std::vector<std::vector<MyEMH_pair*>>>> emhs_pair_uls; // 4D{m, pt, eta, phi}
// std::vector<std::vector<std::vector<std::vector<MyEMH_pair*>>>> emhs_pair_lspp; // 4D{m, pt, eta, phi}
// std::vector<std::vector<std::vector<std::vector<MyEMH_pair*>>>> emhs_pair_lsmm; // 4D{m, pt, eta, phi}
std::map<std::pair<int, int>, uint64_t> map_mixed_eventId_to_globalBC;

std::vector<int> used_trackIds_per_col;
Expand Down Expand Up @@ -1730,19 +1676,6 @@
emh_pair_uls->AddCollisionIdAtLast(key_bin, key_df_collision);
emh_pair_lspp->AddCollisionIdAtLast(key_bin, key_df_collision);
emh_pair_lsmm->AddCollisionIdAtLast(key_bin, key_df_collision);

// for (int im = 0;im<emhs_pair_uls.size();im++) {
// for (int ipt = 0;ipt<emhs_pair_uls[im].size();ipt++) {
// for (int ieta = 0;ieta<emhs_pair_uls[im][ipt].size();ieta++) {
// for (int iphi = 0;iphi<emhs_pair_uls[im][ipt][ieta].size();iphi++) {
// emhs_pair_uls[im][ipt][ieta][iphi]->AddCollisionIdAtLast(key_bin, key_df_collision);
// emhs_pair_lspp[im][ipt][ieta][iphi]->AddCollisionIdAtLast(key_bin, key_df_collision);
// emhs_pair_lsmm[im][ipt][ieta][iphi]->AddCollisionIdAtLast(key_bin, key_df_collision);
// }
// }
// }
// }

} // end of if pair exist

} // end of collision loop
Expand Down
4 changes: 2 additions & 2 deletions PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -181,8 +181,8 @@ struct CreateEMEventDilepton {
q3xft0m = collision.qvecFT0MReVec()[1], q3xft0a = collision.qvecFT0AReVec()[1], q3xft0c = collision.qvecFT0CReVec()[1], q3xfv0a = collision.qvecFV0AReVec()[1], q3xbpos = collision.qvecBPosReVec()[1], q3xbneg = collision.qvecBNegReVec()[1], q3xbtot = collision.qvecBTotReVec()[1];
q3yft0m = collision.qvecFT0MImVec()[1], q3yft0a = collision.qvecFT0AImVec()[1], q3yft0c = collision.qvecFT0CImVec()[1], q3yfv0a = collision.qvecFV0AImVec()[1], q3ybpos = collision.qvecBPosImVec()[1], q3ybneg = collision.qvecBNegImVec()[1], q3ybtot = collision.qvecBTotImVec()[1];
} else if (collision.qvecFT0CReVec().size() >= 1) { // harmonics 2
q2xft0m = collision.qvecFT0MReVec()[0], q2xft0a = collision.qvecFT0AReVec()[0], q2xft0c = collision.qvecFT0CReVec()[0], q2xbpos = collision.qvecBPosReVec()[0], q2xbneg = collision.qvecBNegReVec()[0], q2xbtot = collision.qvecBTotReVec()[0];
q2yft0m = collision.qvecFT0MImVec()[0], q2yft0a = collision.qvecFT0AImVec()[0], q2yft0c = collision.qvecFT0CImVec()[0], q2ybpos = collision.qvecBPosImVec()[0], q2ybneg = collision.qvecBNegImVec()[0], q2ybtot = collision.qvecBTotImVec()[0];
q2xft0m = collision.qvecFT0MReVec()[0], q2xft0a = collision.qvecFT0AReVec()[0], q2xft0c = collision.qvecFT0CReVec()[0], q2xfv0a = collision.qvecFV0AReVec()[0], q2xbpos = collision.qvecBPosReVec()[0], q2xbneg = collision.qvecBNegReVec()[0], q2xbtot = collision.qvecBTotReVec()[0];
q2yft0m = collision.qvecFT0MImVec()[0], q2yft0a = collision.qvecFT0AImVec()[0], q2yft0c = collision.qvecFT0CImVec()[0], q2yfv0a = collision.qvecFV0AImVec()[0], q2ybpos = collision.qvecBPosImVec()[0], q2ybneg = collision.qvecBNegImVec()[0], q2ybtot = collision.qvecBTotImVec()[0];
}
event_qvec(
q2xft0m, q2yft0m, q2xft0a, q2yft0a, q2xft0c, q2yft0c, q2xfv0a, q2yfv0a, q2xbpos, q2ybpos, q2xbneg, q2ybneg, q2xbtot, q2ybtot,
Expand Down
Loading
Loading