Skip to content

[PWGEM/PhotonMeson] Add option for fast analytical propagation of tracks #17288

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
jokonig merged 1 commit into from
Aug 3, 2026
+194 −31
Merged

[PWGEM/PhotonMeson] Add option for fast analytical propagation of tracks #17288

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
142 changes: 111 additions & 31 deletions PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,12 @@
LUT = 2
};

enum TrackPropMode {
kProper = 0,
kFast = 1,
kBoth = 2
};

struct PhotonConversionBuilder {
Produces<aod::V0PhotonsKF> v0photonskf;
Produces<aod::V0Legs> v0legs;
Expand All @@ -114,6 +120,7 @@
// Operation and minimisation criteria
Configurable<double> d_bz_input{"d_bz", -999, "bz field, -999 is automatic"};
Configurable<int> useMatCorrType{"useMatCorrType", 0, "0: none, 1: TGeo, 2: LUT"};
Configurable<int> modeTrackPropagation{"modeTrackPropagation", 0, "0: use real track propagation, including material, 1: use fast approximation using only geometry, 2: Use real track propagation and make comparison to fast propagation (only for debugging and testing)"};

// single track cuts
Configurable<int> min_ncluster_tpc{"min_ncluster_tpc", 0, "min ncluster tpc"};
Expand Down Expand Up @@ -325,6 +332,14 @@
registry.add("V0/hBDTScoreAfterCutVsPt", "BDT score after cut vs pT; pT (GeV/c); BDT score", {HistType::kTH2F, {{1000, 0.0f, 20.0f}, {1000, 0.0f, 1.0f}}});
}
}

// Compare proper propagation and fast geometrical propagation
if (modeTrackPropagation == TrackPropMode::kBoth) {
registry.add("V0Leg/hDCAxyPropagationCompare", "Comparison of DCA_{xy} propagation;DCA_{xy} (cm) proper propagation; DCA_{xy} (cm) geom. propagation", {HistType::kTH2F, {{200, -10., 10.}, {200, -10., 10.}}});
registry.add("V0Leg/hDCAzPropagationCompare", "Comparison of DCA_{z} propagation;DCA_{z} (cm) proper propagation; DCA_{z} (cm) geom. propagation", {HistType::kTH2F, {{200, -10., 10.}, {200, -10., 10.}}});
registry.add("V0/hPhivPropagationCompare", "Comparison of #phi_{v};#phi_{v} proper propagation; #phi_{v} geom. propagation", {HistType::kTH2F, {{100, 0., 1.6}, {100, 0., 1.6}}});
registry.add("V0/hPsiPairPropagationCompare", "Comparison of #Psi_{pair};#Psi_{pair} proper propagation; #Psi_{pair} geom. propagation", {HistType::kTH2F, {{100, 0., 1.6}, {100, 0., 1.6}}});
}
}

void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
Expand All @@ -334,10 +349,10 @@
}

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

Check failure on line 352 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.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 355 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.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 @@ -368,7 +383,7 @@
}
mRunNumber = bc.runNumber();

if (useMatCorrType == 2) {

Check failure on line 386 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.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.
// setMatLUT only after magfield has been initalized (setMatLUT has implicit and problematic init field call if not)
o2::base::Propagator::Instance()->setMatLUT(lut);
}
Expand Down Expand Up @@ -555,8 +570,23 @@
return;
}
auto pTrackC = pTrack;
o2::math_utils::Point3D<float> vtxPrim{
collision.posX(),
collision.posY(),
collision.posZ()};
if (modeTrackPropagation != TrackPropMode::kProper) {
dcaInfo = CalculateDCAFast(pTrackC, vtxPrim, d_bz);
}
pTrackC.setPID(o2::track::PID::Electron);
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, pTrackC, 2.f, matCorr, &dcaInfo);

