3131#include < TLorentzVector.h>
3232#include < TRandom3.h>
3333
34+ #include < array>
3435#include < cmath>
3536#include < cstddef>
3637#include < vector>
3738
38- namespace o2
39- {
40- namespace upgrade
39+ namespace o2 ::upgrade
4140{
4241
4342class Decayer
@@ -47,44 +46,27 @@ class Decayer
4746 Decayer () = default ;
4847
4948 template <typename TDatabase>
50- std::vector<o2::upgrade::OTFParticle> decayParticle (const TDatabase& pdgDB , const OTFParticle& particle )
49+ std::vector<o2::upgrade::OTFParticle> decayParticle (const OTFParticle& particle , const TDatabase& pdgDB )
5150 {
52- const auto & particleInfo = pdgDB->GetParticle (particle.pdgCode ());
51+ auto particleInfo = pdgDB->GetParticle (particle.pdgCode ());
5352 if (!particleInfo) {
5453 return {};
5554 }
5655
5756 const int charge = particleInfo->Charge () / 3 ;
5857 const double mass = particleInfo->Mass ();
59-
60- const double u = mRand3 .Uniform (0.001 , 0.999 );
61- const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime (); // cm
62- const double betaGamma = particle.p () / mass;
63- const double rxyz = -betaGamma * ctau * std::log (1 - u);
64- double px, py, e;
58+ std::array<double , 3 > decayVtx = generateDecayVertex<double >(particle, pdgDB);
59+ mVx = decayVtx[0 ];
60+ mVy = decayVtx[1 ];
61+ mVz = decayVtx[2 ];
62+ double px{}, py{}, e{};
6563
6664 if (!charge) {
67- mVx = particle.vx () + rxyz * (particle.px () / particle.p ());
68- mVy = particle.vy () + rxyz * (particle.py () / particle.p ());
69- mVz = particle.vz () + rxyz * (particle.pz () / particle.p ());
7065 px = particle.px ();
7166 py = particle.py ();
7267 } else {
73- o2::track::TrackParCov track;
74- o2::math_utils::CircleXYf_t circle;
75- o2::upgrade::convertOTFParticleToO2Track (particle, track, pdgDB);
76-
77- float sna{}, csa{};
78- track.getCircleParams (mBz , circle, sna, csa);
79- const double rxy = rxyz / std::sqrt (1 . + track.getTgl () * track.getTgl ());
80- const double theta = rxy / circle.rC ;
81-
82- mVx = ((particle.vx () - circle.xC ) * std::cos (theta) - (particle.vy () - circle.yC ) * std::sin (theta)) + circle.xC ;
83- mVy = ((particle.vy () - circle.yC ) * std::cos (theta) + (particle.vx () - circle.xC ) * std::sin (theta)) + circle.yC ;
84- mVz = particle.vz () + rxyz * (particle.pz () / track.getP ());
85-
86- px = particle.px () * std::cos (theta) - particle.py () * std::sin (theta);
87- py = particle.py () * std::cos (theta) + particle.px () * std::sin (theta);
68+ px = particle.px () * std::cos (mTheta ) - particle.py () * std::sin (mTheta );
69+ py = particle.py () * std::cos (mTheta ) + particle.px () * std::sin (mTheta );
8870 }
8971
9072 double brTotal = 0 .;
@@ -133,6 +115,42 @@ class Decayer
133115 return decayProducts;
134116 }
135117
118+ template <typename T = float , typename TDatabase, typename TParticle>
119+ std::array<T, 3 > generateDecayVertex (const TParticle& particle, const TDatabase& pdgDB)
120+ {
121+ std::array<T, 3 > decayVertex{};
122+ auto particleInfo = pdgDB->GetParticle (particle.pdgCode ());
123+ if (!particleInfo) {
124+ return {};
125+ }
126+
127+ const int charge = particleInfo->Charge () / 3 ;
128+ const double mass = particleInfo->Mass ();
129+ const double u = mRand3 .Uniform (0.001 , 0.999 );
130+ const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime (); // cm
131+ const double betaGamma = particle.p () / mass;
132+ const double rxyz = -betaGamma * ctau * std::log (1 - u);
133+
134+ if (!charge) {
135+ decayVertex[0 ] = particle.vx () + rxyz * (particle.px () / particle.p ());
136+ decayVertex[1 ] = particle.vy () + rxyz * (particle.py () / particle.p ());
137+ decayVertex[2 ] = particle.vz () + rxyz * (particle.pz () / particle.p ());
138+ } else {
139+ o2::math_utils::CircleXYf_t circle;
140+ o2::track::TrackParCov track = o2::upgrade::convertMCParticleToO2Track (particle, pdgDB);
141+
142+ float sna{}, csa{};
143+ track.getCircleParams (mBz , circle, sna, csa);
144+ const double rxy = rxyz / std::sqrt (1 . + track.getTgl () * track.getTgl ());
145+ mTheta = rxy / circle.rC ;
146+
147+ decayVertex[0 ] = ((particle.vx () - circle.xC ) * std::cos (mTheta ) - (particle.vy () - circle.yC ) * std::sin (mTheta )) + circle.xC ;
148+ decayVertex[1 ] = ((particle.vy () - circle.yC ) * std::cos (mTheta ) + (particle.vx () - circle.xC ) * std::sin (mTheta )) + circle.yC ;
149+ decayVertex[2 ] = particle.vz () + rxyz * (particle.pz () / track.getP ());
150+ }
151+ return decayVertex;
152+ }
153+
136154 // Setters
137155 void setBField (const double b) { mBz = b; }
138156 void setSeed (const int seed)
@@ -142,18 +160,18 @@ class Decayer
142160 }
143161
144162 // Getters
145- float getSecondaryVertexX () const { return static_cast <float >(mVx ); }
146- float getSecondaryVertexY () const { return static_cast <float >(mVy ); }
147- float getSecondaryVertexZ () const { return static_cast <float >(mVz ); }
148- float getDecayRadius () const { return static_cast <float >(std::hypot (mVx , mVy )); }
163+ [[nodiscard]] float getSecondaryVertexX () const { return static_cast <float >(mVx ); }
164+ [[nodiscard]] float getSecondaryVertexY () const { return static_cast <float >(mVy ); }
165+ [[nodiscard]] float getSecondaryVertexZ () const { return static_cast <float >(mVz ); }
166+ [[nodiscard]] float getDecayRadius () const { return static_cast <float >(std::hypot (mVx , mVy )); }
149167
150168 private:
151169 double mBz {20 .}; // kG
152170 double mVx {-1 .}, mVy {-1 .}, mVz {-1 .};
153- TRandom3 mRand3 {};
171+ double mTheta {};
172+ TRandom3 mRand3 ;
154173};
155174
156- } // namespace upgrade
157- } // namespace o2
175+ } // namespace o2::upgrade
158176
159177#endif // ALICE3_CORE_DECAYER_H_
0 commit comments