@@ -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 {
@@ -307,10 +307,23 @@ struct OnTheFlyTracker {
307307
308308 if (enablePrimarySmearing) {
309309 auto loadLUT = [&](int icfg, int pdg, const std::vector<std::string>& tables) {
310+ LOG (info) << " Loading LUT for pdg " << pdg << " for config " << icfg << " from provided tables with size " << tables.size ();
311+ if (tables.empty ()) {
312+ LOG (debug) << " No LUT file passed for pdg " << pdg << " , skipping." ;
313+ return false ;
314+ }
310315 const bool foundNewCfg = static_cast <size_t >(icfg) < tables.size ();
311- const std::string& lutFile = foundNewCfg ? tables[icfg] : tables.front ();
316+ std::string lutFile = foundNewCfg ? tables[icfg] : tables.front ();
317+ LOG (info) << " Loading LUT for pdg " << pdg << " from file " << lutFile << " for config " << icfg;
318+ // strip from leading/trailing spaces
319+ lutFile.erase (0 , lutFile.find_first_not_of (" " ));
320+ lutFile.erase (lutFile.find_last_not_of (" " ) + 1 );
321+ if (lutFile.empty ()) {
322+ LOG (debug) << " Empty LUT file name for pdg " << pdg << " , skipping." ;
323+ return false ;
324+ }
312325 bool success = mSmearer [icfg]->loadTable (pdg, lutFile.c_str ());
313- if (!success && !lutFile. empty () ) {
326+ if (!success) {
314327 LOG (fatal) << " Having issue with loading the LUT " << pdg << " " << lutFile;
315328 }
316329
@@ -575,52 +588,52 @@ struct OnTheFlyTracker {
575588 // / \param xiDecayVertex the address of the xi decay vertex
576589 // / \param laDecayVertex the address of the la decay vertex
577590 template <typename McParticleType>
578- void decayParticle (McParticleType particle, o2::track::TrackParCov track, std::vector<TLorentzVector>& decayDaughters, std::vector<double >& xiDecayVertex, std::vector<double >& laDecayVertex)
591+ void decayCascade (McParticleType particle, o2::track::TrackParCov track, std::vector<TLorentzVector>& decayDaughters, std::vector<double >& xiDecayVertex, std::vector<double >& laDecayVertex)
579592 {
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));
593+ const double uXi = rand.Uniform (0 , 1 );
594+ const double ctauXi = 4.91 ; // cm
595+ const double betaGammaXi = particle.p () / o2::constants::physics::MassXiMinus;
596+ const double rxyzXi = (-betaGammaXi * ctauXi * std::log (1 - uXi));
597+
589598 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);
599+ o2::math_utils::CircleXYf_t circleXi;
600+ track.getCircleParams (magneticField, circleXi, sna, csa);
601+ const double rxyXi = rxyzXi / std::sqrt (1 . + track.getTgl () * track.getTgl ());
602+ const double theta = rxyXi / circleXi.rC ;
603+ const double newX = ((particle.vx () - circleXi.xC ) * std::cos (theta) - (particle.vy () - circleXi.yC ) * std::sin (theta)) + circleXi.xC ;
604+ const double newY = ((particle.vy () - circleXi.yC ) * std::cos (theta) + (particle.vx () - circleXi.xC ) * std::sin (theta)) + circleXi.yC ;
605+ const double newPx = particle.px () * std::cos (theta) - particle.py () * std::sin (theta);
606+ const double newPy = particle.py () * std::cos (theta) + particle.px () * std::sin (theta);
607+ const double newE = std::sqrt (o2::constants::physics::MassXiMinus * o2::constants::physics::MassXiMinus + newPx * newPx + newPy * newPy + particle.pz () * particle.pz ());
608+
598609 xiDecayVertex.push_back (newX);
599610 xiDecayVertex.push_back (newY);
600- xiDecayVertex.push_back (particle.vz () + xi_rxyz * (particle.pz () / particle.p ()));
611+ xiDecayVertex.push_back (particle.vz () + rxyzXi * (particle.pz () / particle.p ()));
601612
602- std::vector<double > xiDaughters = {la_mass, pi_mass };
603- TLorentzVector xi (newPx, newPy, particle.pz (), particle. e () );
613+ std::vector<double > xiDaughters = {o2::constants::physics::MassLambda, o2::constants::physics::MassPionCharged };
614+ TLorentzVector xi (newPx, newPy, particle.pz (), newE );
604615 TGenPhaseSpace xiDecay;
605616 xiDecay.SetDecay (xi, 2 , xiDaughters.data ());
606617 xiDecay.Generate ();
607618 decayDaughters.push_back (*xiDecay.GetDecay (1 ));
608619 TLorentzVector la = *xiDecay.GetDecay (0 );
609620
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 ()));
621+ const double uLa = rand. Uniform ( 0 , 1 );
622+ const double ctauLa = 7.845 ; // cm
623+ const double betaGammaLa = la. P () / o2::constants::physics::MassLambda ;
624+ const double rxyzLa = (-betaGammaLa * ctauLa * std::log (1 - uLa ));
625+ laDecayVertex.push_back (xiDecayVertex[0 ] + rxyzLa * (xiDecay.GetDecay (0 )->Px () / xiDecay.GetDecay (0 )->P ()));
626+ laDecayVertex.push_back (xiDecayVertex[1 ] + rxyzLa * (xiDecay.GetDecay (0 )->Py () / xiDecay.GetDecay (0 )->P ()));
627+ laDecayVertex.push_back (xiDecayVertex[2 ] + rxyzLa * (xiDecay.GetDecay (0 )->Pz () / xiDecay.GetDecay (0 )->P ()));
617628
629+ std::vector<double > laDaughters = {o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton};
618630 TGenPhaseSpace laDecay;
619631 laDecay.SetDecay (la, 2 , laDaughters.data ());
620632 laDecay.Generate ();
621633 decayDaughters.push_back (*laDecay.GetDecay (0 ));
622634 decayDaughters.push_back (*laDecay.GetDecay (1 ));
623635 }
636+
624637 // / Function to decay the V0
625638 // / \param particle the particle to decay
626639 // / \param decayDaughters the address of resulting daughters
@@ -629,36 +642,36 @@ struct OnTheFlyTracker {
629642 void decayV0Particle (McParticleType particle, std::vector<TLorentzVector>& decayDaughters, std::vector<double >& v0DecayVertex, int pdgCode)
630643 {
631644 double u = rand.Uniform (0 , 1 );
632- double v0_mass = -1 .;
633- double negDau_mass = -1 .;
634- double posDau_mass = -1 .;
645+ double v0Mass = -1 .;
646+ double negDauMass = -1 .;
647+ double posDauMass = -1 .;
635648 double ctau = -1 .;
649+
636650 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;
651+ v0Mass = o2::constants::physics::MassK0Short;
652+ negDauMass = o2::constants::physics::MassPionCharged;
653+ posDauMass = o2::constants::physics::MassPionCharged;
640654 ctau = 2.68 ;
641655 } 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;
656+ v0Mass = o2::constants::physics::MassLambda;
657+ negDauMass = o2::constants::physics::MassPionCharged;
658+ posDauMass = o2::constants::physics::MassProton;
645659 ctau = 7.845 ;
646660 } 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;
661+ v0Mass = o2::constants::physics::MassLambda;
662+ negDauMass = o2::constants::physics::MassProton;
663+ posDauMass = o2::constants::physics::MassPionCharged;
650664 ctau = 7.845 ;
651665 }
652666
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));
667+ const double v0BetaGamma = particle.p () / v0Mass;
668+ const double v0rxyz = (-v0BetaGamma * ctau * std::log (1 - u));
656669 TLorentzVector v0 (particle.px (), particle.py (), particle.pz (), particle.e ());
657670
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 };
671+ v0DecayVertex.push_back (particle.vx () + v0rxyz * (particle.px () / particle.p ()));
672+ v0DecayVertex.push_back (particle.vy () + v0rxyz * (particle.py () / particle.p ()));
673+ v0DecayVertex.push_back (particle.vz () + v0rxyz * (particle.pz () / particle.p ()));
674+ std::vector<double > v0Daughters = {negDauMass, posDauMass };
662675
663676 TGenPhaseSpace v0Decay;
664677 v0Decay.SetDecay (v0, 2 , v0Daughters.data ());
@@ -744,7 +757,7 @@ struct OnTheFlyTracker {
744757 if (mcParticle.pdgCode () == kXiMinus ) {
745758 o2::track::TrackParCov xiTrackParCov;
746759 o2::upgrade::convertMCParticleToO2Track (mcParticle, xiTrackParCov, pdgDB);
747- decayParticle (mcParticle, xiTrackParCov, decayProducts, xiDecayVertex, laDecayVertex);
760+ decayCascade (mcParticle, xiTrackParCov, decayProducts, xiDecayVertex, laDecayVertex);
748761 xiDecayRadius2D = std::hypot (xiDecayVertex[0 ], xiDecayVertex[1 ]);
749762 laDecayRadius2D = std::hypot (laDecayVertex[0 ], laDecayVertex[1 ]);
750763 }
@@ -1214,22 +1227,28 @@ struct OnTheFlyTracker {
12141227 std::array{negP[0 ], negP[1 ], negP[2 ]}},
12151228 std::array{o2::constants::physics::MassPionCharged,
12161229 o2::constants::physics::MassPionCharged});
1217- } else
1230+ } else {
12181231 thisV0.mK0 = -1 ;
1232+ }
1233+
12191234 if (isLambda) {
12201235 thisV0.mLambda = RecoDecay::m (std::array{std::array{posP[0 ], posP[1 ], posP[2 ]},
12211236 std::array{negP[0 ], negP[1 ], negP[2 ]}},
12221237 std::array{o2::constants::physics::MassPionCharged,
12231238 o2::constants::physics::MassProton});
1224- } else
1239+ } else {
12251240 thisV0.mLambda = -1 ;
1241+ }
1242+
12261243 if (isAntiLambda) {
12271244 thisV0.mAntiLambda = RecoDecay::m (std::array{std::array{posP[0 ], posP[1 ], posP[2 ]},
12281245 std::array{negP[0 ], negP[1 ], negP[2 ]}},
12291246 std::array{o2::constants::physics::MassProton,
12301247 o2::constants::physics::MassPionCharged});
1231- } else
1248+ } else {
12321249 thisV0.mAntiLambda = -1 ;
1250+ }
1251+
12331252 if (v0DecaySettings.doV0QA ) {
12341253 histos.fill (HIST (" V0Building/hV0Building" ), 4 .0f );
12351254
0 commit comments