Skip to content

Commit 0438000

Browse files
sawankumawatSawan Sawan
sawankumawat
and
Sawan Sawan
authored
[PWGLF] "PWGLF: modified Reconstructed and Generated MC" (#11931)
Co-authored-by: Sawan Sawan <sawan.sawan@cern.ch>
1 parent 258a047 commit 0438000

File tree

1 file changed

+63
-87
lines changed

1 file changed

+63
-87
lines changed

‎PWGLF/Tasks/Resonances/higherMassResonances.cxx‎

Lines changed: 63 additions & 87 deletions
Original file line numberDiff line numberDiff line change
@@ -176,7 +176,7 @@ struct HigherMassResonances {
176176
// variables declaration
177177
float multiplicity = 0.0f;
178178
float theta2;
179-
ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, motherRot, fourVecDauCM;
179+
ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, mother1, motherRot, fourVecDauCM, fourVecDauCM1;
180180
ROOT::Math::XYZVector randomVec, beamVec, normalVec;
181181
ROOT::Math::XYZVectorF v1_CM, zaxis_HE, yaxis_HE, xaxis_HE;
182182
// ROOT::Math::XYZVector threeVecDauCM, helicityVec, randomVec, beamVec, normalVec;
@@ -309,18 +309,17 @@ struct HigherMassResonances {
309309
// For MC
310310
if (config.isMC) {
311311
hMChists.add("events_check", "No. of events in the generated MC", kTH1I, {{20, 0, 20}});
312-
hMChists.add("events_checkrec", "No. of events in the reconstructed MC", kTH1I, {{20, 0, 20}});
312+
hMChists.add("events_checkrec", "No. of events in the reconstructed MC", kTH1I, {{25, 0, 25}});
313313
hMChists.add("Genf1710", "Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL});
314314
hMChists.add("Recf1710_pt1", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
315-
hMChists.add("Recf1710_pt2", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
315+
hMChists.add("Recf1710_ptTemp", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
316316
// hMChists.add("Recf1710_p", "Rec f_{0}(1710) p", kTH1F, {ptAxis});
317317
hMChists.add("h1Recsplit", "Rec p_{T}2", kTH1F, {ptAxis});
318318
// hMChists.add("Recf1710_mass", "Rec f_{0}(1710) mass", kTH1F, {glueballMassAxis});
319319
hMChists.add("Genf1710_mass", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
320-
hMChists.add("Genf1710_mass2", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
321-
hMChists.add("Genf1710_pt2", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
320+
hMChists.add("Genf1710_pt2", "Gen f_{0}(1710) p_{T}", kTH1F, {ptAxis});
322321
hMChists.add("GenEta", "Gen Eta", kTHnSparseF, {ptAxis, {100, -1.0f, 1.0f}});
323-
hMChists.add("GenPhi", "Gen Phi", kTH1F, {{70, -3.5f, 3.5f}});
322+
hMChists.add("GenPhi", "Gen Phi", kTH1F, {{70, 0.0, 7.0f}});
324323
hMChists.add("GenRapidity", "Gen Rapidity", kTHnSparseF, {ptAxis, {100, -1.0f, 1.0f}});
325324
hMChists.add("RecEta", "Rec Eta", kTH1F, {{100, -1.0f, 1.0f}});
326325
hMChists.add("RecPhi", "Rec Phi", kTH1F, {{70, 0.0f, 7.0f}});
@@ -590,7 +589,7 @@ struct HigherMassResonances {
590589
// For Monte Carlo
591590
using EventCandidatesMC = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Cs, aod::CentFT0Ms>;
592591
using TrackCandidatesMC = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullKa, aod::pidTOFFullKa, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::McTrackLabels>>;
593-
using V0TrackCandidatesMC = soa::Join<aod::V0Datas, aod::McV0Labels>;
592+
using V0TrackCandidatesMC = soa::Filtered<soa::Join<aod::V0Datas, aod::McV0Labels>>;
594593
// zBeam direction in lab frame
595594

596595
template <typename T>
@@ -1052,13 +1051,6 @@ struct HigherMassResonances {
10521051
// std::cout << "px " << mcParticle.px() << " py " << mcParticle.py() << " pz " << mcParticle.pz() << " y " << mcParticle.y() << std::endl;
10531052
// counter++;
10541053

1055-
hMChists.fill(HIST("GenRapidity"), mcParticle.pt(), mcParticle.y());
1056-
hMChists.fill(HIST("GenPhi"), mcParticle.phi());
1057-
hMChists.fill(HIST("GenEta"), mcParticle.pt(), mcParticle.eta());
1058-
// hMChists.fill(HIST("GenPx"), mcParticle.px());
1059-
// hMChists.fill(HIST("GenPy"), mcParticle.py());
1060-
// hMChists.fill(HIST("GenPz"), mcParticle.pz());
1061-
10621054
auto kDaughters = mcParticle.daughters_as<aod::McParticles>();
10631055
if (kDaughters.size() != 2) {
10641056
continue;
@@ -1083,25 +1075,26 @@ struct HigherMassResonances {
10831075
}
10841076
}
10851077
if (passKs.size() == 2) {
1086-
lResonance_gen = ROOT::Math::PxPyPzMVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e());
1087-
lResonance_gen2 = daughter1 + daughter2; // invariant mass of Kshort pair
1078+
lResonance_gen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e());
1079+
lResonance_gen2 = daughter1 + daughter2;
10881080

10891081
ROOT::Math::Boost boost{lResonance_gen.BoostToCM()};
10901082
fourVecDauCM = boost(daughter1); // boost the frame of daughter to the center of mass frame
10911083

1092-
auto helicity_gen = lResonance_gen.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonance_gen.Vect().Mag2()));
1084+
auto helicity_gen = lResonance_gen2.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonance_gen2.Vect().Mag2()));
10931085

1094-
hMChists.fill(HIST("Genf1710"), multiplicityGen, mcParticle.pt(), helicity_gen);
1095-
hMChists.fill(HIST("Genf1710_mass"), lResonance_gen.M());
1096-
hMChists.fill(HIST("Genf1710_mass2"), lResonance_gen2.M());
1097-
hMChists.fill(HIST("Genf1710_pt2"), lResonance_gen2.Pt());
1086+
hMChists.fill(HIST("Genf1710"), multiplicityGen, lResonance_gen2.pt(), helicity_gen);
1087+
hMChists.fill(HIST("Genf1710_mass"), lResonance_gen2.M());
1088+
hMChists.fill(HIST("Genf1710_pt2"), mcParticle.pt());
1089+
hMChists.fill(HIST("GenRapidity"), lResonance_gen2.Pt(), lResonance_gen2.Y());
1090+
hMChists.fill(HIST("GenPhi"), lResonance_gen2.Phi());
1091+
hMChists.fill(HIST("GenEta"), lResonance_gen2.Pt(), lResonance_gen2.Eta());
10981092
}
10991093
passKs.clear(); // clear the vector for the next iteration
11001094
}
11011095
}
11021096
PROCESS_SWITCH(HigherMassResonances, processGen, "Process Generated", false);
11031097

1104-
int counter2 = 0;
11051098
int eventCounter = 0;
11061099
std::vector<int> gindex1, gindex2;
11071100
void processRec(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const&, V0TrackCandidatesMC const& V0s, aod::McParticles const&, aod::McCollisions const& /*mcCollisions*/)
@@ -1110,7 +1103,6 @@ struct HigherMassResonances {
11101103
return;
11111104
}
11121105

1113-
ROOT::Math::PxPyPzMVector lDecayDaughter1, lDecayDaughter2, lResonance;
11141106
auto multiplicity = collision.centFT0C();
11151107
hMChists.fill(HIST("Rec_Multiplicity"), multiplicity);
11161108

@@ -1125,17 +1117,21 @@ struct HigherMassResonances {
11251117
}
11261118
hMChists.fill(HIST("events_checkrec"), 2.5);
11271119

1128-
if (config.timFrameEvsel && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
1129-
return;
1130-
}
1131-
hMChists.fill(HIST("events_checkrec"), 3.5);
1132-
if (config.cTVXEvsel && (!collision.selection_bit(aod::evsel::kIsTriggerTVX))) {
1120+
// if (config.timFrameEvsel && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
1121+
// return;
1122+
// }
1123+
// hMChists.fill(HIST("events_checkrec"), 3.5);
1124+
// if (config.cTVXEvsel && (!collision.selection_bit(aod::evsel::kIsTriggerTVX))) {
1125+
// return;
1126+
// }
1127+
1128+
if (!collision.sel8()) {
11331129
return;
11341130
}
11351131
hMChists.fill(HIST("events_checkrec"), 4.5);
11361132
hMChists.fill(HIST("MC_mult_after_event_sel"), multiplicity);
11371133
eventCounter++;
1138-
auto oldindex = -999;
1134+
// auto oldindex = -999;
11391135

11401136
for (const auto& v01 : V0s) {
11411137

@@ -1168,7 +1164,6 @@ struct HigherMassResonances {
11681164

11691165
double nTPCSigmaPos1[1]{postrack1.tpcNSigmaPi()};
11701166
double nTPCSigmaNeg1[1]{negtrack1.tpcNSigmaPi()};
1171-
11721167
double nTPCSigmaPos2[1]{postrack2.tpcNSigmaPi()};
11731168
double nTPCSigmaNeg2[1]{negtrack2.tpcNSigmaPi()};
11741169

@@ -1207,92 +1202,73 @@ struct HigherMassResonances {
12071202
continue;
12081203
}
12091204
}
1210-
// if (counter2 < 1e4)
1211-
// std::cout << "Mother1 pdg code: " << motpdgs << " p_{T} " << mothertrack1.pt() << "Global index " << mothertrack1.globalIndex() << " event " << eventCounter << std::endl;
1212-
// counter2++;
1213-
1214-
// int counter_check = 0;
12151205

12161206
for (const auto& mothertrack2 : mctrackv02.mothers_as<aod::McParticles>()) {
12171207

12181208
hMChists.fill(HIST("events_checkrec"), 13.5);
12191209

1220-
if (mothertrack1.pdgCode() != mothertrack2.pdgCode()) {
1210+
if (mothertrack1.pdgCode() != config.pdgCodes[config.selectMCparticles]) {
12211211
continue;
12221212
}
12231213
hMChists.fill(HIST("events_checkrec"), 14.5);
12241214

1225-
// int motpdgs2 = std::abs(mothertrack2.pdgCode());
1215+
if (mothertrack1.pdgCode() != mothertrack2.pdgCode()) {
1216+
continue;
1217+
}
1218+
hMChists.fill(HIST("events_checkrec"), 15.5);
1219+
12261220
gindex2.push_back(mothertrack2.globalIndex());
12271221
if (gindex2.size() > 1) {
12281222
if (std::find(gindex2.begin(), gindex2.end(), mothertrack2.globalIndex()) != gindex2.end()) {
12291223
continue;
12301224
}
12311225
}
1232-
// if (counter2 < 1e4)
1233-
// std::cout << "Mother2 pdg code: " << motpdgs2 << " p_{T} " << mothertrack2.pt() << "Global index " << mothertrack1.globalIndex() << " event " << eventCounter << std::endl;
1234-
1235-
if (mothertrack1.pdgCode() != config.pdgCodes[config.selectMCparticles]) {
1236-
continue;
1237-
}
1238-
hMChists.fill(HIST("events_checkrec"), 15.5);
1239-
1240-
if (mothertrack1.globalIndex() != mothertrack2.globalIndex()) {
1241-
continue;
1242-
}
12431226
hMChists.fill(HIST("events_checkrec"), 16.5);
12441227

1245-
if (!mothertrack1.producedByGenerator()) {
1228+
if (mothertrack1.globalIndex() != mothertrack2.globalIndex()) {
12461229
continue;
12471230
}
12481231
hMChists.fill(HIST("events_checkrec"), 17.5);
12491232

1250-
if (config.apply_rapidityMC && std::abs(mothertrack1.y()) >= 0.5) {
1233+
if (!mothertrack1.producedByGenerator()) {
12511234
continue;
12521235
}
12531236
hMChists.fill(HIST("events_checkrec"), 18.5);
12541237

1255-
if (config.avoidsplitrackMC && oldindex == mothertrack1.globalIndex()) {
1256-
hMChists.fill(HIST("h1Recsplit"), mothertrack1.pt());
1238+
if (config.apply_rapidityMC && std::abs(mothertrack1.y()) >= 0.5) {
12571239
continue;
12581240
}
1259-
oldindex = mothertrack1.globalIndex();
1241+
hMChists.fill(HIST("events_checkrec"), 19.5);
12601242

1261-
// counter_check++;
1262-
// if (counter_check > 1) {
1263-
// std::cout << "Total mothers is " << counter_check << std::endl;
1243+
// if (config.avoidsplitrackMC && oldindex == mothertrack1.globalIndex()) {
1244+
// hMChists.fill(HIST("h1Recsplit"), mothertrack1.pt());
1245+
// continue;
12641246
// }
1265-
// std::cout << "After selection " << " p_{T} " << mothertrack2.pt() << " event " << eventCounter << std::endl;
1266-
1267-
pvec0 = std::array{v01.px(), v01.py(), v01.pz()};
1268-
pvec1 = std::array{v02.px(), v02.py(), v02.pz()};
1269-
auto arrMomrec = std::array{pvec0, pvec1};
1270-
// auto motherP = mothertrack1.p();
1271-
// auto motherE = mothertrack1.e();
1272-
// auto genMass = std::sqrt(motherE * motherE - motherP * motherP);
1273-
auto recMass = RecoDecay::m(arrMomrec, std::array{o2::constants::physics::MassK0Short, o2::constants::physics::MassK0Short});
1274-
// auto recpt = TMath::Sqrt((track1.px() + track2.px()) * (track1.px() + track2.px()) + (track1.py() + track2.py()) * (track1.py() + track2.py()));
1275-
//// Resonance reconstruction
1276-
lDecayDaughter1 = ROOT::Math::PxPyPzMVector(v01.px(), v01.py(), v01.pz(), o2::constants::physics::MassK0Short);
1277-
lDecayDaughter2 = ROOT::Math::PxPyPzMVector(v02.px(), v02.py(), v02.pz(), o2::constants::physics::MassK0Short);
1278-
lResonance = lDecayDaughter1 + lDecayDaughter2;
1279-
if (config.apply_rapidityMC && std::abs(lResonance.Y()) >= 0.5) {
1280-
continue;
1281-
}
1282-
// daughter1, mother, fourVecDauCM
1283-
ROOT::Math::Boost boost{lResonance.BoostToCM()};
1284-
fourVecDauCM = boost(lDecayDaughter1); // boost the frame of daughter to the center of mass frame
1285-
auto helicity_rec = lResonance.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonance.Vect().Mag2()));
1286-
1287-
// hMChists.fill(HIST("Recf1710_p"), motherP);
1288-
// hMChists.fill(HIST("Recf1710_mass"), recMass);
1289-
hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mothertrack1.pt(), recMass, helicity_rec);
1290-
// hMChists.fill(HIST("Genf1710_mass"), genMass);
1291-
hMChists.fill(HIST("Recf1710_pt2"), multiplicity, lResonance.Pt(), recMass, helicity_rec);
1292-
1293-
hMChists.fill(HIST("RecRapidity"), mothertrack1.y());
1294-
hMChists.fill(HIST("RecPhi"), mothertrack1.phi());
1295-
hMChists.fill(HIST("RecEta"), mothertrack1.eta());
1247+
// hMChists.fill(HIST("events_checkrec"), 20.5);
1248+
// oldindex = mothertrack1.globalIndex(); // split tracks is already handled using gindex1 and gindex2
1249+
1250+
daughter1 = ROOT::Math::PxPyPzMVector(v01.px(), v01.py(), v01.pz(), o2::constants::physics::MassK0Short);
1251+
daughter2 = ROOT::Math::PxPyPzMVector(v02.px(), v02.py(), v02.pz(), o2::constants::physics::MassK0Short);
1252+
mother = daughter1 + daughter2;
1253+
mother1 = ROOT::Math::PxPyPzEVector(mothertrack1.px(), mothertrack1.py(), mothertrack1.pz(), mothertrack1.e());
1254+
1255+
ROOT::Math::Boost boost{mother.BoostToCM()};
1256+
ROOT::Math::Boost boost1{mother1.BoostToCM()};
1257+
1258+
fourVecDauCM = boost(daughter1);
1259+
fourVecDauCM1 = boost1(daughter1);
1260+
1261+
auto helicity_rec = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2()));
1262+
1263+
auto helicity_rec2 = mother1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(mother1.Vect().Mag2()));
1264+
1265+
// std::cout << "mother pT is " << mother.Pt() << std::endl;
1266+
1267+
hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mother.Pt(), mother.M(), helicity_rec);
1268+
hMChists.fill(HIST("Recf1710_ptTemp"), multiplicity, mother1.Pt(), mother1.M(), helicity_rec2);
1269+
// hMChists.fill(HIST("RecRapidity"), mother.Y());
1270+
hMChists.fill(HIST("RecPhi"), mother.Phi());
1271+
hMChists.fill(HIST("RecEta"), mother.Eta());
12961272
}
12971273
gindex2.clear();
12981274
}

0 commit comments

Comments
 (0)