From a65e0637228bc4f7fd3240d8a38ad27bcfa822db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alexander=20Tiek=C3=B6tter?= Date: Fri, 22 May 2026 10:19:31 +0200 Subject: [PATCH 1/2] Implement bremsstrahlung in OTF tracker --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 59 ++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index c7da0654e42..d5e72b336df 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -246,6 +246,14 @@ struct OnTheFlyTracker { Configurable 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 radiateBR{"radiateBR", false, "simulate bremsstrahlung"}; + Configurable minBREnergyFraction{"minEnergyFraction", 0.001f, "Minimum energy fraction a bremsstrahlung photon can carry"}; + Configurable maxBREnergyFraction{"maxEnergyFraction", 0.95f, "Maximum energy fraction a bremsstrahlung photon can carry"}; + Configurable radiationStrength{"radiationStrength", 1e-6f, "Strenght of the bremsstrahlung radiation"}; + } brSettings; + using PVertex = o2::dataformats::PrimaryVertex; // for secondary vertex finding @@ -1749,6 +1757,53 @@ struct OnTheFlyTracker { } } + /// 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); + + 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); + } + } + } + } + void processWithLUTs(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const int icfg) { const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; @@ -1811,12 +1866,14 @@ struct OnTheFlyTracker { 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; } @@ -2002,11 +2059,13 @@ struct OnTheFlyTracker { 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) { From 3616a2d7cf0bac71237435dda2eb0fe3c4438710 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alexander=20Tiek=C3=B6tter?= Date: Wed, 27 May 2026 12:53:28 +0200 Subject: [PATCH 2/2] Add basic QA histograms for BR --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 26 ++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index d5e72b336df..21b6b8aa9fa 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -193,6 +193,9 @@ struct OnTheFlyTracker { 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 @@ -249,6 +252,7 @@ struct OnTheFlyTracker { struct : ConfigurableGroup { std::string prefix = "cfgBR"; // configuration for bremsstrahlung Configurable radiateBR{"radiateBR", false, "simulate bremsstrahlung"}; + Configurable doBRQA{"doBRQA", false, "Do QA for bremsstrahlung"}; Configurable minBREnergyFraction{"minEnergyFraction", 0.001f, "Minimum energy fraction a bremsstrahlung photon can carry"}; Configurable maxBREnergyFraction{"maxEnergyFraction", 0.95f, "Maximum energy fraction a bremsstrahlung photon can carry"}; Configurable radiationStrength{"radiationStrength", 1e-6f, "Strenght of the bremsstrahlung radiation"}; @@ -576,6 +580,14 @@ struct OnTheFlyTracker { 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 @@ -1796,10 +1808,24 @@ struct OnTheFlyTracker { 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(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()); + } } } }