Skip to content

Commit e5aefca

Browse files
jesgumjesgum
authored
[ALICE3] First iteration of OTF decayer (#14441)
1 parent 81270aa commit e5aefca

File tree

7 files changed

+633
-0
lines changed

7 files changed

+633
-0
lines changed

‎ALICE3/Core/Decayer.h‎

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
///
13+
/// \file Decayer.h
14+
/// \author Jesper Karlsson Gumprecht
15+
/// \since 15/12/2025
16+
/// \brief Basic class to handle short-lived particle decays in the fast simulation
17+
///
18+
19+
#ifndef ALICE3_CORE_DECAYER_H_
20+
#define ALICE3_CORE_DECAYER_H_
21+
22+
#include "ALICE3/Core/TrackUtilities.h"
23+
24+
#include "ReconstructionDataFormats/Track.h"
25+
26+
#include <TDatabasePDG.h>
27+
#include <TDecayChannel.h>
28+
#include <TGenPhaseSpace.h>
29+
#include <TLorentzVector.h>
30+
#include <TParticlePDG.h>
31+
#include <TRandom3.h>
32+
33+
#include <array>
34+
#include <cmath>
35+
#include <string>
36+
#include <vector>
37+
38+
namespace o2
39+
{
40+
namespace upgrade
41+
{
42+
43+
class Decayer
44+
{
45+
public:
46+
// Default constructor
47+
Decayer() = default;
48+
49+
template <typename TDatabase>
50+
std::vector<o2::upgrade::OTFParticle> decayParticle(const TDatabase& pdgDB, const o2::track::TrackParCov& track, const int pdgCode)
51+
{
52+
const auto& particleInfo = pdgDB->GetParticle(pdgCode);
53+
const int charge = particleInfo->Charge() / 3;
54+
const double mass = particleInfo->Mass();
55+
std::array<float, 3> mom;
56+
std::array<float, 3> pos;
57+
track.getPxPyPzGlo(mom);
58+
track.getXYZGlo(pos);
59+
60+
const double u = mRand3.Uniform(0, 1);
61+
const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime(); // cm
62+
const double betaGamma = track.getP() / mass;
63+
const double rxyz = -betaGamma * ctau * std::log(1 - u);
64+
double vx, vy, vz;
65+
double px, py, e;
66+
67+
if (!charge) {
68+
vx = pos[0] + rxyz * (mom[0] / track.getP());
69+
vy = pos[1] + rxyz * (mom[1] / track.getP());
70+
vz = pos[2] + rxyz * (mom[2] / track.getP());
71+
px = mom[0];
72+
py = mom[2];
73+
} else {
74+
float sna, csa;
75+
o2::math_utils::CircleXYf_t circle;
76+
track.getCircleParams(mBz, circle, sna, csa);
77+
const double rxy = rxyz / std::sqrt(1. + track.getTgl() * track.getTgl());
78+
const double theta = rxy / circle.rC;
79+
80+
vx = ((pos[0] - circle.xC) * std::cos(theta) - (pos[1] - circle.yC) * std::sin(theta)) + circle.xC;
81+
vy = ((pos[1] - circle.yC) * std::cos(theta) + (pos[0] - circle.xC) * std::sin(theta)) + circle.yC;
82+
vz = mom[2] + rxyz * (mom[2] / track.getP());
83+
84+
px = mom[0] * std::cos(theta) - mom[1] * std::sin(theta);
85+
py = mom[1] * std::cos(theta) + mom[0] * std::sin(theta);
86+
}
87+
88+
double brTotal = 0.;
89+
e = std::sqrt(mass * mass + px * px + py * py + mom[2] * mom[2]);
90+
for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) {
91+
brTotal += particleInfo->DecayChannel(ch)->BranchingRatio();
92+
}
93+
94+
double brSum = 0.;
95+
std::vector<double> dauMasses;
96+
std::vector<int> pdgCodesDaughters;
97+
const double randomChannel = mRand3.Uniform(0., brTotal);
98+
for (int ch = 0; ch < particleInfo->NDecayChannels(); ++ch) {
99+
brSum += particleInfo->DecayChannel(ch)->BranchingRatio();
100+
if (randomChannel < brSum) {
101+
for (int dau = 0; dau < particleInfo->DecayChannel(ch)->NDaughters(); ++dau) {
102+
const int pdgDau = particleInfo->DecayChannel(ch)->DaughterPdgCode(dau);
103+
pdgCodesDaughters.push_back(pdgDau);
104+
const auto& dauInfo = pdgDB->GetParticle(pdgDau);
105+
dauMasses.push_back(dauInfo->Mass());
106+
}
107+
break;
108+
}
109+
}
110+
111+
TLorentzVector tlv(px, py, mom[2], e);
112+
TGenPhaseSpace decay;
113+
decay.SetDecay(tlv, dauMasses.size(), dauMasses.data());
114+
decay.Generate();
115+
116+
std::vector<o2::upgrade::OTFParticle> decayProducts;
117+
for (size_t i = 0; i < dauMasses.size(); ++i) {
118+
o2::upgrade::OTFParticle particle;
119+
TLorentzVector dau = *decay.GetDecay(i);
120+
particle.setPDG(pdgCodesDaughters[i]);
121+
particle.setVxVyVz(vx, vy, vz);
122+
particle.setPxPyPzE(dau.Px(), dau.Py(), dau.Pz(), dau.E());
123+
decayProducts.push_back(particle);
124+
}
125+
126+
return decayProducts;
127+
}
128+
129+
void setSeed(const int seed)
130+
{
131+
mRand3.SetSeed(seed); // For decay length sampling
132+
gRandom->SetSeed(seed); // For TGenPhaseSpace
133+
}
134+
135+
void setBField(const double b) { mBz = b; }
136+
137+
private:
138+
TRandom3 mRand3;
139+
double mBz;
140+
};
141+
142+
} // namespace upgrade
143+
} // namespace o2
144+
145+
#endif // ALICE3_CORE_DECAYER_H_

