3030#include < TLorentzVector.h>
3131#include < TRandom3.h>
3232
33- #include < array>
3433#include < cmath>
3534#include < cstddef>
3635#include < vector>
@@ -47,46 +46,49 @@ class Decayer
4746 Decayer () = default ;
4847
4948 template <typename TDatabase>
50- std::vector<o2::upgrade::OTFParticle> decayParticle (const TDatabase& pdgDB, const o2::track::TrackParCov& track, const int pdgCode )
49+ std::vector<o2::upgrade::OTFParticle> decayParticle (const TDatabase& pdgDB, const OTFParticle& particle )
5150 {
52- const auto & particleInfo = pdgDB->GetParticle (pdgCode);
51+ const auto & particleInfo = pdgDB->GetParticle (particle.pdgCode ());
52+ if (!particleInfo) {
53+ return {};
54+ }
55+
5356 const int charge = particleInfo->Charge () / 3 ;
5457 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);
5958
60- const double u = mRand3 .Uniform (0 , 1 );
59+ const double u = mRand3 .Uniform (0.001 , 0.999 );
6160 const double ctau = o2::constants::physics::LightSpeedCm2S * particleInfo->Lifetime (); // cm
62- const double betaGamma = track. getP () / mass;
61+ const double betaGamma = particle. p () / mass;
6362 const double rxyz = -betaGamma * ctau * std::log (1 - u);
6463 double vx, vy, vz;
6564 double px, py, e;
6665
6766 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[ 1 ] ;
67+ vx = particle. vx () + rxyz * (particle. px () / particle. p ());
68+ vy = particle. vy () + rxyz * (particle. py () / particle. p ());
69+ vz = particle. vz () + rxyz * (particle. pz () / particle. p ());
70+ px = particle. px () ;
71+ py = particle. py () ;
7372 } else {
74- float sna, csa ;
73+ o2::track::TrackParCov track ;
7574 o2::math_utils::CircleXYf_t circle;
75+ o2::upgrade::convertOTFParticleToO2Track (particle, track, pdgDB);
76+
77+ float sna, csa;
7678 track.getCircleParams (mBz , circle, sna, csa);
7779 const double rxy = rxyz / std::sqrt (1 . + track.getTgl () * track.getTgl ());
7880 const double theta = rxy / circle.rC ;
7981
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 ());
82+ vx = ((particle. vx () - circle.xC ) * std::cos (theta) - (particle. vy () - circle.yC ) * std::sin (theta)) + circle.xC ;
83+ vy = ((particle. vy () - circle.yC ) * std::cos (theta) + (particle. vx () - circle.xC ) * std::sin (theta)) + circle.yC ;
84+ vz = particle. vz () + rxyz * (particle. pz () / track.getP ());
8385
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+ px = particle. px () * std::cos (theta) - particle. py () * std::sin (theta);
87+ py = particle. py () * std::cos (theta) + particle. px () * std::sin (theta);
8688 }
8789
8890 double brTotal = 0 .;
89- e = std::sqrt (mass * mass + px * px + py * py + mom[ 2 ] * mom[ 2 ] );
91+ e = std::sqrt (mass * mass + px * px + py * py + particle. pz () * particle. pz () );
9092 for (int ch = 0 ; ch < particleInfo->NDecayChannels (); ++ch) {
9193 brTotal += particleInfo->DecayChannel (ch)->BranchingRatio ();
9294 }
@@ -108,7 +110,11 @@ class Decayer
108110 }
109111 }
110112
111- TLorentzVector tlv (px, py, mom[2 ], e);
113+ if (dauMasses.empty ()) {
114+ return {};
115+ }
116+
117+ TLorentzVector tlv (px, py, particle.pz (), e);
112118 TGenPhaseSpace decay;
113119 decay.SetDecay (tlv, dauMasses.size (), dauMasses.data ());
114120 decay.Generate ();
0 commit comments