Skip to content

[ALICE3] Implement basic bremsstrahlung model in OTF tracker #16417

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

Open
axti98 wants to merge 2 commits into
base: master
Choose a base branch
from
+85 −0
Open

[ALICE3] Implement basic bremsstrahlung model in OTF tracker #16417

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a65e063
Implement bremsstrahlung in OTF tracker
axti98 May 22, 2026
3616a2d
Add basic QA histograms for BR
axti98 May 27, 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
85 changes: 85 additions & 0 deletions ALICE3/TableProducer/OTF/onTheFlyTracker.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,7 @@
#include <TGeoGlobalMagField.h>
#include <TH1.h>
#include <TH2.h>
#include <TLorentzVector.h>

Check failure on line 81 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
#include <TMCProcess.h>
#include <TMath.h>
#include <TPDGCode.h>
Expand Down Expand Up @@ -193,6 +193,9 @@

ConfigurableAxis axisRadius{"axisRadius", {2500, 0.0f, +250.0f}, "R (cm)"};
ConfigurableAxis axisZ{"axisZ", {100, -250.0f, +250.0f}, "Z (cm)"};

ConfigurableAxis axisNphotons{"axisNphotons", {10, 0.5f, 10.5f}, "N_{#gamma}"};
ConfigurableAxis axisBRenergyLoss{"axisBRenergyLoss", {500, 0.0f, 1.0f}, "#Delta p / p"};
} axes;

// for topo var QA
Expand Down Expand Up @@ -246,6 +249,15 @@
Configurable<bool> useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"};
} cfgFitter;

struct : ConfigurableGroup {
std::string prefix = "cfgBR"; // configuration for bremsstrahlung
Configurable<bool> radiateBR{"radiateBR", false, "simulate bremsstrahlung"};
Configurable<bool> doBRQA{"doBRQA", false, "Do QA for bremsstrahlung"};
Configurable<float> minBREnergyFraction{"minEnergyFraction", 0.001f, "Minimum energy fraction a bremsstrahlung photon can carry"};
Configurable<float> maxBREnergyFraction{"maxEnergyFraction", 0.95f, "Maximum energy fraction a bremsstrahlung photon can carry"};
Configurable<float> radiationStrength{"radiationStrength", 1e-6f, "Strenght of the bremsstrahlung radiation"};
} brSettings;

using PVertex = o2::dataformats::PrimaryVertex;

// for secondary vertex finding
Expand Down Expand Up @@ -568,6 +580,14 @@
insertHist(histPath + "h2dDCAz", "h2dDCAz;p_{T};DCA_{z}", {kTH2D, {{axes.axisMomentum, axes.axisDCA}}});
}

if (brSettings.doBRQA) {
insertHist(histPath + "h1dNBRPhotons", "h1dNBRPhotons;N_{#gamma};Counts", {kTH1D, {{axes.axisNphotons}}});
insertHist(histPath + "h1dBREnergyLoss", "h1dBREnergyLoss;#Delta p / p;Counts", {kTH1D, {{axes.axisBRenergyLoss}}});

insertHist(histPath + "h2dBRPtRes", "h2dPtRes;Gen p_{T};#Delta p_{T} / Reco p_{T}", {kTH2D, {{axes.axisMomentum, axes.axisPtRes}}});
insertHist(histPath + "h2dBRPtResAbs", "h2dPtResAbs;Gen p_{T};#Delta p_{T}", {kTH2D, {{axes.axisMomentum, axes.axisPtRes}}});
}

} // end config loop

// Basic QA
Expand Down Expand Up @@ -744,7 +764,7 @@
/// \param xiDecayVertex the address of the xi decay vertex
/// \param laDecayVertex the address of the la decay vertex
template <typename McParticleType>
void decayCascade(McParticleType particle, o2::track::TrackParCov track, std::vector<TLorentzVector>& decayDaughters, std::vector<double>& xiDecayVertex, std::vector<double>& laDecayVertex)

Check failure on line 767 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
{
const double uXi = rand.Uniform(0, 1);
const double ctauXi = 4.91; // cm
Expand All @@ -767,12 +787,12 @@
xiDecayVertex.push_back(particle.vz() + rxyzXi * (particle.pz() / particle.p()));

std::vector<double> xiDaughters = {o2::constants::physics::MassLambda, o2::constants::physics::MassPionCharged};
TLorentzVector xi(newPx, newPy, particle.pz(), newE);

Check failure on line 790 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
TGenPhaseSpace xiDecay;
xiDecay.SetDecay(xi, 2, xiDaughters.data());
xiDecay.Generate();
decayDaughters.push_back(*xiDecay.GetDecay(1));
TLorentzVector la = *xiDecay.GetDecay(0);

Check failure on line 795 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.

const double uLa = rand.Uniform(0, 1);
const double ctauLa = 7.845; // cm
Expand All @@ -795,7 +815,7 @@
/// \param decayDaughters the address of resulting daughters
/// \param v0DecayVertex the address of the la decay vertex
template <typename McParticleType>
void decayV0Particle(McParticleType particle, std::vector<TLorentzVector>& decayDaughters, std::vector<double>& v0DecayVertex, int pdgCode)

Check failure on line 818 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
{
double u = rand.Uniform(0, 1);
double v0Mass = -1.;
Expand Down Expand Up @@ -829,7 +849,7 @@

const double v0BetaGamma = particle.p() / v0Mass;
const double v0rxyz = (-v0BetaGamma * ctau * std::log(1 - u));
TLorentzVector v0(particle.px(), particle.py(), particle.pz(), particle.e());

Check failure on line 852 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.

v0DecayVertex.push_back(particle.vx() + v0rxyz * (particle.px() / particle.p()));
v0DecayVertex.push_back(particle.vy() + v0rxyz * (particle.py() / particle.p()));
Expand Down Expand Up @@ -904,7 +924,7 @@
o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB);
const std::string histPath = "Configuration_" + std::to_string(icfg) + "/";

std::vector<TLorentzVector> cascadeDecayProducts;

Check failure on line 927 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
std::vector<double> xiDecayVertex, laDecayVertex;
static constexpr int kCascProngs = 3;
std::array<o2::track::TrackParCov, kCascProngs> xiDaughterTrackParCovsPerfect;
Expand Down Expand Up @@ -934,11 +954,11 @@
getHist(TH1, histPath + "hXiBuilding")->Fill(0.0f);
}