std::array<float, 2> dcaInfoFast = dcaInfo;
if (modeTrackPropagation != TrackPropMode::kFast) {
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, pTrackC, 2.f, matCorr, &dcaInfo);
if (modeTrackPropagation == TrackPropMode::kBoth) {
registry.fill(HIST("V0Leg/hDCAxyPropagationCompare"), dcaInfo[0], dcaInfoFast[0]);
registry.fill(HIST("V0Leg/hDCAzPropagationCompare"), dcaInfo[1], dcaInfoFast[1]);
}
}
auto posdcaXY = dcaInfo[0];
auto posdcaZ = dcaInfo[1];

Expand All @@ -566,8 +596,19 @@
return;
}
auto nTrackC = nTrack;
if (modeTrackPropagation != TrackPropMode::kProper) {
dcaInfo = CalculateDCAFast(nTrackC, vtxPrim, d_bz);
}
nTrackC.setPID(o2::track::PID::Electron);
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, nTrackC, 2.f, matCorr, &dcaInfo);

dcaInfoFast = dcaInfo;
if (modeTrackPropagation != TrackPropMode::kFast) {
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, nTrackC, 2.f, matCorr, &dcaInfo);
if (modeTrackPropagation == TrackPropMode::kBoth) {
registry.fill(HIST("V0Leg/hDCAxyPropagationCompare"), dcaInfo[0], dcaInfoFast[0]);
registry.fill(HIST("V0Leg/hDCAzPropagationCompare"), dcaInfo[1], dcaInfoFast[1]);
}
}
auto eledcaXY = dcaInfo[0];
auto eledcaZ = dcaInfo[1];

Expand All @@ -587,31 +628,70 @@

