@@ -112,15 +112,15 @@ struct OnTheFlyTracker {
112112
113113 struct : ConfigurableGroup {
114114 std::string prefix = " lookUpTables" ; // JSON group name
115- Configurable<std::vector<std::string>> lutEl{" lutEl" , std::vector<std::string>{" lutCovm.el.dat" }, " LUT for electrons (if emtpy no LUT is taken)" };
116- Configurable<std::vector<std::string>> lutMu{" lutMu" , std::vector<std::string>{" lutCovm.mu.dat" }, " LUT for muons (if emtpy no LUT is taken)" };
117- Configurable<std::vector<std::string>> lutPi{" lutPi" , std::vector<std::string>{" lutCovm.pi.dat" }, " LUT for pions (if emtpy no LUT is taken)" };
118- Configurable<std::vector<std::string>> lutKa{" lutKa" , std::vector<std::string>{" lutCovm.ka.dat" }, " LUT for kaons (if emtpy no LUT is taken)" };
119- Configurable<std::vector<std::string>> lutPr{" lutPr" , std::vector<std::string>{" lutCovm.pr.dat" }, " LUT for protons (if emtpy no LUT is taken)" };
120- Configurable<std::vector<std::string>> lutDe{" lutDe" , std::vector<std::string>{" " }, " LUT for deuterons (if emtpy no LUT is taken)" };
121- Configurable<std::vector<std::string>> lutTr{" lutTr" , std::vector<std::string>{" " }, " LUT for tritons (if emtpy no LUT is taken)" };
122- Configurable<std::vector<std::string>> lutHe3{" lutHe3" , std::vector<std::string>{" " }, " LUT for Helium-3 (if emtpy no LUT is taken)" };
123- Configurable<std::vector<std::string>> lutAl{" lutAl" , std::vector<std::string>{" " }, " LUT for Alphas (if emtpy no LUT is taken)" };
115+ Configurable<std::vector<std::string>> lutEl{" lutEl" , std::vector<std::string>{" lutCovm.el.dat " }, " LUT for electrons (if emtpy no LUT is taken)" };
116+ Configurable<std::vector<std::string>> lutMu{" lutMu" , std::vector<std::string>{" lutCovm.mu.dat " }, " LUT for muons (if emtpy no LUT is taken)" };
117+ Configurable<std::vector<std::string>> lutPi{" lutPi" , std::vector<std::string>{" lutCovm.pi.dat " }, " LUT for pions (if emtpy no LUT is taken)" };
118+ Configurable<std::vector<std::string>> lutKa{" lutKa" , std::vector<std::string>{" lutCovm.ka.dat " }, " LUT for kaons (if emtpy no LUT is taken)" };
119+ Configurable<std::vector<std::string>> lutPr{" lutPr" , std::vector<std::string>{" lutCovm.pr.dat " }, " LUT for protons (if emtpy no LUT is taken)" };
120+ Configurable<std::vector<std::string>> lutDe{" lutDe" , std::vector<std::string>{" " }, " LUT for deuterons (if emtpy no LUT is taken)" };
121+ Configurable<std::vector<std::string>> lutTr{" lutTr" , std::vector<std::string>{" " }, " LUT for tritons (if emtpy no LUT is taken)" };
122+ Configurable<std::vector<std::string>> lutHe3{" lutHe3" , std::vector<std::string>{" " }, " LUT for Helium-3 (if emtpy no LUT is taken)" };
123+ Configurable<std::vector<std::string>> lutAl{" lutAl" , std::vector<std::string>{" " }, " LUT for Alphas (if emtpy no LUT is taken)" };
124124 } lookUpTables;
125125
126126 struct : ConfigurableGroup {
@@ -308,9 +308,15 @@ struct OnTheFlyTracker {
308308 if (enablePrimarySmearing) {
309309 auto loadLUT = [&](int icfg, int pdg, const std::vector<std::string>& tables) {
310310 const bool foundNewCfg = static_cast <size_t >(icfg) < tables.size ();
311- const std::string& lutFile = foundNewCfg ? tables[icfg] : tables.front ();
311+ std::string lutFile = foundNewCfg ? tables[icfg] : tables.front ();
312+ // strip from leading/trailing spaces
313+ lutFile.erase (0 , lutFile.find_first_not_of (" " ));
314+ lutFile.erase (lutFile.find_last_not_of (" " ) + 1 );
315+ if (lutFile.empty ()) {
316+ LOG (fatal) << " Empty LUT file passed for pdg " << pdg << " , if you don't want to use a LUT remove the entry from the JSON config." ;
317+ }
312318 bool success = mSmearer [icfg]->loadTable (pdg, lutFile.c_str ());
313- if (!success && !lutFile. empty () ) {
319+ if (!success) {
314320 LOG (fatal) << " Having issue with loading the LUT " << pdg << " " << lutFile;
315321 }
316322
@@ -575,52 +581,52 @@ struct OnTheFlyTracker {
575581 // / \param xiDecayVertex the address of the xi decay vertex
576582 // / \param laDecayVertex the address of the la decay vertex
577583 template <typename McParticleType>
578- void decayParticle (McParticleType particle, o2::track::TrackParCov track, std::vector<TLorentzVector>& decayDaughters, std::vector<double >& xiDecayVertex, std::vector<double >& laDecayVertex)
584+ void decayCascade (McParticleType particle, o2::track::TrackParCov track, std::vector<TLorentzVector>& decayDaughters, std::vector<double >& xiDecayVertex, std::vector<double >& laDecayVertex)
579585 {
580- double u = rand.Uniform (0 , 1 );
581- double xi_mass = o2::constants::physics::MassXiMinus;
582- double la_mass = o2::constants::physics::MassLambda;
583- double pi_mass = o2::constants::physics::MassPionCharged;
584- double pr_mass = o2::constants::physics::MassProton;
585-
586- double xi_gamma = 1 / std::sqrt (1 + (particle.p () * particle.p ()) / (xi_mass * xi_mass));
587- double xi_ctau = 4.91 * xi_gamma;
588- double xi_rxyz = (-xi_ctau * std::log (1 - u));
586+ const double uXi = rand.Uniform (0 , 1 );
587+ const double ctauXi = 4.91 ; // cm
588+ const double betaGammaXi = particle.p () / o2::constants::physics::MassXiMinus;
589+ const double rxyzXi = (-betaGammaXi * ctauXi * std::log (1 - uXi));
590+
589591 float sna, csa;
590- o2::math_utils::CircleXYf_t xi_circle;
591- track.getCircleParams (magneticField, xi_circle, sna, csa);
592- double xi_rxy = xi_rxyz / std::sqrt (1 . + track.getTgl () * track.getTgl ());
593- double theta = xi_rxy / xi_circle.rC ;
594- double newX = ((particle.vx () - xi_circle.xC ) * std::cos (theta) - (particle.vy () - xi_circle.yC ) * std::sin (theta)) + xi_circle.xC ;
595- double newY = ((particle.vy () - xi_circle.yC ) * std::cos (theta) + (particle.vx () - xi_circle.xC ) * std::sin (theta)) + xi_circle.yC ;
596- double newPx = particle.px () * std::cos (theta) - particle.py () * std::sin (theta);
597- double newPy = particle.py () * std::cos (theta) + particle.px () * std::sin (theta);
592+ o2::math_utils::CircleXYf_t circleXi;
593+ track.getCircleParams (magneticField, circleXi, sna, csa);
594+ const double rxyXi = rxyzXi / std::sqrt (1 . + track.getTgl () * track.getTgl ());
595+ const double theta = rxyXi / circleXi.rC ;
596+ const double newX = ((particle.vx () - circleXi.xC ) * std::cos (theta) - (particle.vy () - circleXi.yC ) * std::sin (theta)) + circleXi.xC ;
597+ const double newY = ((particle.vy () - circleXi.yC ) * std::cos (theta) + (particle.vx () - circleXi.xC ) * std::sin (theta)) + circleXi.yC ;
598+ const double newPx = particle.px () * std::cos (theta) - particle.py () * std::sin (theta);
599+ const double newPy = particle.py () * std::cos (theta) + particle.px () * std::sin (theta);
600+ const double newE = std::sqrt (o2::constants::physics::MassXiMinus * o2::constants::physics::MassXiMinus + newPx * newPx + newPy * newPy + particle.pz () * particle.pz ());
601+
598602 xiDecayVertex.push_back (newX);
599603 xiDecayVertex.push_back (newY);
600- xiDecayVertex.push_back (particle.vz () + xi_rxyz * (particle.pz () / particle.p ()));
604+ xiDecayVertex.push_back (particle.vz () + rxyzXi * (particle.pz () / particle.p ()));
601605
602- std::vector<double > xiDaughters = {la_mass, pi_mass };
603- TLorentzVector xi (newPx, newPy, particle.pz (), particle. e () );
606+ std::vector<double > xiDaughters = {o2::constants::physics::MassLambda, o2::constants::physics::MassPionCharged };
607+ TLorentzVector xi (newPx, newPy, particle.pz (), newE );
604608 TGenPhaseSpace xiDecay;
605609 xiDecay.SetDecay (xi, 2 , xiDaughters.data ());
606610 xiDecay.Generate ();
607611 decayDaughters.push_back (*xiDecay.GetDecay (1 ));
608612 TLorentzVector la = *xiDecay.GetDecay (0 );
609613
610- double la_gamma = 1 / std::sqrt ( 1 + (la. P () * la. P ()) / (la_mass * la_mass) );
611- double la_ctau = 7.89 * la_gamma;
612- std::vector< double > laDaughters = {pi_mass, pr_mass} ;
613- double la_rxyz = (-la_ctau * std::log (1 - u ));
614- laDecayVertex.push_back (xiDecayVertex[0 ] + la_rxyz * (xiDecay.GetDecay (0 )->Px () / xiDecay.GetDecay (0 )->P ()));
615- laDecayVertex.push_back (xiDecayVertex[1 ] + la_rxyz * (xiDecay.GetDecay (0 )->Py () / xiDecay.GetDecay (0 )->P ()));
616- laDecayVertex.push_back (xiDecayVertex[2 ] + la_rxyz * (xiDecay.GetDecay (0 )->Pz () / xiDecay.GetDecay (0 )->P ()));
614+ const double uLa = rand. Uniform ( 0 , 1 );
615+ const double ctauLa = 7.845 ; // cm
616+ const double betaGammaLa = la. P () / o2::constants::physics::MassLambda ;
617+ const double rxyzLa = (-betaGammaLa * ctauLa * std::log (1 - uLa ));
618+ laDecayVertex.push_back (xiDecayVertex[0 ] + rxyzLa * (xiDecay.GetDecay (0 )->Px () / xiDecay.GetDecay (0 )->P ()));
619+ laDecayVertex.push_back (xiDecayVertex[1 ] + rxyzLa * (xiDecay.GetDecay (0 )->Py () / xiDecay.GetDecay (0 )->P ()));
620+ laDecayVertex.push_back (xiDecayVertex[2 ] + rxyzLa * (xiDecay.GetDecay (0 )->Pz () / xiDecay.GetDecay (0 )->P ()));
617621
622+ std::vector<double > laDaughters = {o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton};
618623 TGenPhaseSpace laDecay;
619624 laDecay.SetDecay (la, 2 , laDaughters.data ());
620625 laDecay.Generate ();
621626 decayDaughters.push_back (*laDecay.GetDecay (0 ));
622627 decayDaughters.push_back (*laDecay.GetDecay (1 ));
623628 }
629+
624630 // / Function to decay the V0
625631 // / \param particle the particle to decay
626632 // / \param decayDaughters the address of resulting daughters
@@ -629,36 +635,36 @@ struct OnTheFlyTracker {
629635 void decayV0Particle (McParticleType particle, std::vector<TLorentzVector>& decayDaughters, std::vector<double >& v0DecayVertex, int pdgCode)
630636 {
631637 double u = rand.Uniform (0 , 1 );
632- double v0_mass = -1 .;
633- double negDau_mass = -1 .;
634- double posDau_mass = -1 .;
638+ double v0Mass = -1 .;
639+ double negDauMass = -1 .;
640+ double posDauMass = -1 .;
635641 double ctau = -1 .;
642+
636643 if (std::abs (pdgCode) == kK0Short ) {
637- v0_mass = o2::constants::physics::MassK0Short;
638- negDau_mass = o2::constants::physics::MassPionCharged;
639- posDau_mass = o2::constants::physics::MassPionCharged;
644+ v0Mass = o2::constants::physics::MassK0Short;
645+ negDauMass = o2::constants::physics::MassPionCharged;
646+ posDauMass = o2::constants::physics::MassPionCharged;
640647 ctau = 2.68 ;
641648 } else if (std::abs (pdgCode) == kLambda0 ) {
642- v0_mass = o2::constants::physics::MassLambda;
643- negDau_mass = o2::constants::physics::MassPionCharged;
644- posDau_mass = o2::constants::physics::MassProton;
649+ v0Mass = o2::constants::physics::MassLambda;
650+ negDauMass = o2::constants::physics::MassPionCharged;
651+ posDauMass = o2::constants::physics::MassProton;
645652 ctau = 7.845 ;
646653 } else if (std::abs (pdgCode) == kLambda0Bar ) {
647- v0_mass = o2::constants::physics::MassLambda;
648- negDau_mass = o2::constants::physics::MassProton;
649- posDau_mass = o2::constants::physics::MassPionCharged;
654+ v0Mass = o2::constants::physics::MassLambda;
655+ negDauMass = o2::constants::physics::MassProton;
656+ posDauMass = o2::constants::physics::MassPionCharged;
650657 ctau = 7.845 ;
651658 }
652659
653- double v0_gamma = 1 / std::sqrt (1 + (particle.p () * particle.p ()) / (v0_mass * v0_mass));
654- double v0_ctau = ctau * v0_gamma;
655- double v0_rxyz = (-v0_ctau * std::log (1 - u));
660+ const double v0BetaGamma = particle.p () / v0Mass;
661+ const double v0rxyz = (-v0BetaGamma * ctau * std::log (1 - u));
656662 TLorentzVector v0 (particle.px (), particle.py (), particle.pz (), particle.e ());
657663
658- v0DecayVertex.push_back (particle.vx () + v0_rxyz * (particle.px () / particle.p ()));
659- v0DecayVertex.push_back (particle.vy () + v0_rxyz * (particle.py () / particle.p ()));
660- v0DecayVertex.push_back (particle.vz () + v0_rxyz * (particle.pz () / particle.p ()));
661- std::vector<double > v0Daughters = {negDau_mass, posDau_mass };
664+ v0DecayVertex.push_back (particle.vx () + v0rxyz * (particle.px () / particle.p ()));
665+ v0DecayVertex.push_back (particle.vy () + v0rxyz * (particle.py () / particle.p ()));
666+ v0DecayVertex.push_back (particle.vz () + v0rxyz * (particle.pz () / particle.p ()));
667+ std::vector<double > v0Daughters = {negDauMass, posDauMass };
662668
663669 TGenPhaseSpace v0Decay;
664670 v0Decay.SetDecay (v0, 2 , v0Daughters.data ());
@@ -744,7 +750,7 @@ struct OnTheFlyTracker {
744750 if (mcParticle.pdgCode () == kXiMinus ) {
745751 o2::track::TrackParCov xiTrackParCov;
746752 o2::upgrade::convertMCParticleToO2Track (mcParticle, xiTrackParCov, pdgDB);
747- decayParticle (mcParticle, xiTrackParCov, decayProducts, xiDecayVertex, laDecayVertex);
753+ decayCascade (mcParticle, xiTrackParCov, decayProducts, xiDecayVertex, laDecayVertex);
748754 xiDecayRadius2D = std::hypot (xiDecayVertex[0 ], xiDecayVertex[1 ]);
749755 laDecayRadius2D = std::hypot (laDecayVertex[0 ], laDecayVertex[1 ]);
750756 }
@@ -1214,22 +1220,28 @@ struct OnTheFlyTracker {
12141220 std::array{negP[0 ], negP[1 ], negP[2 ]}},
12151221 std::array{o2::constants::physics::MassPionCharged,
12161222 o2::constants::physics::MassPionCharged});
1217- } else
1223+ } else {
12181224 thisV0.mK0 = -1 ;
1225+ }
1226+
12191227 if (isLambda) {
12201228 thisV0.mLambda = RecoDecay::m (std::array{std::array{posP[0 ], posP[1 ], posP[2 ]},
12211229 std::array{negP[0 ], negP[1 ], negP[2 ]}},
12221230 std::array{o2::constants::physics::MassPionCharged,
12231231 o2::constants::physics::MassProton});
1224- } else
1232+ } else {
12251233 thisV0.mLambda = -1 ;
1234+ }
1235+
12261236 if (isAntiLambda) {
12271237 thisV0.mAntiLambda = RecoDecay::m (std::array{std::array{posP[0 ], posP[1 ], posP[2 ]},
12281238 std::array{negP[0 ], negP[1 ], negP[2 ]}},
12291239 std::array{o2::constants::physics::MassProton,
12301240 o2::constants::physics::MassPionCharged});
1231- } else
1241+ } else {
12321242 thisV0.mAntiLambda = -1 ;
1243+ }
1244+
12331245 if (v0DecaySettings.doV0QA ) {
12341246 histos.fill (HIST (" V0Building/hV0Building" ), 4 .0f );
12351247
0 commit comments