o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[0], xiDecayVertex, xiDaughterTrackParCovsPerfect[0], pdgDB);

Check failure on line 957 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
xiDaughterTrackParCovsPerfect[0].setPID(pdgCodeToPID(PDG_t::kPiMinus));
o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kPiMinus, cascadeDecayProducts[1], laDecayVertex, xiDaughterTrackParCovsPerfect[1], pdgDB);

Check failure on line 959 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
xiDaughterTrackParCovsPerfect[1].setPID(pdgCodeToPID(PDG_t::kPiMinus));
o2::upgrade::convertTLorentzVectorToO2Track(PDG_t::kProton, cascadeDecayProducts[2], laDecayVertex, xiDaughterTrackParCovsPerfect[2], pdgDB);

Check failure on line 961 in ALICE3/TableProducer/OTF/onTheFlyTracker.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
xiDaughterTrackParCovsPerfect[2].setPID(pdgCodeToPID(PDG_t::kProton));
o2::track::TrackParCov perfectCascadeTrack;
o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectCascadeTrack, pdgDB);
Expand Down Expand Up @@ -1749,6 +1769,67 @@
}
}

/// Function to compute the bremsstrahlung loss of charged-particles for each layer of the current configuration
/// \param icfg index of the current configuration
/// \param mcParticle true MC particle to identify particle and get the energy
/// \param trackParCov track of the particle to compute bremsstrahlung for
void computeBremsstrahlungLoss(const int icfg, const auto& mcParticle, o2::track::TrackParCov& trackParCov)
{
if (brSettings.radiateBR) {
const o2::fastsim::GeometryEntry geoEntry = mGeoContainer.getEntry(icfg);

for (auto const& layerName : geoEntry.getLayerNames()) {
if (layerName.find("global") != std::string::npos) { // Layers with global tag are skipped
continue;
}

float mass = o2::constants::physics::MassElectron;

switch (std::abs(mcParticle.pdgCode())) {
case kElectron:
mass = o2::constants::physics::MassElectron;
break;
case kMuonMinus:
mass = o2::constants::physics::MassMuon;
break;
case kPiPlus:
mass = o2::constants::physics::MassPionCharged;
break;
case kKPlus:
mass = o2::constants::physics::MassKaonCharged;
break;
case kProton:
mass = o2::constants::physics::MassProton;
break;
default:
break;
}

float lambda = brSettings.radiationStrength * mcParticle.e() * geoEntry.getFloatValue(layerName, "x0") / (mass * mass);
ULong64_t nPhotons = gRandom->Poisson(lambda);

double initialMomentum = trackParCov.getP();

for (ULong64_t photon = 0; photon < nPhotons; ++photon) {
float radiativeLoss = 1.0f - brSettings.minBREnergyFraction * std::pow(brSettings.maxBREnergyFraction / brSettings.minBREnergyFraction, gRandom->Rndm());
trackParCov.setQ2Pt(trackParCov.getQ2Pt() / radiativeLoss);
}

double afterRadiationMomentum = trackParCov.getP();

if (brSettings.doBRQA) {
const std::string histPath = "Configuration_" + std::to_string(icfg) + "/";

getHist(TH1, histPath + "h1dNBRPhotons")->Fill(static_cast<double>(nPhotons));
getHist(TH1, histPath + "h1dBREnergyLoss")->Fill((initialMomentum - afterRadiationMomentum) / afterRadiationMomentum);

getHist(TH2, histPath + "h2dBRPtRes")->Fill(trackParCov.getPt(), (trackParCov.getPt() - mcParticle.pt()) / trackParCov.getPt());
getHist(TH2, histPath + "h2dBRPtResAbs")->Fill(trackParCov.getPt(), trackParCov.getPt() - mcParticle.pt());
}
}
}
}

void processWithLUTs(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const int icfg)
{
const std::string histPath = "Configuration_" + std::to_string(icfg) + "/";
Expand Down Expand Up @@ -1811,12 +1892,14 @@
o2::track::TrackParCov perfectTrackParCov;
o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectTrackParCov, pdgDB);
perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov);
nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta);
if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) {
reconstructed = false;
}
} else {
o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB);
computeBremsstrahlungLoss(icfg, mcParticle, trackParCov);
reconstructed = mSmearer[icfg]->smearTrack(trackParCov, mcParticle.pdgCode(), dNdEta);
nTrkHits = fastTrackerSettings.minSiliconHits;
}
Expand Down Expand Up @@ -2002,11 +2085,13 @@
int nTrkHits = 0;
if (enablePrimarySmearing && mcParticle.isPrimary()) {
o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB);
computeBremsstrahlungLoss(icfg, mcParticle, trackParCov);
reconstructed = mSmearer[icfg]->smearTrack(trackParCov, mcParticle.pdgCode(), dNdEta);
nTrkHits = fastTrackerSettings.minSiliconHits;
} else if (enableSecondarySmearing) {
o2::track::TrackParCov perfectTrackParCov;
o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectTrackParCov, pdgDB);
computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov);
perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode()));
nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta);
if (nTrkHits < fastTrackerSettings.minSiliconHits) {
Expand Down
Loading