Skip to content

[PWGLF] "PWGLF: modified Reconstructed and Generated MC" #11931

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
dmallick2 merged 1 commit into from
Jul 5, 2025
Merged

[PWGLF] "PWGLF: modified Reconstructed and Generated MC" #11931

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
150 changes: 63 additions & 87 deletions PWGLF/Tasks/Resonances/higherMassResonances.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@
#include <TH1F.h>
#include <TH2F.h>
#include <THn.h>
#include <TLorentzVector.h>

Check failure on line 45 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
#include <TMath.h>
#include <TObjArray.h>
#include <TPDGCode.h>
Expand Down Expand Up @@ -176,12 +176,12 @@
// variables declaration
float multiplicity = 0.0f;
float theta2;
ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, motherRot, fourVecDauCM;
ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, mother1, motherRot, fourVecDauCM, fourVecDauCM1;
ROOT::Math::XYZVector randomVec, beamVec, normalVec;
ROOT::Math::XYZVectorF v1_CM, zaxis_HE, yaxis_HE, xaxis_HE;
// ROOT::Math::XYZVector threeVecDauCM, helicityVec, randomVec, beamVec, normalVec;
ROOT::Math::XYZVector zBeam; // ẑ: beam direction in lab frame
double BeamMomentum = TMath::Sqrt(13600 * 13600 / 4 - 0.938 * 0.938); // GeV

Check failure on line 184 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-mass] 

Avoid hard-coded particle masses. Use o2::constants::physics::Mass... instead.

Check failure on line 184 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/entity] 

