Skip to content

[ALICE3] Add smearing to shortlived particles in LUT process function in otf tracker #17737

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
njacazio merged 1 commit into from
Sep 2, 2026
+74 −46
Merged

[ALICE3] Add smearing to shortlived particles in LUT process function in otf tracker #17737

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
90 changes: 54 additions & 36 deletions ALICE3/Core/Decayer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,16 +28,15 @@

#include <TDecayChannel.h> // IWYU pragma: keep
#include <TGenPhaseSpace.h>
#include <TLorentzVector.h>

Check failure on line 31 in ALICE3/Core/Decayer.h

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 <TRandom3.h>

#include <array>
#include <cmath>
#include <cstddef>
#include <vector>

namespace o2
{
namespace upgrade
namespace o2::upgrade
{

class Decayer
Expand All @@ -47,44 +46,27 @@
Decayer() = default;

template <typename TDatabase>
std::vector<o2::upgrade::OTFParticle> decayParticle(const TDatabase& pdgDB, const OTFParticle& particle)
std::vector<o2::upgrade::OTFParticle> 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<double, 3> decayVtx = generateDecayVertex<double>(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.;
Expand Down Expand Up @@ -114,7 +96,7 @@
return {};
}

TLorentzVector tlv(px, py, particle.pz(), e);

Check failure on line 99 in ALICE3/Core/Decayer.h

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 decay;
decay.SetDecay(tlv, dauMasses.size(), dauMasses.data());
decay.Generate();
Expand All @@ -122,7 +104,7 @@
std::vector<o2::upgrade::OTFParticle> decayProducts;
for (size_t i = 0; i < dauMasses.size(); ++i) {
o2::upgrade::OTFParticle particle;
TLorentzVector dau = *decay.GetDecay(i);

Check failure on line 107 in ALICE3/Core/Decayer.h

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.
particle.setPDG(pdgCodesDaughters[i]);
particle.setVxVyVz(mVx, mVy, mVz);
particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E());
Expand All @@ -133,6 +115,42 @@
return decayProducts;
}

template <typename T = float, typename TDatabase, typename TParticle>
std::array<T, 3> generateDecayVertex(const TParticle& particle, const TDatabase& pdgDB)
{
std::array<T, 3> 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)
Expand All @@ -142,18 +160,18 @@
}

// Getters
float getSecondaryVertexX() const { return static_cast<float>(mVx); }
float getSecondaryVertexY() const { return static_cast<float>(mVy); }
float getSecondaryVertexZ() const { return static_cast<float>(mVz); }
float getDecayRadius() const { return static_cast<float>(std::hypot(mVx, mVy)); }
[[nodiscard]] float getSecondaryVertexX() const { return static_cast<float>(mVx); }
[[nodiscard]] float getSecondaryVertexY() const { return static_cast<float>(mVy); }
[[nodiscard]] float getSecondaryVertexZ() const { return static_cast<float>(mVz); }
[[nodiscard]] float getDecayRadius() const { return static_cast<float>(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_
2 changes: 1 addition & 1 deletion ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@
}

particle.setBitOff(o2::upgrade::DecayerBits::IsAlive);
std::vector<o2::upgrade::OTFParticle> decayStack = decayer.decayParticle(pdgDB, particle);
std::vector<o2::upgrade::OTFParticle> decayStack = decayer.decayParticle(particle, pdgDB);
if (decayStack.empty()) {
continue;
}
Expand All @@ -175,7 +175,7 @@
particle.setDecayRadius(std::hypot(decayer.getSecondaryVertexX(), decayer.getSecondaryVertexY()));
const float trackTimeNS = trackLength / trackVelocity * PicoToNano;
particle.setIndicesDaughter(particlesInDataframe - indexOffset + allParticles.size(), particlesInDataframe - indexOffset + allParticles.size() + (decayStack.size() - 1));
for (auto& daughter : decayStack) {

Check failure on line 178 in ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
daughter.setIndicesMother(particlesInDataframe - indexOffset + i, particlesInDataframe - indexOffset + i);
daughter.setCollisionId(particle.collisionId());
daughter.setBitOn(o2::upgrade::DecayerBits::ProducedByDecayer);
Expand Down Expand Up @@ -213,7 +213,7 @@
decayParticles(0, allParticles.size());

// Fill output table
for (auto& otfParticle : allParticles) {

Check failure on line 216 in ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
otfParticle.setIndexOffset(indexOffset);
if (otfParticle.hasNaN()) {
histos.fill(HIST("hNaNBookkeeping"), 1);
Expand Down
28 changes: 19 additions & 9 deletions ALICE3/TableProducer/OTF/onTheFlyTracker.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -77,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 @@ -378,6 +379,9 @@

// Track smearer array, one per geometry
std::vector<std::unique_ptr<o2::delphes::TrackSmearer>> mSmearer;
// Configuration defined at init time
o2::fastsim::GeometryContainer mGeoContainer;
float mMagneticField = 0.0f;

// For processing and vertexing
std::vector<TrackAlice3> recoPrimaries;
Expand All @@ -395,10 +399,8 @@
// For TGenPhaseSpace seed
TRandom3 rand;
Service<o2::ccdb::BasicCCDBManager> 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;
Expand Down Expand Up @@ -438,6 +440,7 @@

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<o2::delphes::TrackSmearer>());
Expand Down Expand Up @@ -809,7 +812,7 @@
/// \param xiDecayVertex the address of the xi decay vertex
/// \param laDecayVertex the address of the la decay vertex
template <typename McParticleType>
void decayCascade(const McParticleType& particle, o2::track::TrackParCov track, std::vector<TLorentzVector>& decayDaughters, std::vector<double>& xiDecayVertex, std::vector<double>& laDecayVertex)

Check failure on line 815 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 @@ -832,12 +835,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 838 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 843 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 @@ -860,7 +863,7 @@
/// \param decayDaughters the address of resulting daughters
/// \param v0DecayVertex the address of the la decay vertex
template <typename McParticleType>
void decayV0Particle(const McParticleType& particle, std::vector<TLorentzVector>& decayDaughters, std::vector<double>& v0DecayVertex, int pdgCode)

Check failure on line 866 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 @@ -1910,7 +1913,6 @@
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;
}
Expand Down Expand Up @@ -1950,15 +1952,24 @@
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<float, 3> 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);
Expand Down Expand Up @@ -2155,8 +2166,7 @@
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());
Expand Down
Loading