‎ALICE3/Core/TrackUtilities.h‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,40 @@
2727
namespace o2::upgrade
2828
{
2929

30+
/// Struct to store mc info for the otf decayer
31+
struct OTFParticle {
32+
int mPdgCode;
33+
float mE;
34+
float mVx, mVy, mVz;
35+
float mPx, mPy, mPz;
36+
37+
// Setters
38+
void setPDG(int pdg) { mPdgCode = pdg; }
39+
void setVxVyVz(float vx, float vy, float vz)
40+
{
41+
mVx = vx;
42+
mVy = vy;
43+
mVz = vz;
44+
}
45+
void setPxPyPzE(float px, float py, float pz, float e)
46+
{
47+
mPx = px;
48+
mPy = py;
49+
mPz = pz;
50+
mE = e;
51+
}
52+
53+
// Getters
54+
int pdgCode() const { return mPdgCode; }
55+
float vx() const { return mVx; }
56+
float vy() const { return mVy; }
57+
float vz() const { return mVz; }
58+
float px() const { return mPx; }
59+
float py() const { return mPy; }
60+
float pz() const { return mPz; }
61+
float e() const { return mE; }
62+
};
63+
3064
/// Function to convert a TLorentzVector into a perfect Track
3165
/// \param charge particle charge (integer)
3266
/// \param particle the particle to convert (TLorentzVector)
@@ -58,6 +92,18 @@ void convertTLorentzVectorToO2Track(int pdgCode,
5892
convertTLorentzVectorToO2Track(charge, particle, productionVertex, o2track);
5993
}
6094

95+
/// Function to convert a OTFParticle into a perfect Track
96+
/// \param particle the particle to convert (OTFParticle)
97+
/// \param o2track the address of the resulting TrackParCov
98+
/// \param pdg the pdg service
99+
template <typename PdgService>
100+
void convertOTFParticleToO2Track(const OTFParticle& particle, o2::track::TrackParCov& o2track, const PdgService& pdg)
101+
{
102+
static TLorentzVector tlv;
103+
tlv.SetPxPyPzE(particle.px(), particle.py(), particle.pz(), particle.e());
104+
convertTLorentzVectorToO2Track(particle.pdgCode(), tlv, {particle.vx(), particle.vy(), particle.vz()}, o2track, pdg);
105+
}
106+
61107
/// Function to convert a McParticle into a perfect Track
62108
/// \param particle the particle to convert (mcParticle)
63109
/// \param o2track the address of the resulting TrackParCov

‎ALICE3/DataModel/OTFMCParticle.h‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
///
13+
/// \file OTFMCParticle.h
14+
/// \author Jesper Karlsson Gumprecht
15+
/// \since 16/12/2025
16+
/// \brief Redefinition of the mcparticles table specifically for the fast sim
17+
///
18+
19+
#ifndef ALICE3_DATAMODEL_OTFMCPARTICLE_H_
20+
#define ALICE3_DATAMODEL_OTFMCPARTICLE_H_
21+
22+
// O2 includes
23+
#include "Framework/AnalysisDataModel.h"
24+
25+
namespace o2::aod
26+
{
27+
28+
namespace otfmcparticle
29+
{
30+
DECLARE_SOA_COLUMN(Phi, phi, float);
31+
DECLARE_SOA_COLUMN(Eta, eta, float);
32+
DECLARE_SOA_COLUMN(Pt, pt, float);
33+
DECLARE_SOA_COLUMN(P, p, float);
34+
DECLARE_SOA_COLUMN(Y, y, float);
35+
} // namespace otfmcparticle
36+
37+
DECLARE_SOA_TABLE_FULL(McParticlesWithDau, "McParticlesWithDau", "AOD", "MCPARTICLEWITHDAU",
38+
o2::soa::Index<>,
39+
mcparticle::McCollisionId,
40+
mcparticle::PdgCode,
41+
mcparticle::StatusCode,
42+
mcparticle::Flags,
43+
mcparticle::MothersIds,
44+
mcparticle::DaughtersIdSlice,
45+
mcparticle::Weight,
46+
mcparticle::Px,
47+
mcparticle::Py,
48+
mcparticle::Pz,
49+
mcparticle::E,
50+
mcparticle::Vx,
51+
mcparticle::Vy,
52+
mcparticle::Vz,
53+
mcparticle::Vt,
54+
otfmcparticle::Phi,
55+
otfmcparticle::Eta,
56+
otfmcparticle::Pt,
57+
otfmcparticle::P,
58+
otfmcparticle::Y,
59+
mcparticle::PVector<mcparticle::Px, mcparticle::Py, mcparticle::Pz>,
60+
mcparticle::ProducedByGenerator<mcparticle::Flags>,
61+
mcparticle::FromBackgroundEvent<mcparticle::Flags>,
62+
mcparticle::GetGenStatusCode<mcparticle::Flags, mcparticle::StatusCode>,
63+
mcparticle::GetHepMCStatusCode<mcparticle::Flags, mcparticle::StatusCode>,
64+
mcparticle::GetProcess<mcparticle::Flags, mcparticle::StatusCode>,
65+
mcparticle::IsPhysicalPrimary<mcparticle::Flags>);
66+
67+
using McParticleWithDau = McParticlesWithDau::iterator;
68+
69+
} // namespace o2::aod
70+
71+
#endif // ALICE3_DATAMODEL_OTFMCPARTICLE_H_

‎ALICE3/TableProducer/OTF/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,11 @@
99
# granted to it by virtue of its status as an Intergovernmental Organization
1010
# or submit itself to any jurisdiction.
1111

12+
o2physics_add_dpl_workflow(onthefly-decayer
13+
SOURCES onTheFlyDecayer.cxx
14+
PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2::DCAFitter O2Physics::ALICE3Core
15+
COMPONENT_NAME Analysis)
16+
1217
o2physics_add_dpl_workflow(onthefly-tracker
1318
SOURCES onTheFlyTracker.cxx
1419
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2::DCAFitter O2Physics::ALICE3Core O2Physics::FastTracker

0 commit comments

Comments
 (0)