Skip to content

Commit 9ec6f95

Browse files
author
Sawan Sawan
committed
rapidity correction
1 parent 7c27e70 commit 9ec6f95

File tree

1 file changed

+55
-56
lines changed

1 file changed

+55
-56
lines changed

‎PWGLF/Tasks/Resonances/higherMassResonances.cxx‎

Lines changed: 55 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,6 @@
4343
#include <TH1F.h>
4444
#include <TH2F.h>
4545
#include <THn.h>
46-
#include <TLorentzVector.h>
4746
#include <TMath.h>
4847
#include <TObjArray.h>
4948
#include <TPDGCode.h>
@@ -143,7 +142,7 @@ struct HigherMassResonances {
143142
Configurable<bool> cTVXEvsel{"cTVXEvsel", true, "Triggger selection"};
144143
Configurable<bool> avoidsplitrackMC{"avoidsplitrackMC", false, "avoid split track in MC"};
145144
Configurable<int> selectMCparticles{"selectMCparticles", 1, "0: f0(1710), 1: f2(1525), 2: a2(1320), 3: f0(1370), 4: f0(1500)"};
146-
Configurable<bool> apply_rapidityMC{"apply_rapidityMC", true, "Apply rapidity cut on generated and reconstructed particles"};
145+
Configurable<bool> applyRapidityMC{"applyRapidityMC", true, "Apply rapidity cut on generated and reconstructed particles"};
147146
std::vector<int> pdgCodes = {10331, 335, 115, 10221, 9030221};
148147

149148
// output THnSparses
@@ -182,13 +181,13 @@ struct HigherMassResonances {
182181
ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, motherRot, fourVecDauCM, fourVecDauCM1;
183182
ROOT::Math::PxPyPzEVector mother1;
184183
ROOT::Math::XYZVector randomVec, beamVec, normalVec;
185-
ROOT::Math::XYZVectorF v1_CM, zaxis_HE, yaxis_HE, xaxis_HE;
184+
ROOT::Math::XYZVectorF v1CM, zaxisHE, yaxisHE, xaxisHE;
186185
// ROOT::Math::XYZVector threeVecDauCM, helicityVec, randomVec, beamVec, normalVec;
187186
ROOT::Math::XYZVector zBeam; // ẑ: beam direction in lab frame
188-
double BeamMomentum = std::sqrt(13600 * 13600 / 4 - 0.938 * 0.938); // GeV
189-
ROOT::Math::PxPyPzEVector Beam1{0., 0., -BeamMomentum, 13600. / 2.};
190-
ROOT::Math::PxPyPzEVector Beam2{0., 0., BeamMomentum, 13600. / 2.};
191-
ROOT::Math::XYZVectorF Beam1_CM, Beam2_CM;
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.};
190+
ROOT::Math::XYZVectorF beam1CM, beam2CM;
192191

193192
// const double massK0s = o2::constants::physics::MassK0Short;
194193
bool isMix = false;
@@ -615,8 +614,8 @@ struct HigherMassResonances {
615614
fourVecDauCM = boost(daughter1); // boost the frame of daughter to the center of mass frame
616615
// threeVecDauCM = fourVecDauCM.Vect(); // get the 3 vector of daughter in the frame of mother
617616

618-
Beam1_CM = ROOT::Math::XYZVectorF((boost(Beam1).Vect()).Unit());
619-
Beam2_CM = ROOT::Math::XYZVectorF((boost(Beam2).Vect()).Unit());
617+
beam1CM = ROOT::Math::XYZVectorF((boost(beam1).Vect()).Unit());
618+
beam2CM = ROOT::Math::XYZVectorF((boost(beam2).Vect()).Unit());
620619

621620
// define y = zBeam x z: Normal to the production plane
622621
// ẑ: mother direction in lab, boosted into mother's rest frame
@@ -634,22 +633,22 @@ struct HigherMassResonances {
634633
// auto p_proj_y = threeVecDauCM.Dot(y_axis);
635634

636635
// // Calculate φ in [-π, π]
637-
// auto angle_phi = std::atan2(p_proj_y, p_proj_x); // φ in radians
636+
// auto anglePhi = std::atan2(p_proj_y, p_proj_x); // φ in radians
638637

639-
v1_CM = ROOT::Math::XYZVectorF(boost(daughter1).Vect()).Unit();
638+
v1CM = ROOT::Math::XYZVectorF(boost(daughter1).Vect()).Unit();
640639
// ROOT::Math::XYZVectorF v2_CM{(boost(daughter1).Vect()).Unit()};
641640
// using positive sign convention for the first track
642-
// ROOT::Math::XYZVectorF v_CM = (t1.sign() > 0 ? v1_CM : v2_CM); // here selected decay daughter momentum is intested. here you can choose one decay daughter no need to check both case as it is neutral particle for our case
641+
// ROOT::Math::XYZVectorF v_CM = (t1.sign() > 0 ? v1CM : v2_CM); // here selected decay daughter momentum is intested. here you can choose one decay daughter no need to check both case as it is neutral particle for our case
643642
// Helicity frame
644-
zaxis_HE = ROOT::Math::XYZVectorF(mother.Vect()).Unit();
645-
yaxis_HE = ROOT::Math::XYZVectorF(Beam1_CM.Cross(Beam2_CM)).Unit();
646-
xaxis_HE = ROOT::Math::XYZVectorF(yaxis_HE.Cross(zaxis_HE)).Unit();
643+
zaxisHE = ROOT::Math::XYZVectorF(mother.Vect()).Unit();
644+
yaxisHE = ROOT::Math::XYZVectorF(beam1CM.Cross(beam2CM)).Unit();
645+
xaxisHE = ROOT::Math::XYZVectorF(yaxisHE.Cross(zaxisHE)).Unit();
647646

648-
// CosThetaHE = zaxis_HE.Dot(v_CM);
647+
// CosThetaHE = zaxisHE.Dot(v_CM);
649648

650-
auto angle_phi = TMath::ATan2(yaxis_HE.Dot(v1_CM), xaxis_HE.Dot(v1_CM));
651-
if (angle_phi < 0) {
652-
angle_phi += 2 * TMath::Pi(); // ensure phi is in [0, 2pi]
649+
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]
653652
}
654653

655654
// if (std::abs(mother.Rapidity()) < 0.5) {
@@ -659,7 +658,7 @@ struct HigherMassResonances {
659658
auto cosThetaStarHelicity = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2()));
660659
if (!isMix) {
661660
if (std::abs(mother.Rapidity()) < 0.5) {
662-
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, angle_phi);
661+
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi);
663662
}
664663

