1818
1919#include " Common/Core/TableHelper.h"
2020
21+ #include " CommonConstants/MathConstants.h"
2122#include " CommonConstants/PhysicsConstants.h"
2223#include " Framework/InitContext.h"
2324
@@ -136,14 +137,15 @@ float qn(T const& col)
136137// / Recalculate pT for Kinks (Sigmas) using kinematic constraints
137138inline float calcPtnew (float pxMother, float pyMother, float pzMother, float pxDaughter, float pyDaughter, float pzDaughter)
138139{
140+ float almost0 = 1e-6f ;
139141 // Particle masses in GeV/c^2
140- const float massPion = 0 . 13957f ;
141- const float massNeutron = 0 . 93957f ;
142- const float massSigmaMinus = 1 . 19745f ;
142+ auto massPion = o2::constants::physics::MassPionCharged ;
143+ auto massNeutron = o2::constants::physics::MassNeutron ;
144+ auto massSigmaMinus = o2::constants::physics::MassSigmaMinus ;
143145
144146 // Calculate mother momentum and direction versor
145147 float pMother = std::sqrt (pxMother * pxMother + pyMother * pyMother + pzMother * pzMother);
146- if (pMother < 1e- 6f )
148+ if (pMother < almost0 )
147149 return -999 .f ;
148150
149151 float versorX = pxMother / pMother;
@@ -154,39 +156,39 @@ inline float calcPtnew(float pxMother, float pyMother, float pzMother, float pxD
154156 float ePi = std::sqrt (massPion * massPion + pxDaughter * pxDaughter + pyDaughter * pyDaughter + pzDaughter * pzDaughter);
155157
156158 // Scalar product of versor with daughter momentum
157- float a = versorX * pxDaughter + versorY * pyDaughter + versorZ * pzDaughter;
159+ float scalarProduct = versorX * pxDaughter + versorY * pyDaughter + versorZ * pzDaughter;
158160
159161 // Solve quadratic equation for momentum magnitude
160- float K = massSigmaMinus * massSigmaMinus + massPion * massPion - massNeutron * massNeutron;
161- float A = 4 .f * (ePi * ePi - a * a );
162- float B = -4 .f * a * K ;
163- float C = 4 .f * ePi * ePi * massSigmaMinus * massSigmaMinus - K * K ;
162+ float k = massSigmaMinus * massSigmaMinus + massPion * massPion - massNeutron * massNeutron;
163+ float a = 4 .f * (ePi * ePi - scalarProduct * scalarProduct );
164+ float b = -4 .f * scalarProduct * k ;
165+ float c = 4 .f * ePi * ePi * massSigmaMinus * massSigmaMinus - k * k ;
164166
165- if (std::abs (A ) < 1e- 6f )
167+ if (std::abs (a ) < almost0 )
166168 return -999 .f ;
167169
168- float D = B * B - 4 .f * A * C ;
169- if (D < 0 .f )
170+ float d = b * b - 4 .f * a * c ;
171+ if (d < 0 .f )
170172 return -999 .f ;
171173
172- float sqrtD = std::sqrt (D );
173- float P1 = (-B + sqrtD) / (2 .f * A );
174- float P2 = (-B - sqrtD) / (2 .f * A );
174+ float sqrtD = std::sqrt (d );
175+ float p1 = (-b + sqrtD) / (2 .f * a );
176+ float p2 = (-b - sqrtD) / (2 .f * a );
175177
176178 // Pick physical solution: prefer P2 if positive, otherwise P1
177- if (P2 < 0 .f && P1 < 0 .f )
179+ if (p2 < 0 .f && p1 < 0 .f )
178180 return -999 .f ;
179- if (P2 < 0 .f )
180- return P1 ;
181+ if (p2 < 0 .f )
182+ return p1 ;
181183
182184 // Choose solution closest to original momentum
183- float p1Diff = std::abs (P1 - pMother);
184- float p2Diff = std::abs (P2 - pMother);
185- float P = (p1Diff < p2Diff) ? P1 : P2 ;
185+ float p1Diff = std::abs (p1 - pMother);
186+ float p2Diff = std::abs (p2 - pMother);
187+ float p = (p1Diff < p2Diff) ? p1 : p2 ;
186188
187189 // Calculate pT from recalibrated momentum
188- float pxS = versorX * P ;
189- float pyS = versorY * P ;
190+ float pxS = versorX * p ;
191+ float pyS = versorY * p ;
190192 return std::sqrt (pxS * pxS + pyS * pyS);
191193}
192194
0 commit comments