2424#include " Framework/runDataProcessing.h"
2525
2626#include < random>
27- #include " TLorentzVector .h"
27+ #include " Math/Vector4D .h"
2828
2929#include " Common/DataModel/PIDResponse.h"
3030
@@ -116,7 +116,9 @@ struct UpcRhoAnalysis {
116116 Produces<o2::aod::RecoTree> recoTree;
117117 Produces<o2::aod::McTree> mcTree;
118118
119- float pcEtaCut = 0.9 ; // physics coordination recommendation
119+ float pcEtaCut = 0.9 ; // physics coordination recommendation
120+ const int gapSide = 2 ; // required gap side
121+ const int piPDG = 211 ; // PDG code for pion
120122 Configurable<int > numPions{" numPions" , 2 , " required number of pions in the event" };
121123 Configurable<bool > requireTof{" requireTof" , false , " require TOF signal" };
122124 Configurable<bool > onlyGoldenRuns{" onlyGoldenRuns" , false , " process only golden runs" };
@@ -387,18 +389,12 @@ struct UpcRhoAnalysis {
387389 template <typename C>
388390 bool collisionPassesCuts (const C& collision) // collision cuts
389391 {
392+ if (!collision.vtxITSTPC () || !collision.sbp () || !collision.itsROFb () || !collision.tfb ()) // not applied automatically in pass5
393+ return false ;
390394 if (std::abs (collision.posZ ()) > collisionsPosZMaxCut)
391395 return false ;
392396 if (collision.numContrib () > collisionsNumContribsMaxCut)
393397 return false ;
394- if (!collision.vtxITSTPC ())
395- return false ;
396- if (!collision.sbp ())
397- return false ;
398- if (!collision.itsROFb ())
399- return false ;
400- if (!collision.tfb ())
401- return false ;
402398 return true ;
403399 }
404400
@@ -505,7 +501,7 @@ struct UpcRhoAnalysis {
505501 return charge;
506502 }
507503
508- bool systemPassesCuts (const TLorentzVector & system) // system cuts
504+ bool systemPassesCuts (const ROOT::Math::PxPyPzMVector & system) // system cuts
509505 {
510506 if (system.M () < systemMassMinCut || system.M () > systemMassMaxCut)
511507 return false ;
@@ -516,47 +512,57 @@ struct UpcRhoAnalysis {
516512 return true ;
517513 }
518514
519- TLorentzVector reconstructSystem (const std::vector<TLorentzVector >& cutTracksLVs) // reconstruct system from 4-vectors
515+ ROOT::Math::PxPyPzMVector reconstructSystem (const std::vector<ROOT::Math::PxPyPzMVector >& cutTracksLVs) // reconstruct system from 4-vectors
520516 {
521- TLorentzVector system;
517+ ROOT::Math::PxPyPzMVector system;
522518 for (const auto & trackLV : cutTracksLVs)
523519 system += trackLV;
524520 return system;
525521 }
526522
527- float getPhiRandom (const std::vector<TLorentzVector>& cutTracksLVs) // decay phi anisotropy
528- { // two possible definitions of phi: randomize the tracks
523+ double deltaPhi (const ROOT::Math::PxPyPzMVector& p1, const ROOT::Math::PxPyPzMVector& p2)
524+ {
525+ double dPhi = p1.Phi () - p2.Phi ();
526+ if (dPhi > o2::constants::math::PI)
527+ dPhi -= o2::constants::math::TwoPI;
528+ else if (dPhi < -o2::constants::math::PI)
529+ dPhi += o2::constants::math::TwoPI;
530+ return dPhi; // calculate delta phi in (-pi, pi)
531+ }
532+
533+ float getPhiRandom (const std::vector<ROOT::Math::PxPyPzMVector>& cutTracksLVs) // decay phi anisotropy
534+ { // two possible definitions of phi: randomize the tracks
529535 int indices[2 ] = {0 , 1 };
530- unsigned seed = std::chrono::system_clock::now ().time_since_epoch ().count (); // get time-based seed
536+ unsigned seed = std::chrono::system_clock::now ().time_since_epoch ().count (); // get time-based seed
531537 std::shuffle (std::begin (indices), std::end (indices), std::default_random_engine (seed)); // shuffle indices
532538 // calculate phi
533- TLorentzVector pOne = cutTracksLVs[indices[0 ]];
534- TLorentzVector pTwo = cutTracksLVs[indices[1 ]];
535- TLorentzVector pPlus = pOne + pTwo;
536- TLorentzVector pMinus = pOne - pTwo;
537- return pPlus. DeltaPhi ( pMinus);
539+ ROOT::Math::PxPyPzMVector pOne = cutTracksLVs[indices[0 ]];
540+ ROOT::Math::PxPyPzMVector pTwo = cutTracksLVs[indices[1 ]];
541+ ROOT::Math::PxPyPzMVector pPlus = pOne + pTwo;
542+ ROOT::Math::PxPyPzMVector pMinus = pOne - pTwo;
543+ return deltaPhi (pPlus, pMinus);
538544 }
539545
540546 template <typename T>
541- float getPhiCharge (const T& cutTracks, const std::vector<TLorentzVector >& cutTracksLVs)
547+ float getPhiCharge (const T& cutTracks, const std::vector<ROOT::Math::PxPyPzMVector >& cutTracksLVs)
542548 { // two possible definitions of phi: charge-based assignment
543- TLorentzVector pOne, pTwo;
549+ ROOT::Math::PxPyPzMVector pOne, pTwo;
544550 pOne = (cutTracks[0 ].sign () > 0 ) ? cutTracksLVs[0 ] : cutTracksLVs[1 ];
545551 pTwo = (cutTracks[0 ].sign () > 0 ) ? cutTracksLVs[1 ] : cutTracksLVs[0 ];
546- TLorentzVector pPlus = pOne + pTwo;
547- TLorentzVector pMinus = pOne - pTwo;
548- return pPlus. DeltaPhi ( pMinus);
552+ ROOT::Math::PxPyPzMVector pPlus = pOne + pTwo;
553+ ROOT::Math::PxPyPzMVector pMinus = pOne - pTwo;
554+ return deltaPhi (pPlus, pMinus);
549555 }
550556
551557 template <typename T>
552- float getPhiChargeMC (const T& cutTracks, const std::vector<TLorentzVector >& cutTracksLVs)
558+ float getPhiChargeMC (const T& cutTracks, const std::vector<ROOT::Math::PxPyPzMVector >& cutTracksLVs)
553559 { // the same as for data but using pdg code instead of charge
554- TLorentzVector pOne, pTwo;
560+ ROOT::Math::PxPyPzMVector pOne, pTwo;
555561 pOne = (cutTracks[0 ].pdgCode () > 0 ) ? cutTracksLVs[0 ] : cutTracksLVs[1 ];
556562 pTwo = (cutTracks[0 ].pdgCode () > 0 ) ? cutTracksLVs[1 ] : cutTracksLVs[0 ];
557- TLorentzVector pPlus = pOne + pTwo;
558- TLorentzVector pMinus = pOne - pTwo;
559- return pPlus. DeltaPhi ( pMinus);
563+ ROOT::Math::PxPyPzMVector pPlus = pOne + pTwo;
564+ ROOT::Math::PxPyPzMVector pMinus = pOne - pTwo;
565+ return deltaPhi (pPlus, pMinus);
560566 }
561567
562568 template <typename C, typename T>
@@ -583,20 +589,20 @@ struct UpcRhoAnalysis {
583589 if (std::isinf (timeZNC))
584590 timeZNC = -999 ;
585591
586- if (energyCommonZNA <= znCommonEnergyCut && energyCommonZNC <= znCommonEnergyCut){
592+ if (energyCommonZNA <= znCommonEnergyCut && energyCommonZNC <= znCommonEnergyCut) {
587593 onon = true ;
588594 neutronClass = 0 ;
589595 }
590- if (energyCommonZNA > znCommonEnergyCut && std::abs (timeZNA) <= znTimeCut && energyCommonZNC <= znCommonEnergyCut){
596+ if (energyCommonZNA > znCommonEnergyCut && std::abs (timeZNA) <= znTimeCut && energyCommonZNC <= znCommonEnergyCut) {
591597 xnon = true ;
592598 neutronClass = 1 ;
593599 }
594- if (energyCommonZNA <= znCommonEnergyCut && energyCommonZNC > znCommonEnergyCut && std::abs (timeZNC) <= znTimeCut){
600+ if (energyCommonZNA <= znCommonEnergyCut && energyCommonZNC > znCommonEnergyCut && std::abs (timeZNC) <= znTimeCut) {
595601 onxn = true ;
596602 neutronClass = 2 ;
597603 }
598604 if (energyCommonZNA > znCommonEnergyCut && std::abs (timeZNA) <= znTimeCut &&
599- energyCommonZNC > znCommonEnergyCut && std::abs (timeZNC) <= znTimeCut){
605+ energyCommonZNC > znCommonEnergyCut && std::abs (timeZNC) <= znTimeCut) {
600606 xnxn = true ;
601607 neutronClass = 3 ;
602608 }
@@ -624,11 +630,9 @@ struct UpcRhoAnalysis {
624630 rQC.fill (HIST (" QC/tracks/trackSelections/hTpcNSigmaKa2D" ), cutTracks[0 ].tpcNSigmaKa (), cutTracks[1 ].tpcNSigmaKa ());
625631
626632 // create a vector of 4-vectors for selected tracks
627- std::vector<TLorentzVector > cutTracksLVs;
633+ std::vector<ROOT::Math::PxPyPzMVector > cutTracksLVs;
628634 for (const auto & cutTrack : cutTracks) {
629- TLorentzVector cutTrackLV;
630- cutTrackLV.SetXYZM (cutTrack.px (), cutTrack.py (), cutTrack.pz (), o2::constants::physics::MassPionCharged); // apriori assume pion mass
631- cutTracksLVs.push_back (cutTrackLV);
635+ cutTracksLVs.push_back (ROOT::Math::PxPyPzMVector (cutTrack.px (), cutTrack.py (), cutTrack.pz (), o2::constants::physics::MassPionCharged)); // apriori assume pion mass
632636 }
633637
634638 // differentiate leading- and subleading-momentum tracks
@@ -659,7 +663,7 @@ struct UpcRhoAnalysis {
659663 float trackDcaXYs[2 ] = {positiveTrack.dcaXY (), negativeTrack.dcaXY ()};
660664 float trackDcaZs[2 ] = {positiveTrack.dcaZ (), negativeTrack.dcaZ ()};
661665 float trackTpcSignals[2 ] = {positiveTrack.tpcSignal (), negativeTrack.tpcSignal ()};
662- recoTree (collision.flags (),collision.runNumber (), localBc, collision.numContrib (), collision.posX (), collision.posY (), collision.posZ (),
666+ recoTree (collision.flags (), collision.runNumber (), localBc, collision.numContrib (), collision.posX (), collision.posY (), collision.posZ (),
663667 collision.totalFT0AmplitudeA (), collision.totalFT0AmplitudeC (), collision.totalFV0AmplitudeA (), collision.totalFDDAmplitudeA (), collision.totalFDDAmplitudeC (),
664668 collision.timeFT0A (), collision.timeFT0C (), collision.timeFV0A (), collision.timeFDDA (), collision.timeFDDC (),
665669 energyCommonZNA, energyCommonZNC, timeZNA, timeZNC, neutronClass,
@@ -676,7 +680,7 @@ struct UpcRhoAnalysis {
676680 rQC.fill (HIST (" QC/tracks/hTofHitCheck" ), leadingMomentumTrack.hasTOF (), subleadingMomentumTrack.hasTOF ());
677681 fillCollisionQcHistos<1 >(collision); // fill QC histograms after track selections
678682
679- TLorentzVector system = reconstructSystem (cutTracksLVs);
683+ ROOT::Math::PxPyPzMVector system = reconstructSystem (cutTracksLVs);
680684 int totalCharge = tracksTotalCharge (cutTracks);
681685 float mass = system.M ();
682686 float pT = system.Pt ();
@@ -764,7 +768,7 @@ struct UpcRhoAnalysis {
764768 rMC.fill (HIST (" MC/collisions/hPosZ" ), mcCollision.posZ ());
765769
766770 std::vector<decltype (mcParticles.begin ())> cutMcParticles;
767- std::vector<TLorentzVector > mcParticlesLVs;
771+ std::vector<ROOT::Math::PxPyPzMVector > mcParticlesLVs;
768772
769773 for (auto const & mcParticle : mcParticles) {
770774 rMC.fill (HIST (" MC/tracks/all/hPdgCode" ), mcParticle.pdgCode ());
@@ -773,10 +777,10 @@ struct UpcRhoAnalysis {
773777 rMC.fill (HIST (" MC/tracks/all/hPt" ), pt (mcParticle.px (), mcParticle.py ()));
774778 rMC.fill (HIST (" MC/tracks/all/hEta" ), eta (mcParticle.px (), mcParticle.py (), mcParticle.pz ()));
775779 rMC.fill (HIST (" MC/tracks/all/hPhi" ), phi (mcParticle.px (), mcParticle.py ()));
776- if (!mcParticle.isPhysicalPrimary () || std::abs (mcParticle.pdgCode ()) != 211 )
780+ if (!mcParticle.isPhysicalPrimary () || std::abs (mcParticle.pdgCode ()) != piPDG )
777781 continue ;
778782 cutMcParticles.push_back (mcParticle);
779- TLorentzVector pionLV;
783+ ROOT::Math::PxPyPzMVector pionLV;
780784 pionLV.SetPxPyPzE (mcParticle.px (), mcParticle.py (), mcParticle.pz (), mcParticle.e ());
781785 mcParticlesLVs.push_back (pionLV);
782786 }
@@ -789,7 +793,7 @@ struct UpcRhoAnalysis {
789793 if (tracksTotalChargeMC (cutMcParticles) != 0 ) // shouldn't happen in theory
790794 return ;
791795
792- TLorentzVector system = reconstructSystem (mcParticlesLVs);
796+ ROOT::Math::PxPyPzMVector system = reconstructSystem (mcParticlesLVs);
793797 float mass = system.M ();
794798 float pT = system.Pt ();
795799 float rapidity = system.Rapidity ();
@@ -802,7 +806,7 @@ struct UpcRhoAnalysis {
802806 rMC.fill (HIST (" MC/tracks/hPt" ), pt (leadingMomentumPion.px (), leadingMomentumPion.py ()), pt (subleadingMomentumPion.px (), subleadingMomentumPion.py ()));
803807 rMC.fill (HIST (" MC/tracks/hEta" ), eta (leadingMomentumPion.px (), leadingMomentumPion.py (), leadingMomentumPion.pz ()), eta (subleadingMomentumPion.px (), subleadingMomentumPion.py (), subleadingMomentumPion.pz ()));
804808 rMC.fill (HIST (" MC/tracks/hPhi" ), phi (leadingMomentumPion.px (), leadingMomentumPion.py ()), phi (subleadingMomentumPion.px (), subleadingMomentumPion.py ()));
805-
809+
806810 rMC.fill (HIST (" MC/system/hM" ), mass);
807811 rMC.fill (HIST (" MC/system/hPt" ), pT);
808812 rMC.fill (HIST (" MC/system/hPtVsM" ), mass, pT);
@@ -813,7 +817,7 @@ struct UpcRhoAnalysis {
813817 rMC.fill (HIST (" MC/system/hPhiCharge" ), phiCharge);
814818 rMC.fill (HIST (" MC/system/hPhiRandomVsM" ), mass, phiRandom);
815819 rMC.fill (HIST (" MC/system/hPhiChargeVsM" ), mass, phiCharge);
816-
820+
817821 if (systemPassesCuts (system)) {
818822 rMC.fill (HIST (" MC/system/selected/hM" ), mass);
819823 rMC.fill (HIST (" MC/system/selected/hPt" ), pT);
@@ -826,7 +830,7 @@ struct UpcRhoAnalysis {
826830 rMC.fill (HIST (" MC/system/selected/hPhiRandomVsM" ), mass, phiRandom);
827831 rMC.fill (HIST (" MC/system/selected/hPhiChargeVsM" ), mass, phiCharge);
828832 }
829-
833+
830834 // fill mcTree
831835 auto positivePion = cutMcParticles[0 ].pdgCode () > 0 ? cutMcParticles[0 ] : cutMcParticles[1 ];
832836 auto negativePion = cutMcParticles[0 ].pdgCode () > 0 ? cutMcParticles[1 ] : cutMcParticles[0 ];
@@ -849,7 +853,7 @@ struct UpcRhoAnalysis {
849853
850854 void processSGdata (FullUdSgCollision const & collision, FullUdTracks const & tracks)
851855 {
852- if (collision.gapSide () != 2 )
856+ if (collision.gapSide () != gapSide )
853857 return ;
854858 processReco (collision, tracks);
855859 }
0 commit comments