665664
for (int i = 0; i < config.cRotations; i++) {
@@ -674,49 +673,49 @@ struct HigherMassResonances {
674673

675674
auto cosThetaStarHelicityRot = motherRot.Vect().Dot(daughterRotCM.Vect()) / (std::sqrt(daughterRotCM.Vect().Mag2()) * std::sqrt(motherRot.Vect().Mag2()));
676675
if (motherRot.Rapidity() < 0.5)
677-
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, angle_phi);
676+
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, anglePhi);
678677
}
679678
} else {
680679
if (std::abs(mother.Rapidity()) < 0.5) {
681-
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, angle_phi);
680+
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi);
682681
}
683682
}
684683
} else if (config.activateTHnSparseCosThStarProduction) {
685684
normalVec = ROOT::Math::XYZVector(mother.Py(), -mother.Px(), 0.f);
686685
auto cosThetaStarProduction = normalVec.Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(normalVec.Mag2()));
687686
if (!isMix) {
688687
if (std::abs(mother.Rapidity()) < 0.5) {
689-
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, angle_phi);
688+
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi);
690689
}
691690
for (int i = 0; i < config.cRotations; i++) {
692691
theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
693692
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());
694693
if (std::abs(motherRot.Rapidity()) < 0.5) {
695-
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarProduction, angle_phi);
694+
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarProduction, anglePhi);
696695
}
697696
}
698697
} else {
699698
if (std::abs(mother.Rapidity()) < 0.5) {
700-
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, angle_phi);
699+
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi);
701700
}
702701
}
703702
} else if (config.activateTHnSparseCosThStarBeam) {
704703
beamVec = ROOT::Math::XYZVector(0.f, 0.f, 1.f);
705704
auto cosThetaStarBeam = beamVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2());
706705
if (!isMix) {
707706
if (std::abs(mother.Rapidity()) < 0.5) {
708-
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, angle_phi);
707+
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi);
709708
}
710709
for (int i = 0; i < config.cRotations; i++) {
711710
theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
712711
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());
713712
if (std::abs(motherRot.Rapidity()) < 0.5) {
714-
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam, angle_phi);
713+
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam, anglePhi);
715714
}
716715
}
717716
} else {
718717
if (std::abs(mother.Rapidity()) < 0.5) {
719-
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, angle_phi);
718+
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi);
720719
}
721720
}
722721
} else if (config.activateTHnSparseCosThStarRandom) {
@@ -727,18 +726,18 @@ struct HigherMassResonances {
727726
auto cosThetaStarRandom = randomVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2());
728727
if (!isMix) {
729728
if (std::abs(mother.Rapidity()) < 0.5) {
730-
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, angle_phi);
729+
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi);
731730
}
732731
for (int i = 0; i < config.cRotations; i++) {
733732
theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
734733
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());
735734
if (std::abs(motherRot.Rapidity()) < 0.5) {
736-
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, angle_phi);
735+
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, anglePhi);
737736
}
738737
}
739738
} else {
740739
if (std::abs(mother.Rapidity()) < 0.5) {
741-
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, angle_phi);
740+
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi);
742741
}
743742
}
744743
}
@@ -864,7 +863,7 @@ struct HigherMassResonances {
864863

865864
using EventCandidatesDerivedData = soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels, aod::StraStamps, aod::StraZDCSP>;
866865
using V0CandidatesDerivedData = soa::Join<aod::V0CollRefs, aod::V0Cores, aod::V0Extras, aod::V0TOFPIDs, aod::V0TOFNSigmas>;
867-
using dauTracks = soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>;
866+
using DauTracks = soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>;
868867

869868
void processSEderived(EventCandidatesDerivedData::iterator const& collision, TrackCandidates const& /*tracks*/, aod::V0Datas const& V0s)
870869
{
@@ -1290,8 +1289,8 @@ struct HigherMassResonances {
12901289
int counter = 0;
12911290
float multiplicityGen = 0.0;
12921291
std::vector<bool> passKs;
1293-
ROOT::Math::PxPyPzMVector lResonance_gen1;
1294-
ROOT::Math::PxPyPzEVector lResonance_gen;
1292+
ROOT::Math::PxPyPzMVector lResonanceGen1;
1293+
ROOT::Math::PxPyPzEVector lResonanceGen;
12951294

12961295
void processGen(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups<EventCandidatesMC>& collisions)
12971296
{
@@ -1349,7 +1348,7 @@ struct HigherMassResonances {
13491348
}
13501349
hMChists.fill(HIST("events_check"), 5.5);
13511350

1352-
if (config.apply_rapidityMC && std::abs(mcParticle.y()) >= 0.5) {
1351+
if (config.applyRapidityMC && std::abs(mcParticle.y()) >= 0.5) {
13531352
continue;
13541353
}
13551354
hMChists.fill(HIST("events_check"), 6.5);
@@ -1382,33 +1381,33 @@ struct HigherMassResonances {
13821381
}
13831382
}
13841383
if (passKs.size() == 2) {
1385-
lResonance_gen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e());
1386-
lResonance_gen1 = daughter1 + daughter2;
1384+
lResonanceGen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e());
1385+
lResonanceGen1 = daughter1 + daughter2;
13871386

1388-
ROOT::Math::Boost boost{lResonance_gen.BoostToCM()};
1389-
ROOT::Math::Boost boost1{lResonance_gen1.BoostToCM()};
1387+
ROOT::Math::Boost boost{lResonanceGen.BoostToCM()};
1388+
ROOT::Math::Boost boost1{lResonanceGen1.BoostToCM()};
13901389

13911390
fourVecDauCM = boost(daughter1);
13921391
fourVecDauCM1 = boost1(daughter1);
13931392

1394-
auto helicity_gen = lResonance_gen.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonance_gen.Vect().Mag2()));
1395-
auto helicity_gen1 = lResonance_gen1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(lResonance_gen1.Vect().Mag2()));
1393+
auto helicityGen = lResonanceGen.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonanceGen.Vect().Mag2()));
1394+
auto helicityGen1 = lResonanceGen1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(lResonanceGen1.Vect().Mag2()));
13961395