Replace ROOT entities with equivalents from standard C++ or from O2.
ROOT::Math::PxPyPzEVector Beam1{0., 0., -BeamMomentum, 13600. / 2.};
ROOT::Math::PxPyPzEVector Beam2{0., 0., BeamMomentum, 13600. / 2.};
ROOT::Math::XYZVectorF Beam1_CM, Beam2_CM;
Expand Down Expand Up @@ -309,18 +309,17 @@
// For MC
if (config.isMC) {
hMChists.add("events_check", "No. of events in the generated MC", kTH1I, {{20, 0, 20}});
hMChists.add("events_checkrec", "No. of events in the reconstructed MC", kTH1I, {{20, 0, 20}});
hMChists.add("events_checkrec", "No. of events in the reconstructed MC", kTH1I, {{25, 0, 25}});
hMChists.add("Genf1710", "Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL});
hMChists.add("Recf1710_pt1", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
hMChists.add("Recf1710_pt2", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
hMChists.add("Recf1710_ptTemp", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
// hMChists.add("Recf1710_p", "Rec f_{0}(1710) p", kTH1F, {ptAxis});
hMChists.add("h1Recsplit", "Rec p_{T}2", kTH1F, {ptAxis});
// hMChists.add("Recf1710_mass", "Rec f_{0}(1710) mass", kTH1F, {glueballMassAxis});
hMChists.add("Genf1710_mass", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
hMChists.add("Genf1710_mass2", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
hMChists.add("Genf1710_pt2", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
hMChists.add("Genf1710_pt2", "Gen f_{0}(1710) p_{T}", kTH1F, {ptAxis});
hMChists.add("GenEta", "Gen Eta", kTHnSparseF, {ptAxis, {100, -1.0f, 1.0f}});
hMChists.add("GenPhi", "Gen Phi", kTH1F, {{70, -3.5f, 3.5f}});
hMChists.add("GenPhi", "Gen Phi", kTH1F, {{70, 0.0, 7.0f}});
hMChists.add("GenRapidity", "Gen Rapidity", kTHnSparseF, {ptAxis, {100, -1.0f, 1.0f}});
hMChists.add("RecEta", "Rec Eta", kTH1F, {{100, -1.0f, 1.0f}});
hMChists.add("RecPhi", "Rec Phi", kTH1F, {{70, 0.0f, 7.0f}});
Expand Down Expand Up @@ -378,7 +377,7 @@
const float cpav0 = candidate.v0cosPA();

float ctauK0s = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassK0Short;
float lowmasscutks0 = 0.497 - config.cWidthKs0 * config.cSigmaMassKs0;

Check failure on line 380 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-mass] 

Avoid hard-coded particle masses. Use o2::constants::physics::Mass... instead.
float highmasscutks0 = 0.497 + config.cWidthKs0 * config.cSigmaMassKs0;
// float decayLength = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * RecoDecay::sqrtSumOfSquares(candidate.px(), candidate.py(), candidate.pz());

Expand Down Expand Up @@ -590,7 +589,7 @@
// For Monte Carlo
using EventCandidatesMC = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Cs, aod::CentFT0Ms>;
using TrackCandidatesMC = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullKa, aod::pidTOFFullKa, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::McTrackLabels>>;
using V0TrackCandidatesMC = soa::Join<aod::V0Datas, aod::McV0Labels>;
using V0TrackCandidatesMC = soa::Filtered<soa::Join<aod::V0Datas, aod::McV0Labels>>;
// zBeam direction in lab frame

template <typename T>
Expand Down Expand Up @@ -636,9 +635,9 @@

// CosThetaHE = zaxis_HE.Dot(v_CM);

auto angle_phi = TMath::ATan2(yaxis_HE.Dot(v1_CM), xaxis_HE.Dot(v1_CM));

Check failure on line 638 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/entity] 

Replace ROOT entities with equivalents from standard C++ or from O2.
if (angle_phi < 0) {
angle_phi += 2 * TMath::Pi(); // ensure phi is in [0, 2pi]

Check failure on line 640 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pi-multiple-fraction] 

Use multiples/fractions of PI defined in o2::constants::math.

Check failure on line 640 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[two-pi-add-subtract] 

Use RecoDecay::constrainAngle to restrict angle to a given range.

Check failure on line 640 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[external-pi] 

Use the PI constant (and its multiples and fractions) defined in o2::constants::math.
}

if (std::abs(mother.Rapidity()) < 0.5) {
Expand Down Expand Up @@ -1052,13 +1051,6 @@
// std::cout << "px " << mcParticle.px() << " py " << mcParticle.py() << " pz " << mcParticle.pz() << " y " << mcParticle.y() << std::endl;
// counter++;

hMChists.fill(HIST("GenRapidity"), mcParticle.pt(), mcParticle.y());
hMChists.fill(HIST("GenPhi"), mcParticle.phi());
hMChists.fill(HIST("GenEta"), mcParticle.pt(), mcParticle.eta());
// hMChists.fill(HIST("GenPx"), mcParticle.px());
// hMChists.fill(HIST("GenPy"), mcParticle.py());
// hMChists.fill(HIST("GenPz"), mcParticle.pz());

auto kDaughters = mcParticle.daughters_as<aod::McParticles>();
if (kDaughters.size() != 2) {
continue;
Expand All @@ -1072,7 +1064,7 @@
continue;
}
hMChists.fill(HIST("events_check"), 8.5);
if (std::abs(kCurrentDaughter.pdgCode()) == 310) {

Check failure on line 1067 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
passKs.push_back(true);
hMChists.fill(HIST("events_check"), 9.5);
if (passKs.size() == 1) {
Expand All @@ -1083,25 +1075,26 @@
}
}
if (passKs.size() == 2) {
lResonance_gen = ROOT::Math::PxPyPzMVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e());
lResonance_gen2 = daughter1 + daughter2; // invariant mass of Kshort pair
lResonance_gen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e());
lResonance_gen2 = daughter1 + daughter2;

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

auto helicity_gen = lResonance_gen.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonance_gen.Vect().Mag2()));
auto helicity_gen = lResonance_gen2.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonance_gen2.Vect().Mag2()));

hMChists.fill(HIST("Genf1710"), multiplicityGen, mcParticle.pt(), helicity_gen);
hMChists.fill(HIST("Genf1710_mass"), lResonance_gen.M());
hMChists.fill(HIST("Genf1710_mass2"), lResonance_gen2.M());
hMChists.fill(HIST("Genf1710_pt2"), lResonance_gen2.Pt());
hMChists.fill(HIST("Genf1710"), multiplicityGen, lResonance_gen2.pt(), helicity_gen);
hMChists.fill(HIST("Genf1710_mass"), lResonance_gen2.M());
hMChists.fill(HIST("Genf1710_pt2"), mcParticle.pt());
hMChists.fill(HIST("GenRapidity"), lResonance_gen2.Pt(), lResonance_gen2.Y());
hMChists.fill(HIST("GenPhi"), lResonance_gen2.Phi());
hMChists.fill(HIST("GenEta"), lResonance_gen2.Pt(), lResonance_gen2.Eta());
}
passKs.clear(); // clear the vector for the next iteration
}
}
PROCESS_SWITCH(HigherMassResonances, processGen, "Process Generated", false);

