2525#include " Framework/runDataProcessing.h"
2626#include " ReconstructionDataFormats/PID.h"
2727
28+ #include < Math/GenVector/Boost.h>
29+ #include < Math/Vector4D.h>
30+ #include < TLorentzVector.h>
31+ #include < TMath.h>
32+
2833using namespace o2 ;
2934using namespace o2 ::framework;
3035using namespace o2 ::framework::expressions;
@@ -54,6 +59,7 @@ struct sigmaProtonCand {
5459 float sigmaDauPy = -1 ; // Py of the daughter track from Sigma decay
5560 float sigmaDauPz = -1 ; // Pz of the daughter track from Sigma decay
5661 float sigmaDecRadius = -1 ; // Decay radius of the Sigma candidate
62+ float sigmaCosPA = -1 ; // Cosine of pointing angle of the Sigma candidate
5763
5864 int chargePr = 0 ; // Charge of the proton candidate
5965 float pxPr = -1 ; // Px of the proton candidate
@@ -84,7 +90,7 @@ struct sigmaProtonCorrTask {
8490 Configurable<float > cutDCAtoPVSigma{" cutDCAtoPVSigma" , 0 .1f , " Max DCA to primary vertex for Sigma candidates (cm)" };
8591 Configurable<bool > doSigmaMinus{" doSigmaMinus" , true , " If true, pair Sigma- candidates, else Sigma+" };
8692 Configurable<float > cutSigmaRadius{" cutSigmaRadius" , 20 .f , " Minimum radius for Sigma candidates (cm)" };
87- Configurable<float > cutSigmaMass{" cutSigmaMass" , 0.3 , " Sigma mass window (MeV /c^2)" };
93+ Configurable<float > cutSigmaMass{" cutSigmaMass" , 0.3 , " Sigma mass window (GeV /c^2)" };
8894 Configurable<float > alphaAPCut{" alphaAPCut" , 0 ., " Alpha AP cut for Sigma candidates" };
8995 Configurable<float > qtAPCutLow{" qtAPCutLow" , 0.15 , " Lower qT AP cut for Sigma candidates (GeV/c)" };
9096 Configurable<float > qtAPCutHigh{" qtAPCutHigh" , 0.2 , " Upper qT AP cut for Sigma candidates (GeV/c)" };
@@ -94,6 +100,8 @@ struct sigmaProtonCorrTask {
94100 Configurable<float > cutNSigmaTPC{" cutNSigmaTPC" , 3 , " NSigmaTPCPr" };
95101 Configurable<float > cutNSigmaTOF{" cutNSigmaTOF" , 3 , " NSigmaTOFPr" };
96102
103+ Configurable<float > cutMaxKStar{" cutMaxKStar" , 1.5 , " Maximum k* for Sigma-Proton pairs (GeV/c)" };
104+
97105 Configurable<bool > fillOutputTree{" fillOutputTree" , true , " If true, fill the output tree with Sigma-Proton candidates" };
98106
99107 ConfigurableAxis CfgVtxBins{" CfgVtxBins" , {10 , -10 , 10 }, " Mixing bins - z-vertex" };
@@ -142,6 +150,36 @@ struct sigmaProtonCorrTask {
142150 return std::sqrt (p2A - dp * dp / p2V0);
143151 }
144152
153+ float getCosPA (const std::array<float , 3 >& momMother, const std::array<float , 3 >& decayVertex, const std::array<float , 3 >& primaryVertex)
154+ {
155+ std::array<float , 3 > decayVec = {decayVertex[0 ] - primaryVertex[0 ], decayVertex[1 ] - primaryVertex[1 ], decayVertex[2 ] - primaryVertex[2 ]};
156+ float dotProduct = std::inner_product (momMother.begin (), momMother.end (), decayVec.begin (), 0 .f );
157+ float momMotherMag = std::sqrt (std::inner_product (momMother.begin (), momMother.end (), momMother.begin (), 0 .f ));
158+ float decayVecMag = std::sqrt (std::inner_product (decayVec.begin (), decayVec.end (), decayVec.begin (), 0 .f ));
159+ return dotProduct / (momMotherMag * decayVecMag);
160+ }
161+
162+ TLorentzVector trackSum, PartOneCMS, PartTwoCMS, trackRelK;
163+ float getKStar (const sigmaProtonCand& candidate)
164+ {
165+ TLorentzVector part1; // Sigma
166+ TLorentzVector part2; // Proton
167+ part1.SetXYZM (candidate.sigmaPx , candidate.sigmaPy , candidate.sigmaPz , o2::constants::physics::MassSigmaMinus);
168+ part2.SetXYZM (candidate.pxPr , candidate.pyPr , candidate.pzPr , o2::constants::physics::MassProton);
169+ trackSum = part1 + part2;
170+ const float beta = trackSum.Beta ();
171+ const float betax = beta * std::cos (trackSum.Phi ()) * std::sin (trackSum.Theta ());
172+ const float betay = beta * std::sin (trackSum.Phi ()) * std::sin (trackSum.Theta ());
173+ const float betaz = beta * std::cos (trackSum.Theta ());
174+ PartOneCMS.SetXYZM (part1.Px (), part1.Py (), part1.Pz (), part1.M ());
175+ PartTwoCMS.SetXYZM (part2.Px (), part2.Py (), part2.Pz (), part2.M ());
176+ const ROOT::Math::Boost boostPRF = ROOT::Math::Boost (-betax, -betay, -betaz);
177+ PartOneCMS = boostPRF (PartOneCMS);
178+ PartTwoCMS = boostPRF (PartTwoCMS);
179+ trackRelK = PartOneCMS - PartTwoCMS;
180+ return 0.5 * trackRelK.P ();
181+ }
182+
145183 template <typename Ttrack>
146184 bool selectPrTrack (const Ttrack& candidate)
147185 {
@@ -211,8 +249,8 @@ struct sigmaProtonCorrTask {
211249 return true ;
212250 }
213251
214- template <typename Ttrack>
215- void fillTreeAndHistograms (aod::KinkCands const & kinkCands, Ttrack const & tracksDauSigma, Ttrack const & tracks)
252+ template <typename Ttrack, typename Tcollision >
253+ void fillTreeAndHistograms (aod::KinkCands const & kinkCands, Ttrack const & tracksDauSigma, Ttrack const & tracks, Tcollision const & collision )
216254 {
217255 for (const auto & sigmaCand : kinkCands) {
218256 if (selectSigma (sigmaCand, tracksDauSigma)) {
@@ -237,6 +275,11 @@ struct sigmaProtonCorrTask {
237275 candidate.sigmaDauPz = sigmaCand.pzDaug ();
238276 candidate.sigmaDecRadius = std::hypot (sigmaCand.xDecVtx (), sigmaCand.yDecVtx ());
239277
278+ std::array<float , 3 > momMoth = {sigmaCand.pxMoth (), sigmaCand.pyMoth (), sigmaCand.pzMoth ()};
279+ std::array<float , 3 > decayVtx = {sigmaCand.xDecVtx (), sigmaCand.yDecVtx (), sigmaCand.zDecVtx ()};
280+ std::array<float , 3 > primaryVtx = {collision.posX (), collision.posY (), collision.posZ ()};
281+ candidate.sigmaCosPA = getCosPA (momMoth, decayVtx, primaryVtx);
282+
240283 candidate.chargePr = prTrack.sign ();
241284 candidate.pxPr = prTrack.px ();
242285 candidate.pyPr = prTrack.py ();
@@ -249,6 +292,10 @@ struct sigmaProtonCorrTask {
249292 candidate.kinkDauID = sigmaCand.trackDaugId ();
250293 candidate.prID = prTrack.globalIndex ();
251294
295+ if (getKStar (candidate) > cutMaxKStar) {
296+ continue ;
297+ }
298+
252299 rSigmaProton.fill (HIST (" h2PtPrNSigma" ), candidate.ptPr (), candidate.nSigmaTPCPr );
253300 if (prTrack.hasTOF ()) {
254301 rSigmaProton.fill (HIST (" h2PtPrNSigmaTOF" ), candidate.ptPr (), candidate.nSigmaTOFPr );
@@ -270,7 +317,7 @@ struct sigmaProtonCorrTask {
270317 continue ;
271318 }
272319 rEventSelection.fill (HIST (" hVertexZRec" ), collision.posZ ());
273- fillTreeAndHistograms (kinkCands_c, tracks_c, tracks_c);
320+ fillTreeAndHistograms (kinkCands_c, tracks_c, tracks_c, collision );
274321 if (fillOutputTree) {
275322 // Fill output table
276323 for (const auto & candidate : sigmaProtonCandidates) {
@@ -282,6 +329,7 @@ struct sigmaProtonCorrTask {
282329 candidate.sigmaDauPy ,
283330 candidate.sigmaDauPz ,
284331 candidate.sigmaDecRadius ,
332+ candidate.sigmaCosPA ,
285333 candidate.chargePr ,
286334 candidate.pxPr ,
287335 candidate.pyPr ,
@@ -316,10 +364,7 @@ struct sigmaProtonCorrTask {
316364 auto kinkCands_c1 = kinkCands.sliceBy (kinkCandsPerCollisionPreslice, collision1.globalIndex ());
317365 auto tracks_c1 = tracks.sliceBy (tracksPerCollisionPreslice, collision1.globalIndex ());
318366 auto tracks_c2 = tracks.sliceBy (tracksPerCollisionPreslice, collision2.globalIndex ());
319- fillTreeAndHistograms (kinkCands_c1, tracks_c1, tracks_c2);
320-
321- auto kinkCands_c2 = kinkCands.sliceBy (kinkCandsPerCollisionPreslice, collision2.globalIndex ());
322- fillTreeAndHistograms (kinkCands_c2, tracks_c2, tracks_c1);
367+ fillTreeAndHistograms (kinkCands_c1, tracks_c1, tracks_c2, collision1);
323368
324369 if (fillOutputTree) {
325370 // Fill output table
@@ -332,6 +377,7 @@ struct sigmaProtonCorrTask {
332377 candidate.sigmaDauPy ,
333378 candidate.sigmaDauPz ,
334379 candidate.sigmaDecRadius ,
380+ candidate.sigmaCosPA ,
335381 candidate.chargePr ,
336382 candidate.pxPr ,
337383 candidate.pyPr ,
@@ -357,7 +403,7 @@ struct sigmaProtonCorrTask {
357403 continue ;
358404 }
359405 rEventSelection.fill (HIST (" hVertexZRec" ), collision.posZ ());
360- fillTreeAndHistograms (kinkCands_c, tracks_c, tracks_c);
406+ fillTreeAndHistograms (kinkCands_c, tracks_c, tracks_c, collision );
361407 if (fillOutputTree) {
362408 // Fill output table
363409 for (const auto & candidate : sigmaProtonCandidates) {
@@ -375,6 +421,7 @@ struct sigmaProtonCorrTask {
375421 candidate.sigmaDauPy ,
376422 candidate.sigmaDauPz ,
377423 candidate.sigmaDecRadius ,
424+ candidate.sigmaCosPA ,
378425 candidate.chargePr ,
379426 candidate.pxPr ,
380427 candidate.pyPr ,
@@ -407,10 +454,7 @@ struct sigmaProtonCorrTask {
407454 auto kinkCands_c1 = kinkCands.sliceBy (kinkCandsPerCollisionPreslice, collision1.globalIndex ());
408455 auto tracks_c1 = tracks.sliceBy (tracksPerCollisionPreslice, collision1.globalIndex ());
409456 auto tracks_c2 = tracks.sliceBy (tracksPerCollisionPreslice, collision2.globalIndex ());
410- fillTreeAndHistograms (kinkCands_c1, tracks_c1, tracks_c2);
411-
412- auto kinkCands_c2 = kinkCands.sliceBy (kinkCandsPerCollisionPreslice, collision2.globalIndex ());
413- fillTreeAndHistograms (kinkCands_c2, tracks_c2, tracks_c1);
457+ fillTreeAndHistograms (kinkCands_c1, tracks_c1, tracks_c2, collision1);
414458
415459 if (fillOutputTree) {
416460 // Fill output table
@@ -429,6 +473,7 @@ struct sigmaProtonCorrTask {
429473 candidate.sigmaDauPy ,
430474 candidate.sigmaDauPz ,
431475 candidate.sigmaDecRadius ,
476+ candidate.sigmaCosPA ,
432477 candidate.chargePr ,
433478 candidate.pxPr ,
434479 candidate.pyPr ,
0 commit comments