float phiv = 999.f;
float psipair = 999.f;
float phivFast = 999.f;
float psipairFast = 999.f;
float baseR = std::hypot(xyz[0], xyz[1]);
std::array<float, 3> offsetsR = {propV0LegsRadius, 30.f, 10.f};
bool pPropagatedSuccess = false;
bool nPropagatedSuccess = false;
auto pTrackProp = pTrack;
auto nTrackProp = nTrack;
for (const float& offsetR : offsetsR) {
pTrackProp = pTrack;
pTrackProp.setPID(o2::track::PID::Electron);
nTrackProp = nTrack;
nTrackProp.setPID(o2::track::PID::Electron);
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, pTrackProp, 2.f, matCorr, &dcaInfo);
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, nTrackProp, 2.f, matCorr, &dcaInfo);
pPropagatedSuccess = o2::base::Propagator::Instance()->propagateToR(pTrackProp, baseR + offsetR);
nPropagatedSuccess = o2::base::Propagator::Instance()->propagateToR(nTrackProp, baseR + offsetR);
if (pPropagatedSuccess && nPropagatedSuccess) {
KFPTrack kfp_track_posProp = createKFPTrackFromTrackParCov(pTrackProp, pos.sign(), pos.tpcNClsFound(), pos.tpcChi2NCl());
KFPTrack kfp_track_eleProp = createKFPTrackFromTrackParCov(nTrackProp, ele.sign(), ele.tpcNClsFound(), ele.tpcChi2NCl());
phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(kfp_track_posProp.GetPx(), kfp_track_posProp.GetPy(), kfp_track_posProp.GetPz(), kfp_track_eleProp.GetPx(), kfp_track_eleProp.GetPy(), kfp_track_eleProp.GetPz(), pos.sign(), ele.sign(), d_bz);
psipair = o2::aod::pwgem::dilepton::utils::pairutil::getPsiPair(kfp_track_posProp.GetPx(), kfp_track_posProp.GetPy(), kfp_track_posProp.GetPz(), kfp_track_eleProp.GetPx(), kfp_track_eleProp.GetPy(), kfp_track_eleProp.GetPz());
break;
// This method uses the track Helix instead of the full propagation.
// Hence, it is only an approximation but much faster
if (modeTrackPropagation != TrackPropMode::kProper) {

o2::track::TrackAuxPar helixPosEle(nTrack, d_bz);
o2::track::TrackAuxPar helixPosPos(pTrack, d_bz);

float diffX = helixPosEle.xC - helixPosPos.xC;
float diffY = helixPosEle.yC - helixPosPos.yC;
auto phiHelix = RecoDecay::constrainAngle<float, float>(std::atan2(diffY, diffX) - o2::constants::math::PI / 2.);

// Electron
float arcLenghtEle = helixPosEle.rC * 0.9 > propV0LegsRadius ? std::asin(propV0LegsRadius / helixPosEle.rC) * helixPosEle.rC : o2::constants::math::PI / 2.2 * helixPosEle.rC; // This assumes that the photon momentum vector is a tangent of the circle

Check failure on line 646 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.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.
auto propTrackEle = getPropMomentumFromTrackHelix(arcLenghtEle, ele, helixPosEle, d_bz / 10., phiHelix - ele.phi());
// Positron
float arcLenghtPos = helixPosPos.rC * 0.9 > propV0LegsRadius ? std::asin(propV0LegsRadius / helixPosPos.rC) * helixPosPos.rC : o2::constants::math::PI / 2.2 * helixPosPos.rC; // This assumes that the photon momentum vector is a tangent of the circle

Check failure on line 649 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.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.
auto propTrackPos = getPropMomentumFromTrackHelix(arcLenghtPos, pos, helixPosPos, d_bz / 10., phiHelix - pos.phi());

phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(propTrackPos[0], propTrackPos[1], propTrackPos[2], propTrackEle[0], propTrackEle[1], propTrackEle[2], pos.sign(), ele.sign(), d_bz);
psipair = o2::aod::pwgem::dilepton::utils::pairutil::getPsiPair(propTrackPos[0], propTrackPos[1], propTrackPos[2], propTrackEle[0], propTrackEle[1], propTrackEle[2]);

// Store values for later comparison
phivFast = phiv;
psipairFast = psipair;
}
// This uses the full propagation including material effects
if (modeTrackPropagation != TrackPropMode::kFast) {
std::array<float, 3> offsetsR = {propV0LegsRadius, 30.f, 10.f};
bool pPropagatedSuccess = false;
bool nPropagatedSuccess = false;
auto pTrackProp = pTrack;
auto nTrackProp = nTrack;
for (const auto& offsetR : offsetsR) {
pTrackProp = pTrack;
pTrackProp.setPID(o2::track::PID::Electron);
nTrackProp = nTrack;
nTrackProp.setPID(o2::track::PID::Electron);

o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, pTrackProp, 2.f, matCorr, &dcaInfo);
o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, nTrackProp, 2.f, matCorr, &dcaInfo);

pPropagatedSuccess = o2::base::Propagator::Instance()->propagateToR(pTrackProp, baseR + offsetR);
nPropagatedSuccess = o2::base::Propagator::Instance()->propagateToR(nTrackProp, baseR + offsetR);

if (pPropagatedSuccess && nPropagatedSuccess) {
KFPTrack kfp_track_posProp = createKFPTrackFromTrackParCov(pTrackProp, pos.sign(), pos.tpcNClsFound(), pos.tpcChi2NCl());
KFPTrack kfp_track_eleProp = createKFPTrackFromTrackParCov(nTrackProp, ele.sign(), ele.tpcNClsFound(), ele.tpcChi2NCl());
phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(kfp_track_posProp.GetPx(), kfp_track_posProp.GetPy(), kfp_track_posProp.GetPz(), kfp_track_eleProp.GetPx(), kfp_track_eleProp.GetPy(), kfp_track_eleProp.GetPz(), pos.sign(), ele.sign(), d_bz);
psipair = o2::aod::pwgem::dilepton::utils::pairutil::getPsiPair(kfp_track_posProp.GetPx(), kfp_track_posProp.GetPy(), kfp_track_posProp.GetPz(), kfp_track_eleProp.GetPx(), kfp_track_eleProp.GetPy(), kfp_track_eleProp.GetPz());
break;
} else {
LOG(debug) << "Propagation to offset" << offsetR << " cm failed for " << (pPropagatedSuccess ? "negative" : "positive") << " track. Trying smaller offset.";
}
}
if (modeTrackPropagation == TrackPropMode::kBoth) {
registry.fill(HIST("V0/hPhivPropagationCompare"), phiv, phivFast);
registry.fill(HIST("V0/hPsiPairPropagationCompare"), psipair, psipairFast);
}
LOG(debug) << "Propagation to offset" << offsetR << " cm failed for " << (pPropagatedSuccess ? "negative" : "positive") << " track. Trying smaller offset.";
}