1397-
hMChists.fill(HIST("Genf1710"), multiplicityGen, lResonance_gen.pt(), helicity_gen);
1398-
hMChists.fill(HIST("Genf1710_mass"), lResonance_gen.M());
1396+
hMChists.fill(HIST("Genf1710"), multiplicityGen, lResonanceGen.pt(), helicityGen);
1397+
hMChists.fill(HIST("Genf1710_mass"), lResonanceGen.M());
13991398
hMChists.fill(HIST("GenRapidity"), mcParticle.y());
14001399
hMChists.fill(HIST("GenEta"), mcParticle.eta());
14011400
hMChists.fill(HIST("GenPhi"), mcParticle.phi());
14021401

1403-
if (config.applyPairRapidityGen && std::abs(lResonance_gen1.Y()) >= 0.5) {
1402+
if (config.applyPairRapidityGen && std::abs(lResonanceGen1.Rapidity()) >= 0.5) {
14041403
continue;
14051404
}
14061405

1407-
hMChists.fill(HIST("Genf17102"), multiplicityGen, lResonance_gen1.pt(), helicity_gen1);
1408-
hMChists.fill(HIST("Genf1710_mass2"), lResonance_gen1.M());
1409-
hMChists.fill(HIST("GenRapidity2"), lResonance_gen1.Y());
1410-
hMChists.fill(HIST("GenEta2"), lResonance_gen1.Eta());
1411-
hMChists.fill(HIST("GenPhi2"), lResonance_gen1.Phi());
1406+
hMChists.fill(HIST("Genf17102"), multiplicityGen, lResonanceGen1.pt(), helicityGen1);
1407+
hMChists.fill(HIST("Genf1710_mass2"), lResonanceGen1.M());
1408+
hMChists.fill(HIST("GenRapidity2"), lResonanceGen1.Rapidity());
1409+
hMChists.fill(HIST("GenEta2"), lResonanceGen1.Eta());
1410+
hMChists.fill(HIST("GenPhi2"), lResonanceGen1.Phi());
14121411
}
14131412
passKs.clear(); // clear the vector for the next iteration
14141413
}
@@ -1555,7 +1554,7 @@ struct HigherMassResonances {
15551554
}
15561555
hMChists.fill(HIST("events_checkrec"), 18.5);
15571556

