Skip to content

Commit f9e8a71

Browse files
author
Sawan Sawan
committed
modified MC
1 parent 38e18b0 commit f9e8a71

File tree

1 file changed

+87
-63
lines changed

1 file changed

+87
-63
lines changed

‎PWGLF/Tasks/Resonances/higherMassResonances.cxx‎

Lines changed: 87 additions & 63 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, mother1, motherRot, fourVecDauCM, fourVecDauCM1;
179+
ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, motherRot, fourVecDauCM;
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,17 +309,18 @@ 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, {{25, 0, 25}});
312+
hMChists.add("events_checkrec", "No. of events in the reconstructed MC", kTH1I, {{20, 0, 20}});
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_ptTemp", "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});
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_pt2", "Gen f_{0}(1710) p_{T}", kTH1F, {ptAxis});
320+
hMChists.add("Genf1710_mass2", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
321+
hMChists.add("Genf1710_pt2", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
321322
hMChists.add("GenEta", "Gen Eta", kTHnSparseF, {ptAxis, {100, -1.0f, 1.0f}});
322-
hMChists.add("GenPhi", "Gen Phi", kTH1F, {{70, 0.0, 7.0f}});
323+
hMChists.add("GenPhi", "Gen Phi", kTH1F, {{70, -3.5f, 3.5f}});
323324
hMChists.add("GenRapidity", "Gen Rapidity", kTHnSparseF, {ptAxis, {100, -1.0f, 1.0f}});
324325
hMChists.add("RecEta", "Rec Eta", kTH1F, {{100, -1.0f, 1.0f}});
325326
hMChists.add("RecPhi", "Rec Phi", kTH1F, {{70, 0.0f, 7.0f}});
@@ -589,7 +590,7 @@ struct HigherMassResonances {
589590
// For Monte Carlo
590591
using EventCandidatesMC = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Cs, aod::CentFT0Ms>;
591592
using TrackCandidatesMC = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullKa, aod::pidTOFFullKa, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::McTrackLabels>>;
592-
using V0TrackCandidatesMC = soa::Filtered<soa::Join<aod::V0Datas, aod::McV0Labels>>;
593+
using V0TrackCandidatesMC = soa::Join<aod::V0Datas, aod::McV0Labels>;
593594
// zBeam direction in lab frame
594595

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

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+
10541062
auto kDaughters = mcParticle.daughters_as<aod::McParticles>();
10551063
if (kDaughters.size() != 2) {
10561064
continue;
@@ -1075,26 +1083,25 @@ struct HigherMassResonances {
10751083
}
10761084
}
10771085
if (passKs.size() == 2) {
1078-
lResonance_gen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e());
1079-
lResonance_gen2 = daughter1 + daughter2;
1086+
lResonance_gen = ROOT::Math::PxPyPzMVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e());
1087+
lResonance_gen2 = daughter1 + daughter2; // invariant mass of Kshort pair
10801088

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

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

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());
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());
10921098
}
10931099
passKs.clear(); // clear the vector for the next iteration
10941100
}
10951101
}
10961102
PROCESS_SWITCH(HigherMassResonances, processGen, "Process Generated", false);
10971103

1104+
int counter2 = 0;
10981105
int eventCounter = 0;
10991106
std::vector<int> gindex1, gindex2;
11001107
void processRec(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const&, V0TrackCandidatesMC const& V0s, aod::McParticles const&, aod::McCollisions const& /*mcCollisions*/)
@@ -1103,6 +1110,7 @@ struct HigherMassResonances {
11031110
return;
11041111
}
11051112

1113+
ROOT::Math::PxPyPzMVector lDecayDaughter1, lDecayDaughter2, lResonance;
11061114
auto multiplicity = collision.centFT0C();
11071115
hMChists.fill(HIST("Rec_Multiplicity"), multiplicity);
11081116