if (phiv == 999.f || psipair == 999.f) {

Check failure on line 694 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.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.
LOG(debug) << "Propagation failed for all radii (" << propV0LegsRadius << ", 30, 10 cm). Using default values for phiv and psipair (999.f).";
}

Expand All @@ -624,7 +704,7 @@
KFParticle gammaKF;
gammaKF.SetConstructMethod(2);
gammaKF.Construct(GammaDaughters.data(), 2);
if (kfMassConstrain > -0.1) {

Check failure on line 707 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.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.
gammaKF.SetNonlinearMassConstraint(kfMassConstrain);
}
KFPVertex kfpVertex = createKFPVertexFromCollision(collision);
Expand Down Expand Up @@ -760,7 +840,7 @@
return;
}

if (v0photoncandidate.getChi2NDF() > 6e+3) { // protection for uint16.

Check failure on line 843 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.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.
return;
}

Expand Down Expand Up @@ -976,7 +1056,7 @@
auto v0 = v0s.rawIteratorAt(v0Id);
if constexpr (enableFilter) {
auto collision_tmp = v0.template collision_as<TCollisions>(); // collision where this v0 belongs.
if (!(collision_tmp.neeuls() >= 1 || collision_tmp.neeuls() + nv0_map[collision_tmp.globalIndex()] >= 2)) {

Check failure on line 1059 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.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.
continue;
}
// LOGF(info, "collision_tmp.globalIndex() = %d, collision_tmp.neeuls() = %d, nv0_map = %d", collision_tmp.globalIndex(), collision_tmp.neeuls(), nv0_map[collision_tmp.globalIndex()]);
Expand All @@ -985,14 +1065,14 @@
fillV0Table<isMC, TBCs, TCollisions, TTracks>(v0, true);
} // end of fullv0Id loop

for (const auto& collision : collisions) {
if constexpr (isMC) {
if (!collision.has_mcCollision()) {
continue;
}
}
// events_ngpcm(nv0_map[collision.globalIndex()]);
} // end of collision loop
// for (const auto& collision : collisions) {
// if constexpr (isMC) {
// if (!collision.has_mcCollision()) {
// continue;
// }
// }
// // events_ngpcm(nv0_map[collision.globalIndex()]);
// } // end of collision loop

pca_map.clear();
cospa_map.clear();
Expand Down Expand Up @@ -1041,5 +1121,5 @@
WorkflowSpec defineDataProcessing(ConfigContext const& context)
{
return WorkflowSpec{
adaptAnalysisTask<PhotonConversionBuilder>(context, TaskName{"photon-conversion-builder"})};

Check failure on line 1124 in PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/o2-task] 

Specified task name photon-conversion-builder and the struct name PhotonConversionBuilder produce the same device name photon-conversion-builder. TaskName is redundant.
}
83 changes: 83 additions & 0 deletions PWGEM/PhotonMeson/Utils/PCMUtilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@

#include <CommonConstants/MathConstants.h>
#include <DetectorsBase/Propagator.h>
#include <Framework/ASoA.h>
#include <Framework/Concepts.h>
#include <ReconstructionDataFormats/HelixHelper.h>
#include <ReconstructionDataFormats/TrackParametrizationWithError.h>
Expand Down Expand Up @@ -101,10 +102,92 @@ inline void Vtx_recalculationParCov(o2::base::Propagator* prop, const o2::track:

xyz[2] = (trackPosInformationCopy.getZ() * helixNeg.rC + trackNegInformationCopy.getZ() * helixPos.rC) / (helixPos.rC + helixNeg.rC);
}

