From 483988d3b729a7190ee098f7d72e26718c4f22b4 Mon Sep 17 00:00:00 2001 From: jesgum Date: Wed, 2 Sep 2026 18:36:18 +0200 Subject: [PATCH] Add smearing to shortlived particles in LUT process function in otf tracker --- ALICE3/Core/Decayer.h | 90 ++++++++++++-------- ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx | 2 +- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 28 ++++-- 3 files changed, 74 insertions(+), 46 deletions(-) diff --git a/ALICE3/Core/Decayer.h b/ALICE3/Core/Decayer.h index 52ca0a75adb..8f88e4e311b 100644 --- a/ALICE3/Core/Decayer.h +++ b/ALICE3/Core/Decayer.h @@ -31,13 +31,12 @@ #include #include +#include #include #include #include -namespace o2 -{ -namespace upgrade +namespace o2::upgrade { class Decayer @@ -47,44 +46,27 @@ class Decayer Decayer() = default; template - std::vector decayParticle(const TDatabase& pdgDB, const OTFParticle& particle) + std::vector decayParticle(const OTFParticle& particle, const TDatabase& pdgDB) { - const auto& particleInfo = pdgDB->GetParticle(particle.pdgCode()); + auto particleInfo = pdgDB->GetParticle(particle.pdgCode()); if (!particleInfo) { return {}; } const int charge = particleInfo->Charge() / 3; const double mass = particleInfo->Mass(); - - const double u = mRand3.Uniform(0.001, 0.999); - const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm - const double betaGamma = particle.p() / mass; - const double rxyz = -betaGamma * ctau * std::log(1 - u); - double px, py, e; + std::array decayVtx = generateDecayVertex(particle, pdgDB); + mVx = decayVtx[0]; + mVy = decayVtx[1]; + mVz = decayVtx[2]; + double px{}, py{}, e{}; if (!charge) { - mVx = particle.vx() + rxyz * (particle.px() / particle.p()); - mVy = particle.vy() + rxyz * (particle.py() / particle.p()); - mVz = particle.vz() + rxyz * (particle.pz() / particle.p()); px = particle.px(); py = particle.py(); } else { - o2::track::TrackParCov track; - o2::math_utils::CircleXYf_t circle; - o2::upgrade::convertOTFParticleToO2Track(particle, track, pdgDB); - - float sna{}, csa{}; - track.getCircleParams(mBz, circle, sna, csa); - const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl()); - const double theta = rxy / circle.rC; - - mVx = ((particle.vx() - circle.xC) * std::cos(theta) - (particle.vy() - circle.yC) * std::sin(theta)) + circle.xC; - mVy = ((particle.vy() - circle.yC) * std::cos(theta) + (particle.vx() - circle.xC) * std::sin(theta)) + circle.yC; - mVz = particle.vz() + rxyz * (particle.pz() / track.getP()); - - px = particle.px() * std::cos(theta) - particle.py() * std::sin(theta); - py = particle.py() * std::cos(theta) + particle.px() * std::sin(theta); + px = particle.px() * std::cos(mTheta) - particle.py() * std::sin(mTheta); + py = particle.py() * std::cos(mTheta) + particle.px() * std::sin(mTheta); } double brTotal = 0.; @@ -133,6 +115,42 @@ class Decayer return decayProducts; } + template + std::array generateDecayVertex(const TParticle& particle, const TDatabase& pdgDB) + { + std::array decayVertex{}; + auto particleInfo = pdgDB->GetParticle(particle.pdgCode()); + if (!particleInfo) { + return {}; + } + + const int charge = particleInfo->Charge() / 3; + const double mass = particleInfo->Mass(); + const double u = mRand3.Uniform(0.001, 0.999); + const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm + const double betaGamma = particle.p() / mass; + const double rxyz = -betaGamma * ctau * std::log(1 - u); + + if (!charge) { + decayVertex[0] = particle.vx() + rxyz * (particle.px() / particle.p()); + decayVertex[1] = particle.vy() + rxyz * (particle.py() / particle.p()); + decayVertex[2] = particle.vz() + rxyz * (particle.pz() / particle.p()); + } else { + o2::math_utils::CircleXYf_t circle; + o2::track::TrackParCov track = o2::upgrade::convertMCParticleToO2Track(particle, pdgDB); + + float sna{}, csa{}; + track.getCircleParams(mBz, circle, sna, csa); + const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl()); + mTheta = rxy / circle.rC; + + decayVertex[0] = ((particle.vx() - circle.xC) * std::cos(mTheta) - (particle.vy() - circle.yC) * std::sin(mTheta)) + circle.xC; + decayVertex[1] = ((particle.vy() - circle.yC) * std::cos(mTheta) + (particle.vx() - circle.xC) * std::sin(mTheta)) + circle.yC; + decayVertex[2] = particle.vz() + rxyz * (particle.pz() / track.getP()); + } + return decayVertex; + } + // Setters void setBField(const double b) { mBz = b; } void setSeed(const int seed) @@ -142,18 +160,18 @@ class Decayer } // Getters - float getSecondaryVertexX() const { return static_cast(mVx); } - float getSecondaryVertexY() const { return static_cast(mVy); } - float getSecondaryVertexZ() const { return static_cast(mVz); } - float getDecayRadius() const { return static_cast(std::hypot(mVx, mVy)); } + [[nodiscard]] float getSecondaryVertexX() const { return static_cast(mVx); } + [[nodiscard]] float getSecondaryVertexY() const { return static_cast(mVy); } + [[nodiscard]] float getSecondaryVertexZ() const { return static_cast(mVz); } + [[nodiscard]] float getDecayRadius() const { return static_cast(std::hypot(mVx, mVy)); } private: double mBz{20.}; // kG double mVx{-1.}, mVy{-1.}, mVz{-1.}; - TRandom3 mRand3{}; + double mTheta{}; + TRandom3 mRand3; }; -} // namespace upgrade -} // namespace o2 +} // namespace o2::upgrade #endif // ALICE3_CORE_DECAYER_H_ diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx index c03447ecb75..2c8c4dea0c9 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx @@ -152,7 +152,7 @@ struct OnTheFlyDecayer { } particle.setBitOff(o2::upgrade::DecayerBits::IsAlive); - std::vector decayStack = decayer.decayParticle(pdgDB, particle); + std::vector decayStack = decayer.decayParticle(particle, pdgDB); if (decayStack.empty()) { continue; } diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index b2d55f13015..949d707928d 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -23,6 +23,7 @@ /// \author Roberto Preghenella preghenella@bo.infn.it /// +#include "ALICE3/Core/Decayer.h" #include "ALICE3/Core/DetLayer.h" #include "ALICE3/Core/FastTracker.h" #include "ALICE3/Core/FlatTrackSmearer.h" @@ -378,6 +379,9 @@ struct OnTheFlyTracker { // Track smearer array, one per geometry std::vector> mSmearer; + // Configuration defined at init time + o2::fastsim::GeometryContainer mGeoContainer; + float mMagneticField = 0.0f; // For processing and vertexing std::vector recoPrimaries; @@ -395,10 +399,8 @@ struct OnTheFlyTracker { // For TGenPhaseSpace seed TRandom3 rand; Service ccdb{}; + o2::upgrade::Decayer decayer; - // Configuration defined at init time - o2::fastsim::GeometryContainer mGeoContainer; - float mMagneticField = 0.0f; // Time resolution constants static constexpr float timeResolutionNs = 100.f; // ns static constexpr float nsToMus = 1e-3f; @@ -438,6 +440,7 @@ struct OnTheFlyTracker { const int nGeometries = mGeoContainer.getNumberOfConfigurations(); mMagneticField = mGeoContainer.getFloatValue(0, "global", "magneticfield"); + decayer.setBField(mMagneticField); for (int icfg = 0; icfg < nGeometries; ++icfg) { const std::string histPath = "Configuration_" + std::to_string(icfg) + "/"; mSmearer.emplace_back(std::make_unique()); @@ -1910,7 +1913,6 @@ struct OnTheFlyTracker { uint32_t multiplicityCounter = 0; // Now that the multiplicity is known, we can process the particles to smear them for (const auto& mcParticle : mcParticles) { - if (!mcParticle.isPhysicalPrimary()) { continue; } @@ -1950,15 +1952,24 @@ struct OnTheFlyTracker { bool reconstructed = true; int nTrkHits = 0; if (enablePrimarySmearing) { - if (fastPrimaryTrackerSettings.fastTrackPrimaries || fastPrimaryTrackerSettings.fastTrackShortLivedParticles) { - o2::track::TrackParCov perfectTrackParCov; - o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectTrackParCov, pdgDB); + if (fastPrimaryTrackerSettings.fastTrackPrimaries && longLivedToBeHandled) { + o2::track::TrackParCov perfectTrackParCov = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdgDB); perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode())); computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov); nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta); if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) { reconstructed = false; } + } else if (fastPrimaryTrackerSettings.fastTrackShortLivedParticles && shortLivedToBeHandled) { + o2::track::TrackParCov perfectTrackParCov = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdgDB); + perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode())); + computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov); + const std::array decayVtx = decayer.generateDecayVertex(mcParticle, pdgDB); + const float decayRadius2D = std::hypot(decayVtx[0], decayVtx[1]); + nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta, decayRadius2D); + if (nTrkHits < fastPrimaryTrackerSettings.minSiliconHits) { + reconstructed = false; + } } else { o2::upgrade::convertMCParticleToO2Track(mcParticle, trackParCov, pdgDB); computeBremsstrahlungLoss(icfg, mcParticle, trackParCov); @@ -2155,8 +2166,7 @@ struct OnTheFlyTracker { computeBremsstrahlungLoss(icfg, mcParticle, trackParCov); reconstructed = mSmearer[icfg]->smearTrack(trackParCov, mcParticle.pdgCode(), dNdEta); } else if (shortLivedToBeHandled && fastPrimaryTrackerSettings.fastTrackShortLivedParticles) { - o2::track::TrackParCov perfectTrackParCov; - o2::upgrade::convertMCParticleToO2Track(mcParticle, perfectTrackParCov, pdgDB); + o2::track::TrackParCov perfectTrackParCov = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdgDB); perfectTrackParCov.setPID(pdgCodeToPID(mcParticle.pdgCode())); computeBremsstrahlungLoss(icfg, mcParticle, perfectTrackParCov); nTrkHits = fastTracker[icfg]->FastTrack(perfectTrackParCov, trackParCov, dNdEta, mcParticle.decayRadius());