@@ -1117,21 +1125,17 @@ struct HigherMassResonances {
11171125
}
11181126
hMChists.fill(HIST("events_checkrec"), 2.5);
11191127

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()) {
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))) {
11291133
return;
11301134
}
11311135
hMChists.fill(HIST("events_checkrec"), 4.5);
11321136
hMChists.fill(HIST("MC_mult_after_event_sel"), multiplicity);
11331137
eventCounter++;
1134-
// auto oldindex = -999;
1138+
auto oldindex = -999;
11351139

11361140
for (const auto& v01 : V0s) {
11371141

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

11651169
double nTPCSigmaPos1[1]{postrack1.tpcNSigmaPi()};
11661170
double nTPCSigmaNeg1[1]{negtrack1.tpcNSigmaPi()};
1171+
11671172
double nTPCSigmaPos2[1]{postrack2.tpcNSigmaPi()};
11681173
double nTPCSigmaNeg2[1]{negtrack2.tpcNSigmaPi()};
11691174

@@ -1202,73 +1207,92 @@ struct HigherMassResonances {
12021207
continue;
12031208
}
12041209
}
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;
12051215

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

12081218
hMChists.fill(HIST("events_checkrec"), 13.5);
12091219

1210-
if (mothertrack1.pdgCode() != config.pdgCodes[config.selectMCparticles]) {
1211-
continue;
1212-
}
1213-
hMChists.fill(HIST("events_checkrec"), 14.5);
1214-
12151220
if (mothertrack1.pdgCode() != mothertrack2.pdgCode()) {
12161221
continue;
12171222
}
1218-
hMChists.fill(HIST("events_checkrec"), 15.5);
1223+
hMChists.fill(HIST("events_checkrec"), 14.5);
12191224

1225+
// int motpdgs2 = std::abs(mothertrack2.pdgCode());
12201226
gindex2.push_back(mothertrack2.globalIndex());
12211227
if (gindex2.size() > 1) {
12221228
if (std::find(gindex2.begin(), gindex2.end(), mothertrack2.globalIndex()) != gindex2.end()) {
12231229
continue;
12241230
}
12251231
}
1226-
hMChists.fill(HIST("events_checkrec"), 16.5);
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);
12271239

12281240
if (mothertrack1.globalIndex() != mothertrack2.globalIndex()) {
12291241
continue;
12301242
}
1231-
hMChists.fill(HIST("events_checkrec"), 17.5);
1243+
hMChists.fill(HIST("events_checkrec"), 16.5);
12321244

12331245
if (!mothertrack1.producedByGenerator()) {
12341246
continue;
12351247
}
1236-
hMChists.fill(HIST("events_checkrec"), 18.5);
1248+
hMChists.fill(HIST("events_checkrec"), 17.5);
12371249

12381250
if (config.apply_rapidityMC && std::abs(mothertrack1.y()) >= 0.5) {
12391251
continue;
12401252
}
1241-
hMChists.fill(HIST("events_checkrec"), 19.5);
1242-
1243-
// if (config.avoidsplitrackMC && oldindex == mothertrack1.globalIndex()) {
1244-
// hMChists.fill(HIST("h1Recsplit"), mothertrack1.pt());
1245-
// continue;
1246-
// }
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()));
1253+
hMChists.fill(HIST("events_checkrec"), 18.5);
12641254

1265-
// std::cout << "mother pT is " << mother.Pt() << std::endl;
1255+
if (config.avoidsplitrackMC && oldindex == mothertrack1.globalIndex()) {
1256+
hMChists.fill(HIST("h1Recsplit"), mothertrack1.pt());
1257+
continue;
1258+
}
1259+
oldindex = mothertrack1.globalIndex();
12661260

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());
1261+
// counter_check++;
1262+
// if (counter_check > 1) {
1263+
// std::cout << "Total mothers is " << counter_check << std::endl;
1264+
// }
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());
12721296
}
12731297
gindex2.clear();
12741298
}

0 commit comments

Comments
 (0)