int counter2 = 0;
int eventCounter = 0;
std::vector<int> gindex1, gindex2;
void processRec(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const&, V0TrackCandidatesMC const& V0s, aod::McParticles const&, aod::McCollisions const& /*mcCollisions*/)
Expand All @@ -1110,7 +1103,6 @@
return;
}

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

Expand All @@ -1125,17 +1117,21 @@
}
hMChists.fill(HIST("events_checkrec"), 2.5);

if (config.timFrameEvsel && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
return;
}
hMChists.fill(HIST("events_checkrec"), 3.5);
if (config.cTVXEvsel && (!collision.selection_bit(aod::evsel::kIsTriggerTVX))) {
// if (config.timFrameEvsel && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
// return;
// }
// hMChists.fill(HIST("events_checkrec"), 3.5);
// if (config.cTVXEvsel && (!collision.selection_bit(aod::evsel::kIsTriggerTVX))) {
// return;
// }

if (!collision.sel8()) {
return;
}
hMChists.fill(HIST("events_checkrec"), 4.5);
hMChists.fill(HIST("MC_mult_after_event_sel"), multiplicity);
eventCounter++;
auto oldindex = -999;
// auto oldindex = -999;

for (const auto& v01 : V0s) {

Expand Down Expand Up @@ -1168,7 +1164,6 @@

double nTPCSigmaPos1[1]{postrack1.tpcNSigmaPi()};
double nTPCSigmaNeg1[1]{negtrack1.tpcNSigmaPi()};

double nTPCSigmaPos2[1]{postrack2.tpcNSigmaPi()};
double nTPCSigmaNeg2[1]{negtrack2.tpcNSigmaPi()};

Expand All @@ -1193,7 +1188,7 @@
int trackv0PDG1 = std::abs(mctrackv01.pdgCode());
int trackv0PDG2 = std::abs(mctrackv02.pdgCode());

if (std::abs(trackv0PDG1) != 310 || std::abs(trackv0PDG2) != 310) {

Check failure on line 1191 in PWGLF/Tasks/Resonances/higherMassResonances.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[pdg/explicit-code] 

Avoid hard-coded PDG codes. Use named values from PDG_t or o2::constants::physics::Pdg instead.
continue;
}
hMChists.fill(HIST("events_checkrec"), 12.5);
Expand All @@ -1207,92 +1202,73 @@
continue;
}
}
// if (counter2 < 1e4)
// std::cout << "Mother1 pdg code: " << motpdgs << " p_{T} " << mothertrack1.pt() << "Global index " << mothertrack1.globalIndex() << " event " << eventCounter << std::endl;
// counter2++;

// int counter_check = 0;

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

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

