@@ -133,6 +133,96 @@ float qn(T const& col)
133133 return qn;
134134}
135135
136+ // / Recalculate pT for Kinks (Sigmas) using kinematic constraints
137+ inline float calcPtnew (float pxMother, float pyMother, float pzMother, float pxDaughter, float pyDaughter, float pzDaughter)
138+ {
139+ // Particle masses in GeV/c^2
140+ const float massPion = 0 .13957f ;
141+ const float massNeutron = 0 .93957f ;
142+ const float massSigmaMinus = 1 .19745f ;
143+
144+ // Calculate mother momentum and direction versor
145+ float pMother = std::sqrt (pxMother * pxMother + pyMother * pyMother + pzMother * pzMother);
146+ if (pMother < 1e-6f )
147+ return -999 .f ;
148+
149+ float versorX = pxMother / pMother;
150+ float versorY = pyMother / pMother;
151+ float versorZ = pzMother / pMother;
152+
153+ // Calculate daughter energy
154+ float ePi = std::sqrt (massPion * massPion + pxDaughter * pxDaughter + pyDaughter * pyDaughter + pzDaughter * pzDaughter);
155+
156+ // Scalar product of versor with daughter momentum
157+ float a = versorX * pxDaughter + versorY * pyDaughter + versorZ * pzDaughter;
158+
159+ // 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;
164+
165+ if (std::abs (A) < 1e-6f )
166+ return -999 .f ;
167+
168+ float D = B * B - 4 .f * A * C;
169+ if (D < 0 .f )
170+ return -999 .f ;
171+
172+ float sqrtD = std::sqrt (D);
173+ float P1 = (-B + sqrtD) / (2 .f * A);
174+ float P2 = (-B - sqrtD) / (2 .f * A);
175+
176+ // Pick physical solution: prefer P2 if positive, otherwise P1
177+ if (P2 < 0 .f && P1 < 0 .f )
178+ return -999 .f ;
179+ if (P2 < 0 .f )
180+ return P1;
181+
182+ // 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;
186+
187+ // Calculate pT from recalibrated momentum
188+ float pxS = versorX * P;
189+ float pyS = versorY * P;
190+ return std::sqrt (pxS * pxS + pyS * pyS);
191+ }
192+
193+ // / Helper function to calculate recalculated pT for kink particles (Sigma/SigmaPlus)
194+ template <typename TParticle, typename TTrackTable>
195+ float getRecalculatedPtForKink (const TParticle& particle, const TTrackTable& trackTable)
196+ {
197+ // Check if particle has chaDau index
198+ if constexpr (requires { particle.has_chaDau (); }) {
199+ if (particle.has_chaDau ()) {
200+ try {
201+ auto chaDaughter = trackTable.rawIteratorAt (particle.chaDauId () - trackTable.offset ());
202+
203+ // Extract momentum components directly from dynamic columns
204+ float pxDaug = chaDaughter.px ();
205+ float pyDaug = chaDaughter.py ();
206+ float pzDaug = chaDaughter.pz ();
207+
208+ // Get momentum components from dynamic columns
209+ float pxMoth = particle.px ();
210+ float pyMoth = particle.py ();
211+ float pzMoth = particle.pz ();
212+
213+ // Recalculate pT using kinematic constraints
214+ float ptRecalc = calcPtnew (pxMoth, pyMoth, pzMoth, pxDaug, pyDaug, pzDaug);
215+ if (ptRecalc > 0 ) {
216+ return ptRecalc;
217+ }
218+ } catch (const std::exception& e) {
219+ return -1 .0f ;
220+ }
221+ }
222+ }
223+ return -1 .0f ;
224+ }
225+
136226inline bool enableTable (const char * tableName, int userSetting, o2::framework::InitContext& initContext)
137227{
138228 if (userSetting == 1 ) {
0 commit comments