@@ -170,6 +170,11 @@ struct HigherMassResonances {
170170 // ConfigurableAxis axisdEdx{"axisdEdx", {20000, 0.0f, 200.0f}, "dE/dx (a.u.)"};
171171 // ConfigurableAxis axisPtfordEbydx{"axisPtfordEbydx", {2000, 0, 20}, "pT (GeV/c)"};
172172 // ConfigurableAxis axisMultdist{"axisMultdist", {3500, 0, 70000}, "Multiplicity distribution"};
173+
174+ // fixed variables
175+ float rapidityMotherData = 0.5 ;
176+ std::array<int , 6 > numbers = {0 , 1 , 2 , 3 , 4 , 5 };
177+ double beamMomentum = std::sqrt(13600 * 13600 / 4 - o2::constants::physics::MassProton * o2::constants::physics::MassProton); // GeV
173178 } config;
174179
175180 // Service<o2::framework::O2DatabasePDG> PDGdatabase;
@@ -183,10 +188,9 @@ struct HigherMassResonances {
183188 ROOT::Math::XYZVector randomVec, beamVec, normalVec;
184189 ROOT::Math::XYZVectorF v1CM, zaxisHE, yaxisHE, xaxisHE;
185190 // ROOT::Math::XYZVector threeVecDauCM, helicityVec, randomVec, beamVec, normalVec;
186- ROOT::Math::XYZVector zBeam; // ẑ: beam direction in lab frame
187- double beamMomentum = std::sqrt(13600 * 13600 / 4 - 0.938 * 0.938 ); // GeV
188- ROOT::Math::PxPyPzEVector beam1{0 ., 0 ., -beamMomentum, 13600 . / 2 .};
189- ROOT::Math::PxPyPzEVector beam2{0 ., 0 ., beamMomentum, 13600 . / 2 .};
191+ ROOT::Math::XYZVector zBeam; // ẑ: beam direction in lab frame
192+ ROOT::Math::PxPyPzEVector beam1{0 ., 0 ., -config.beamMomentum , 13600 . / 2 .};
193+ ROOT::Math::PxPyPzEVector beam2{0 ., 0 ., config.beamMomentum , 13600 . / 2 .};
190194 ROOT::Math::XYZVectorF beam1CM, beam2CM;
191195
192196 // const double massK0s = o2::constants::physics::MassK0Short;
@@ -388,8 +392,8 @@ struct HigherMassResonances {
388392 const float cpav0 = candidate.v0cosPA ();
389393
390394 float ctauK0s = candidate.distovertotmom (collision.posX (), collision.posY (), collision.posZ ()) * o2::constants::physics::MassK0Short;
391- float lowmasscutks0 = 0.497 - config.cWidthKs0 * config.cSigmaMassKs0 ;
392- float highmasscutks0 = 0.497 + config.cWidthKs0 * config.cSigmaMassKs0 ;
395+ float lowmasscutks0 = o2::constants::physics::MassKPlus - config.cWidthKs0 * config.cSigmaMassKs0 ;
396+ float highmasscutks0 = o2::constants::physics::MassKPlus + config.cWidthKs0 * config.cSigmaMassKs0 ;
393397 // float decayLength = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * RecoDecay::sqrtSumOfSquares(candidate.px(), candidate.py(), candidate.pz());
394398
395399 if (config.qAv0 ) {
@@ -647,17 +651,18 @@ struct HigherMassResonances {
647651 // CosThetaHE = zaxisHE.Dot(v_CM);
648652
649653 auto anglePhi = std::atan2 (yaxisHE.Dot (v1CM), xaxisHE.Dot (v1CM));
650- if (anglePhi < 0 ) {
651- anglePhi += 2 * o2::constants::math::PI; // ensure phi is in [0, 2pi]
652- }
654+ anglePhi = RecoDecay::constrainAngle (anglePhi, 0.0 );
655+ // if (anglePhi < 0) {
656+ // anglePhi += o2::constants::math::TwoPI; // ensure phi is in [0, 2pi]
657+ // }
653658
654- // if (std::abs(mother.Rapidity()) < 0.5 ) {
659+ // if (std::abs(mother.Rapidity()) < config.rapidityMotherData ) {
655660 if (config.activateTHnSparseCosThStarHelicity ) {
656661 // helicityVec = mother.Vect(); // 3 vector of mother in COM frame
657662 // auto cosThetaStarHelicity = helicityVec.Dot(threeVecDauCM) / (std::sqrt(threeVecDauCM.Mag2()) * std::sqrt(helicityVec.Mag2()));
658663 auto cosThetaStarHelicity = mother.Vect ().Dot (fourVecDauCM.Vect ()) / (std::sqrt (fourVecDauCM.Vect ().Mag2 ()) * std::sqrt (mother.Vect ().Mag2 ()));
659664 if (!isMix) {
660- if (std::abs (mother.Rapidity ()) < 0.5 ) {
665+ if (std::abs (mother.Rapidity ()) < config. rapidityMotherData ) {
661666 hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, mother.Pt (), mother.M (), cosThetaStarHelicity, anglePhi);
662667 }
663668
@@ -672,49 +677,49 @@ struct HigherMassResonances {
672677 daughterRotCM = boost2 (daughterRot);
673678
674679 auto cosThetaStarHelicityRot = motherRot.Vect ().Dot (daughterRotCM.Vect ()) / (std::sqrt (daughterRotCM.Vect ().Mag2 ()) * std::sqrt (motherRot.Vect ().Mag2 ()));
675- if (motherRot.Rapidity () < 0.5 )
680+ if (motherRot.Rapidity () < config. rapidityMotherData )
676681 hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, motherRot.Pt (), motherRot.M (), cosThetaStarHelicityRot, anglePhi);
677682 }
678683 } else {
679- if (std::abs (mother.Rapidity ()) < 0.5 ) {
684+ if (std::abs (mother.Rapidity ()) < config. rapidityMotherData ) {
680685 hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, mother.Pt (), mother.M (), cosThetaStarHelicity, anglePhi);
681686 }
682687 }
683688 } else if (config.activateTHnSparseCosThStarProduction ) {
684689 normalVec = ROOT::Math::XYZVector (mother.Py (), -mother.Px (), 0 .f );
685690 auto cosThetaStarProduction = normalVec.Dot (fourVecDauCM.Vect ()) / (std::sqrt (fourVecDauCM.Vect ().Mag2 ()) * std::sqrt (normalVec.Mag2 ()));
686691 if (!isMix) {
687- if (std::abs (mother.Rapidity ()) < 0.5 ) {
692+ if (std::abs (mother.Rapidity ()) < config. rapidityMotherData ) {
688693 hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, mother.Pt (), mother.M (), cosThetaStarProduction, anglePhi);
689694 }
690695 for (int i = 0 ; i < config.cRotations ; i++) {
691696 theta2 = rn->Uniform (o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut , o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut );
692697 motherRot = ROOT::Math::PxPyPzMVector (mother.Px () * std::cos (theta2) - mother.Py () * std::sin (theta2), mother.Px () * std::sin (theta2) + mother.Py () * std::cos (theta2), mother.Pz (), mother.M ());
693- if (std::abs (motherRot.Rapidity ()) < 0.5 ) {
698+ if (std::abs (motherRot.Rapidity ()) < config. rapidityMotherData ) {
694699 hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, motherRot.Pt (), motherRot.M (), cosThetaStarProduction, anglePhi);
695700 }
696701 }
697702 } else {
698- if (std::abs (mother.Rapidity ()) < 0.5 ) {
703+ if (std::abs (mother.Rapidity ()) < config. rapidityMotherData ) {
699704 hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, mother.Pt (), mother.M (), cosThetaStarProduction, anglePhi);
700705 }
701706 }
702707 } else if (config.activateTHnSparseCosThStarBeam ) {
703708 beamVec = ROOT::Math::XYZVector (0 .f , 0 .f , 1 .f );
704709 auto cosThetaStarBeam = beamVec.Dot (fourVecDauCM.Vect ()) / std::sqrt (fourVecDauCM.Vect ().Mag2 ());
705710 if (!isMix) {
706- if (std::abs (mother.Rapidity ()) < 0.5 ) {
711+ if (std::abs (mother.Rapidity ()) < config. rapidityMotherData ) {
707712 hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, mother.Pt (), mother.M (), cosThetaStarBeam, anglePhi);
708713 }
709714 for (int i = 0 ; i < config.cRotations ; i++) {
710715 theta2 = rn->Uniform (o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut , o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut );
711716 motherRot = ROOT::Math::PxPyPzMVector (mother.Px () * std::cos (theta2) - mother.Py () * std::sin (theta2), mother.Px () * std::sin (theta2) + mother.Py () * std::cos (theta2), mother.Pz (), mother.M ());
712- if (std::abs (motherRot.Rapidity ()) < 0.5 ) {
717+ if (std::abs (motherRot.Rapidity ()) < config. rapidityMotherData ) {
713718 hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, motherRot.Pt (), motherRot.M (), cosThetaStarBeam, anglePhi);
714719 }
715720 }
716721 } else {
717- if (std::abs (mother.Rapidity ()) < 0.5 ) {
722+ if (std::abs (mother.Rapidity ()) < config. rapidityMotherData ) {
718723 hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, mother.Pt (), mother.M (), cosThetaStarBeam, anglePhi);
719724 }
720725 }
@@ -725,18 +730,18 @@ struct HigherMassResonances {
725730 randomVec = ROOT::Math::XYZVector (std::sin (thetaRandom) * std::cos (phiRandom), std::sin (thetaRandom) * std::sin (phiRandom), std::cos (thetaRandom));
726731 auto cosThetaStarRandom = randomVec.Dot (fourVecDauCM.Vect ()) / std::sqrt (fourVecDauCM.Vect ().Mag2 ());
727732 if (!isMix) {
728- if (std::abs (mother.Rapidity ()) < 0.5 ) {
733+ if (std::abs (mother.Rapidity ()) < config. rapidityMotherData ) {
729734 hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, mother.Pt (), mother.M (), cosThetaStarRandom, anglePhi);
730735 }
731736 for (int i = 0 ; i < config.cRotations ; i++) {
732737 theta2 = rn->Uniform (o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut , o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut );
733738 motherRot = ROOT::Math::PxPyPzMVector (mother.Px () * std::cos (theta2) - mother.Py () * std::sin (theta2), mother.Px () * std::sin (theta2) + mother.Py () * std::cos (theta2), mother.Pz (), mother.M ());
734- if (std::abs (motherRot.Rapidity ()) < 0.5 ) {
739+ if (std::abs (motherRot.Rapidity ()) < config. rapidityMotherData ) {
735740 hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, motherRot.Pt (), motherRot.M (), cosThetaStarRandom, anglePhi);
736741 }
737742 }
738743 } else {
739- if (std::abs (mother.Rapidity ()) < 0.5 ) {
744+ if (std::abs (mother.Rapidity ()) < config. rapidityMotherData ) {
740745 hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, mother.Pt (), mother.M (), cosThetaStarRandom, anglePhi);
741746 }
742747 }
@@ -854,7 +859,7 @@ struct HigherMassResonances {
854859 }
855860 int sizeofv0indexes = v0indexes.size ();
856861 rKzeroShort.fill (HIST (" NksProduced" ), sizeofv0indexes);
857- if (config.selectTWOKsOnly && sizeofv0indexes == 2 && allConditionsMet) {
862+ if (config.selectTWOKsOnly && sizeofv0indexes == config. numbers [ 2 ] && allConditionsMet) {
858863 fillInvMass (mother, multiplicity, daughter1, daughter2, false );
859864 }
860865 v0indexes.clear ();
@@ -975,7 +980,7 @@ struct HigherMassResonances {
975980 }
976981 int sizeofv0indexes = v0indexes.size ();
977982 rKzeroShort.fill (HIST (" NksProduced" ), sizeofv0indexes);
978- if (config.selectTWOKsOnly && sizeofv0indexes == 2 && allConditionsMet) {
983+ if (config.selectTWOKsOnly && sizeofv0indexes == config. numbers [ 2 ] && allConditionsMet) {
979984 fillInvMass (mother, multiplicity, daughter1, daughter2, false );
980985 }
981986 v0indexes.clear ();
@@ -1348,7 +1353,7 @@ struct HigherMassResonances {
13481353 }
13491354 hMChists.fill (HIST (" events_check" ), 5.5 );
13501355
1351- if (config.applyRapidityMC && std::abs (mcParticle.y ()) >= 0.5 ) {
1356+ if (config.applyRapidityMC && std::abs (mcParticle.y ()) >= config. rapidityMotherData ) {
13521357 continue ;
13531358 }
13541359 hMChists.fill (HIST (" events_check" ), 6.5 );
@@ -1358,7 +1363,7 @@ struct HigherMassResonances {
13581363 // counter++;
13591364
13601365 auto kDaughters = mcParticle.daughters_as <aod::McParticles>();
1361- if (kDaughters .size () != 2 ) {
1366+ if (kDaughters .size () != config. numbers [ 2 ] ) {
13621367 continue ;
13631368 }
13641369 hMChists.fill (HIST (" events_check" ), 7.5 );
@@ -1370,17 +1375,17 @@ struct HigherMassResonances {
13701375 continue ;
13711376 }
13721377 hMChists.fill (HIST (" events_check" ), 8.5 );
1373- if (std::abs (kCurrentDaughter .pdgCode ()) == 310 ) {
1378+ if (std::abs (kCurrentDaughter .pdgCode ()) == PDG_t:: kK0Short ) {
13741379 passKs.push_back (true );
13751380 hMChists.fill (HIST (" events_check" ), 9.5 );
13761381 if (passKs.size () == 1 ) {
13771382 daughter1 = ROOT::Math::PxPyPzMVector (kCurrentDaughter .px (), kCurrentDaughter .py (), kCurrentDaughter .pz (), o2::constants::physics::MassK0Short);
1378- } else if (passKs.size () == 2 ) {
1383+ } else if (passKs.size () == config. numbers [ 2 ] ) {
13791384 daughter2 = ROOT::Math::PxPyPzMVector (kCurrentDaughter .px (), kCurrentDaughter .py (), kCurrentDaughter .pz (), o2::constants::physics::MassK0Short);
13801385 }
13811386 }
13821387 }
1383- if (passKs.size () == 2 ) {
1388+ if (passKs.size () == config. numbers [ 2 ] ) {
13841389 lResonanceGen = ROOT::Math::PxPyPzEVector (mcParticle.pt (), mcParticle.eta (), mcParticle.phi (), mcParticle.e ());
13851390 lResonanceGen1 = daughter1 + daughter2;
13861391
@@ -1399,7 +1404,7 @@ struct HigherMassResonances {
13991404 hMChists.fill (HIST (" GenEta" ), mcParticle.eta ());
14001405 hMChists.fill (HIST (" GenPhi" ), mcParticle.phi ());
14011406
1402- if (config.applyPairRapidityGen && std::abs (lResonanceGen1.Rapidity ()) >= 0.5 ) {
1407+ if (config.applyPairRapidityGen && std::abs (lResonanceGen1.Rapidity ()) >= config. rapidityMotherData ) {
14031408 continue ;
14041409 }
14051410
@@ -1507,7 +1512,7 @@ struct HigherMassResonances {
15071512 int trackv0PDG1 = std::abs (mctrackv01.pdgCode ());
15081513 int trackv0PDG2 = std::abs (mctrackv02.pdgCode ());
15091514
1510- if (std::abs (trackv0PDG1) != 310 || std::abs (trackv0PDG2) != 310 ) {
1515+ if (std::abs (trackv0PDG1) != PDG_t:: kK0Short || std::abs (trackv0PDG2) != PDG_t:: kK0Short ) {
15111516 continue ;
15121517 }
15131518 hMChists.fill (HIST (" events_checkrec" ), 12.5 );
@@ -1554,7 +1559,7 @@ struct HigherMassResonances {
15541559 }
15551560 hMChists.fill (HIST (" events_checkrec" ), 18.5 );
15561561
1557- if (config.applyRapidityMC && std::abs (mothertrack1.y ()) >= 0.5 ) {
1562+ if (config.applyRapidityMC && std::abs (mothertrack1.y ()) >= config. rapidityMotherData ) {
15581563 continue ;
15591564 }
15601565 hMChists.fill (HIST (" events_checkrec" ), 19.5 );
@@ -1586,7 +1591,7 @@ struct HigherMassResonances {
15861591 hMChists.fill (HIST (" RecPhi" ), mothertrack1.phi ());
15871592 hMChists.fill (HIST (" RecEta" ), mothertrack1.eta ());
15881593
1589- if (config.applyPairRapidityRec && std::abs (mother.Rapidity ()) >= 0.5 ) {
1594+ if (config.applyPairRapidityRec && std::abs (mother.Rapidity ()) >= config. rapidityMotherData ) {
15901595 continue ;
15911596 }
15921597
0 commit comments