if (mothertrack1.pdgCode() != mothertrack2.pdgCode()) {
if (mothertrack1.pdgCode() != config.pdgCodes[config.selectMCparticles]) {
continue;
}
hMChists.fill(HIST("events_checkrec"), 14.5);

// int motpdgs2 = std::abs(mothertrack2.pdgCode());
if (mothertrack1.pdgCode() != mothertrack2.pdgCode()) {
continue;
}
hMChists.fill(HIST("events_checkrec"), 15.5);

gindex2.push_back(mothertrack2.globalIndex());
if (gindex2.size() > 1) {
if (std::find(gindex2.begin(), gindex2.end(), mothertrack2.globalIndex()) != gindex2.end()) {
continue;
}
}
// if (counter2 < 1e4)
// std::cout << "Mother2 pdg code: " << motpdgs2 << " p_{T} " << mothertrack2.pt() << "Global index " << mothertrack1.globalIndex() << " event " << eventCounter << std::endl;

if (mothertrack1.pdgCode() != config.pdgCodes[config.selectMCparticles]) {
continue;
}
hMChists.fill(HIST("events_checkrec"), 15.5);

if (mothertrack1.globalIndex() != mothertrack2.globalIndex()) {
continue;
}
hMChists.fill(HIST("events_checkrec"), 16.5);

if (!mothertrack1.producedByGenerator()) {
if (mothertrack1.globalIndex() != mothertrack2.globalIndex()) {
continue;
}
hMChists.fill(HIST("events_checkrec"), 17.5);

if (config.apply_rapidityMC && std::abs(mothertrack1.y()) >= 0.5) {
if (!mothertrack1.producedByGenerator()) {
continue;
}
hMChists.fill(HIST("events_checkrec"), 18.5);

if (config.avoidsplitrackMC && oldindex == mothertrack1.globalIndex()) {
hMChists.fill(HIST("h1Recsplit"), mothertrack1.pt());
if (config.apply_rapidityMC && std::abs(mothertrack1.y()) >= 0.5) {
continue;
}
oldindex = mothertrack1.globalIndex();
hMChists.fill(HIST("events_checkrec"), 19.5);

// counter_check++;
// if (counter_check > 1) {
// std::cout << "Total mothers is " << counter_check << std::endl;
// if (config.avoidsplitrackMC && oldindex == mothertrack1.globalIndex()) {
// hMChists.fill(HIST("h1Recsplit"), mothertrack1.pt());
// continue;
// }
// std::cout << "After selection " << " p_{T} " << mothertrack2.pt() << " event " << eventCounter << std::endl;

pvec0 = std::array{v01.px(), v01.py(), v01.pz()};
pvec1 = std::array{v02.px(), v02.py(), v02.pz()};
auto arrMomrec = std::array{pvec0, pvec1};
// auto motherP = mothertrack1.p();
// auto motherE = mothertrack1.e();
// auto genMass = std::sqrt(motherE * motherE - motherP * motherP);
auto recMass = RecoDecay::m(arrMomrec, std::array{o2::constants::physics::MassK0Short, o2::constants::physics::MassK0Short});
// auto recpt = TMath::Sqrt((track1.px() + track2.px()) * (track1.px() + track2.px()) + (track1.py() + track2.py()) * (track1.py() + track2.py()));
//// Resonance reconstruction
lDecayDaughter1 = ROOT::Math::PxPyPzMVector(v01.px(), v01.py(), v01.pz(), o2::constants::physics::MassK0Short);
lDecayDaughter2 = ROOT::Math::PxPyPzMVector(v02.px(), v02.py(), v02.pz(), o2::constants::physics::MassK0Short);
lResonance = lDecayDaughter1 + lDecayDaughter2;
if (config.apply_rapidityMC && std::abs(lResonance.Y()) >= 0.5) {
continue;
}
// daughter1, mother, fourVecDauCM
ROOT::Math::Boost boost{lResonance.BoostToCM()};
fourVecDauCM = boost(lDecayDaughter1); // boost the frame of daughter to the center of mass frame
auto helicity_rec = lResonance.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonance.Vect().Mag2()));

// hMChists.fill(HIST("Recf1710_p"), motherP);
// hMChists.fill(HIST("Recf1710_mass"), recMass);
hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mothertrack1.pt(), recMass, helicity_rec);
// hMChists.fill(HIST("Genf1710_mass"), genMass);
hMChists.fill(HIST("Recf1710_pt2"), multiplicity, lResonance.Pt(), recMass, helicity_rec);

hMChists.fill(HIST("RecRapidity"), mothertrack1.y());
hMChists.fill(HIST("RecPhi"), mothertrack1.phi());
hMChists.fill(HIST("RecEta"), mothertrack1.eta());
// hMChists.fill(HIST("events_checkrec"), 20.5);
// oldindex = mothertrack1.globalIndex(); // split tracks is already handled using gindex1 and gindex2

daughter1 = ROOT::Math::PxPyPzMVector(v01.px(), v01.py(), v01.pz(), o2::constants::physics::MassK0Short);
daughter2 = ROOT::Math::PxPyPzMVector(v02.px(), v02.py(), v02.pz(), o2::constants::physics::MassK0Short);
mother = daughter1 + daughter2;
mother1 = ROOT::Math::PxPyPzEVector(mothertrack1.px(), mothertrack1.py(), mothertrack1.pz(), mothertrack1.e());

ROOT::Math::Boost boost{mother.BoostToCM()};
ROOT::Math::Boost boost1{mother1.BoostToCM()};

fourVecDauCM = boost(daughter1);
fourVecDauCM1 = boost1(daughter1);

auto helicity_rec = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2()));

auto helicity_rec2 = mother1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(mother1.Vect().Mag2()));

// std::cout << "mother pT is " << mother.Pt() << std::endl;

hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mother.Pt(), mother.M(), helicity_rec);
hMChists.fill(HIST("Recf1710_ptTemp"), multiplicity, mother1.Pt(), mother1.M(), helicity_rec2);
// hMChists.fill(HIST("RecRapidity"), mother.Y());
hMChists.fill(HIST("RecPhi"), mother.Phi());
hMChists.fill(HIST("RecEta"), mother.Eta());
}
gindex2.clear();
}
Expand Down
Loading