1558-
if (config.apply_rapidityMC && std::abs(mothertrack1.y()) >= 0.5) {
1557+
if (config.applyRapidityMC && std::abs(mothertrack1.y()) >= 0.5) {
15591558
continue;
15601559
}
15611560
hMChists.fill(HIST("events_checkrec"), 19.5);
@@ -1578,21 +1577,21 @@ struct HigherMassResonances {
15781577
fourVecDauCM = boost(daughter1);
15791578
fourVecDauCM1 = boost1(daughter1);
15801579

1581-
auto helicity_rec = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2()));
1580+
auto helicityRec = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2()));
15821581

1583-
auto helicity_rec2 = mother1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(mother1.Vect().Mag2()));
1582+
auto helicityRec2 = mother1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(mother1.Vect().Mag2()));
15841583

1585-
hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mothertrack1.pt(), mother1.M(), helicity_rec2);
1584+
hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mothertrack1.pt(), mother1.M(), helicityRec2);
15861585
hMChists.fill(HIST("RecRapidity"), mothertrack1.y());
15871586
hMChists.fill(HIST("RecPhi"), mothertrack1.phi());
15881587
hMChists.fill(HIST("RecEta"), mothertrack1.eta());
15891588

1590-
if (config.applyPairRapidityRec && std::abs(mother.Y()) >= 0.5) {
1589+
if (config.applyPairRapidityRec && std::abs(mother.Rapidity()) >= 0.5) {
15911590
continue;
15921591
}
15931592

1594-
hMChists.fill(HIST("Recf1710_pt2"), multiplicity, mother.Pt(), mother.M(), helicity_rec);
1595-
hMChists.fill(HIST("RecRapidity2"), mother.Y());
1593+
hMChists.fill(HIST("Recf1710_pt2"), multiplicity, mother.Pt(), mother.M(), helicityRec);
1594+
hMChists.fill(HIST("RecRapidity2"), mother.Rapidity());
15961595
hMChists.fill(HIST("RecPhi2"), mother.Phi());
15971596
hMChists.fill(HIST("RecEta2"), mother.Eta());
15981597
}

0 commit comments

Comments
 (0)