//_______________________________________________________________________
/// \brief Calculate DCA for tracks using the track helix and the inclination angl tan(lambda)
/// \param trk track parameterization to obtain helix parameters
/// \param vtx primary vertex position
/// \param magField magnetic field strenght of L3
/// \return DCAxy, DCAz
template <typename TrackPrecision = float>
std::array<float, 2> CalculateDCAFast(const o2::track::TrackParametrizationWithError<TrackPrecision>& trk, const o2::math_utils::Point3D<float>& vtx, const float magField)
{

std::array<float, 2> dca;

// obtain circle from track in x-y plane
const o2::track::TrackAuxPar helixPos(trk, magField);

// obtain position in global coordinates and tan(lambda)
const auto posTrack = trk.getXYZGlo();
const float trX = posTrack.X();
const float trY = posTrack.Y();
const float trZ = posTrack.Z();
const float tangentLambda = trk.getTgl();

// Calculate DCAxy
// Use distance in x and y between circle center and vtx, afterwards subtract radius of circle
const float distX = helixPos.xC - vtx.X();
const float distY = helixPos.yC - vtx.Y();
const float trackCircCenter = std::sqrt(distX * distX + distY * distY);
dca[0] = helixPos.rC - trackCircCenter;

// Calculate DCAz
// First step: Calculate arc lenght of circle between current position and primary vertex in x-y
const float theta0 = std::atan2(trY - helixPos.yC, trX - helixPos.xC);
const float thetav = std::atan2(vtx.Y() - helixPos.yC, vtx.X() - helixPos.xC);

// Make sure angle is between -pi and pi
const auto dtheta = RecoDecay::constrainAngle<float>(thetav - theta0, -o2::constants::math::PI);

// arc-lenght along helix
const float arcLenght = std::fabs(helixPos.rC * dtheta);

// get global z-position at DCA
const float ZPosGlo = trZ - arcLenght * tangentLambda;

// DCA calculated from
dca[1] = ZPosGlo - vtx.Z();

return dca;
}

//_______________________________________________________________________
/// \brief Calculate the track momentum at a different place of the track Helix.
/// \param s arc-lenght where track should be propagated to
/// \param track particle track
/// \param trHelix track helix param. for circle approximation
/// \param bz magnetic field strenght of L3
/// \param addPhi optional additional rotation of the track in phi-direction
/// \return track momentum vector at new position
template <o2::soa::is_iterator TTrack>
inline std::array<float, 3> getPropMomentumFromTrackHelix(const float s, const TTrack& track, const o2::track::TrackAuxPar& trHelix, const float bz, float addPhi = 0.f)
{

// Calculate the change in phi considering the track radius and the track arc lenght s
const float dphi = -track.sign() * 0.3f * bz * s / 100. / track.pt(); // s in cm
const float dotProd = std::cos(track.phi() + addPhi) * track.pt() * trHelix.xC + std::sin(track.phi() + addPhi) * track.pt() * trHelix.yC;
if (dotProd < 0) {
addPhi -= o2::constants::math::PI;
}

// Calculate the phi at the secondary vertex
const auto phi = RecoDecay::constrainAngle<float>(track.phi() + dphi + addPhi);

// Calculate px,y,z at the new propagated vertex
std::array<float, 3> trackP;
trackP[0] = std::cos(phi) * track.pt();
trackP[1] = std::sin(phi) * track.pt();
trackP[2] = track.tgl() * track.pt();

return trackP;
}

//_______________________________________________________________________
template <typename TrackPrecision = float, o2::soa::is_iterator T1, o2::soa::is_iterator T2>
inline void Vtx_recalculation(o2::base::Propagator* prop, T1 lTrackPos, T2 lTrackNeg, std::array<float, 3>& xyz, o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE)
{

// o2::track::TrackParametrizationWithError<TrackPrecision> = TrackParCov, I use the full version to have control over the data type
o2::track::TrackParametrizationWithError<TrackPrecision> trackPosInformation = getTrackParCov(lTrackPos); // first get an object that stores Track information (positive)
o2::track::TrackParametrizationWithError<TrackPrecision> trackNegInformation = getTrackParCov(lTrackNeg); // first get an object that stores Track information (negative)
Expand Down
Loading