From 9fe4f4ae2e25ad95b4537df33f0205596a5d3dec Mon Sep 17 00:00:00 2001 From: Sawan Sawan Date: Mon, 28 Jul 2025 20:15:00 +0530 Subject: [PATCH 1/5] added signal/event loss corrections --- PWGLF/Tasks/Resonances/kstarqa.cxx | 106 ++++++++++++++++++----------- 1 file changed, 68 insertions(+), 38 deletions(-) diff --git a/PWGLF/Tasks/Resonances/kstarqa.cxx b/PWGLF/Tasks/Resonances/kstarqa.cxx index 078fc6916bf..1bddc1a379b 100644 --- a/PWGLF/Tasks/Resonances/kstarqa.cxx +++ b/PWGLF/Tasks/Resonances/kstarqa.cxx @@ -147,6 +147,7 @@ struct Kstarqa { ConfigurableAxis configThnAxisPOL{"configThnAxisPOL", {20, -1.0, 1.0}, "Costheta axis"}; ConfigurableAxis invMassKstarAxis{"invMassKstarAxis", {300, 0.7f, 1.3f}, "Kstar invariant mass axis"}; ConfigurableAxis ptAxisKstar{"ptAxisKstar", {200, 0.0f, 20.0f}, "Kstar pT axis"}; + ConfigurableAxis binsImpactPar{"binsImpactPar", {100, 0, 25}, "Binning of the impact parameter axis"}; // Event plane configurables Configurable boostDaugter1{"boostDaugter1", false, "Boost daughter Kaon in the COM frame"}; @@ -167,6 +168,7 @@ struct Kstarqa { AxisSpec invmassAxis = {invMassKstarAxis, "Invariant mass (GeV/#it{c}^{2})"}; AxisSpec thnAxisPOL{configThnAxisPOL, "cos(#theta)"}; AxisSpec multiplicityAxis = {binsMultPlot, "Multiplicity Axis"}; + AxisSpec impactParAxis = {binsImpactPar, "Impact Parameter (cm)"}; // Histograms // Event selection @@ -256,12 +258,22 @@ struct Kstarqa { hInvMass.add("h1genmass", "Invariant mass of generated kstar meson", kTH1F, {invmassAxis}); hInvMass.add("h1GenMult", "Multiplicity generated", kTH1F, {multiplicityAxis}); hInvMass.add("h1RecMult", "Multiplicity reconstructed", kTH1F, {multiplicityAxis}); + hInvMass.add("h1KSRecsplit", "KS meson Rec split", kTH1F, {{100, 0.0f, 10.0f}}); + hInvMass.add("MCcorrections/hSignalLossDenominator", "Kstar generated before event selection", kTH2F, {{ptAxis}, {impactParAxis}}); + hInvMass.add("MCcorrections/hSignalLossNumerator", "Kstar generated after event selection", kTH2F, {{ptAxis}, {impactParAxis}}); + // hInvMass.add("hAllGenCollisionsImpact", "All generated collisions vs impact parameter", kTH1F, {multiplicityAxis}); + hInvMass.add("hAllGenCollisions", "All generated events", kTH1F, {multiplicityAxis}); + hInvMass.add("hAllGenCollisions1Rec", "All gen events with at least one rec event", kTH1F, {multiplicityAxis}); + hInvMass.add("hAllKstarGenCollisisons", "All generated Kstar in events with rapidity in 0.5", kTH2F, {{multiplicityAxis}, {ptAxis}}); + hInvMass.add("hAllKstarGenCollisisons1Rec", "All generated Kstar in events with at least one rec event in rapidity in 0.5", kTH2F, {{multiplicityAxis}, {ptAxis}}); + hInvMass.add("hAllRecCollisions", "All reconstructed events", kTH2F, {{multiplicityAxis}}); + hInvMass.add("MCcorrections/hImpactParameterRec", "Impact parameter in reconstructed MC", kTH1F, {{impactParAxis}}); + hInvMass.add("MCcorrections/hImpactParameterGen", "Impact parameter in generated MC", kTH1F, {{impactParAxis}}); + hInvMass.add("MCcorrections/hImpactParametervsMultiplicity", "Impact parameter vs multiplicity in reconstructed MC", kTH1F, {{impactParAxis}, {multiplicityAxis}}); rEventSelection.add("events_check_data", "No. of events in the data", kTH1I, {{20, 0, 20}}); rEventSelection.add("events_check", "No. of events in the generated MC", kTH1I, {{20, 0, 20}}); rEventSelection.add("events_checkrec", "No. of events in the reconstructed MC", kTH1I, {{20, 0, 20}}); - hInvMass.add("h1KSRecsplit", "KS meson Rec split", kTH1F, {{100, 0.0f, 10.0f}}); - hInvMass.add("kstargenBeforeEvtSel", "Kstar generated before event selection", kTH1F, {ptAxis}); - hInvMass.add("kstargenAfterEvtSel", "Kstar generated after event selection", kTH1F, {ptAxis}); + rEventSelection.add("hOccupancy", "Occupancy distribution", kTH1F, {{1000, 0, 15000}}); // Multplicity distribution if (cQAevents) { @@ -750,9 +762,11 @@ struct Kstarqa { if (!selectionEvent(collision, true)) { return; } - rEventSelection.fill(HIST("events_check_data"), 3.5); + int occupancy = collision.trackOccupancyInTimeRange(); + rEventSelection.fill(HIST("hOccupancy"), occupancy); + multiplicity = -1; if (cSelectMultEstimator == 0) { @@ -993,6 +1007,12 @@ struct Kstarqa { std::vector selectedEvents(collisions.size()); int nevts = 0; multiplicity = -1.0; + // float impactParameter = mcCollision.impactParameter(); + + // if (mcCollision.isInelGt0()) { + // return; + // } + for (const auto& collision : collisions) { // if (!collision.sel8() || std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex) { if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex) { @@ -1019,15 +1039,24 @@ struct Kstarqa { multiplicity = collision.centFT0M(); hInvMass.fill(HIST("h1GenMult"), multiplicity); selectedEvents[nevts++] = collision.mcCollision_as().globalIndex(); + int occupancy = collision.trackOccupancyInTimeRange(); + rEventSelection.fill(HIST("hOccupancy"), occupancy); } selectedEvents.resize(nevts); rEventSelection.fill(HIST("events_check"), 3.5); - const auto evtReconstructedAndSelected = std::find(selectedEvents.begin(), selectedEvents.end(), mcCollision.globalIndex()) != selectedEvents.end(); + for (const auto& mcParticle : mcParticles) { + if (std::abs(mcParticle.y()) < 0.5 && std::abs(mcParticle.pdgCode()) == o2::constants::physics::kK0Star892) { + hInvMass.fill(HIST("hAllKstarGenCollisisons"), multiplicity, mcParticle.pt()); + } + } + const auto evtReconstructedAndSelected = std::find(selectedEvents.begin(), selectedEvents.end(), mcCollision.globalIndex()) != selectedEvents.end(); + hInvMass.fill(HIST("hAllGenCollisions"), multiplicity); if (!cAllGenCollisions && !evtReconstructedAndSelected) { // Check that the event is reconstructed and that the reconstructed events pass the selection return; } + hInvMass.fill(HIST("hAllGenCollisions1Rec"), multiplicity); rEventSelection.fill(HIST("events_check"), 4.5); for (const auto& mcParticle : mcParticles) { @@ -1040,6 +1069,7 @@ struct Kstarqa { continue; } rEventSelection.fill(HIST("events_check"), 6.5); + hInvMass.fill(HIST("hAllKstarGenCollisisons1Rec"), multiplicity, mcParticle.pt()); auto kDaughters = mcParticle.daughters_as(); if (kDaughters.size() != 2) { @@ -1073,43 +1103,41 @@ struct Kstarqa { } } PROCESS_SWITCH(Kstarqa, processGen, "Process Generated", false); - /* - void processEvtLossSigLossMC(aod::McCollisions::iterator const&, aod::McParticles const& mcParticles, const soa::SmallGroups& recCollisions) - { - - bool isSel = false; - // auto multiplicity1 = -999.; - for (const auto& RecCollision : recCollisions) { - if (!selectionEvent(RecCollision, false)) - continue; - // if (cSelectMultEstimator == 0) { - // multiplicity1 = RecCollision.centFT0M(); - // } else if (cSelectMultEstimator == 1) { - // multiplicity1 = RecCollision.centFT0A(); - // } else if (cSelectMultEstimator == 2) { - // multiplicity1 = RecCollision.centFT0C(); - // } else { - // multiplicity1 = RecCollision.centFT0M(); - // } - - isSel = true; - } + void processEvtLossSigLossMC(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& recCollisions) + { + auto impactPar = mcCollision.impactParameter(); + hInvMass.fill(HIST("MCcorrections/hImpactParameterGen"), impactPar); - // Generated MC - for (const auto& mcPart : mcParticles) { - if (std::abs(mcPart.y()) >= 0.5 || std::abs(mcPart.pdgCode()) != o2::constants::physics::kK0Star892) - continue; + bool isSelectedEvent = false; + auto multiplicity1 = -999.; + for (const auto& RecCollision : recCollisions) { + if (!selectionEvent(RecCollision, false)) + continue; + multiplicity1 = RecCollision.centFT0M(); + isSelectedEvent = true; + } - // signal loss estimation - hInvMass.fill(HIST("kstargenBeforeEvtSel"), mcPart.pt()); - if (isSel) { - hInvMass.fill(HIST("kstargenAfterEvtSel"), mcPart.pt()); - } - } // end loop on gen particles + // Event loss + if (isSelectedEvent) { + hInvMass.fill(HIST("MCcorrections/hImpactParameterRec"), impactPar); + hInvMass.fill(HIST("MCcorrections/hImpactParametervsMultiplicity"), impactPar, multiplicity1); } - PROCESS_SWITCH(Kstarqa, processEvtLossSigLossMC, "Process Signal Loss, Event Loss", false); - */ + + // Generated MC + for (const auto& mcPart : mcParticles) { + if (std::abs(mcPart.y()) >= 0.5 || std::abs(mcPart.pdgCode()) != o2::constants::physics::kK0Star892) + continue; + + // signal loss estimation + hInvMass.fill(HIST("MCcorrections/hSignalLossDenominator"), mcPart.pt(), impactPar); + if (isSelectedEvent) { + hInvMass.fill(HIST("MCcorrections/hSignalLossNumerator"), mcPart.pt(), impactPar); + } + } // end loop on gen particles + } + PROCESS_SWITCH(Kstarqa, processEvtLossSigLossMC, "Process Signal Loss, Event Loss", false); + void processRec(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&, aod::McCollisions const& /*mcCollisions*/) { @@ -1123,6 +1151,8 @@ struct Kstarqa { if (selectionConfig.isINELgt0 && !collision.isInelGt0()) { return; } + multiplicity = collision.centFT0M(); + hInvMass.fill(HIST("hAllRecCollisions"), multiplicity); // if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex || !collision.sel8()) { if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex) { From 7c27e70701c6fd23f2ad094e55ecc7d4a53dc0ec Mon Sep 17 00:00:00 2001 From: Sawan Sawan Date: Mon, 28 Jul 2025 20:19:14 +0530 Subject: [PATCH 2/5] some minor corrections --- PWGLF/Tasks/Resonances/kstarqa.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Resonances/kstarqa.cxx b/PWGLF/Tasks/Resonances/kstarqa.cxx index 1bddc1a379b..7c01c7b63ec 100644 --- a/PWGLF/Tasks/Resonances/kstarqa.cxx +++ b/PWGLF/Tasks/Resonances/kstarqa.cxx @@ -1327,7 +1327,7 @@ struct Kstarqa { hInvMass.fill(HIST("h2KstarRecpt2"), mothertrack1.pt(), multiplicity, std::sqrt(mothertrack1.e() * mothertrack1.e() - mothertrack1.p() * mothertrack1.p())); - if (applyRecMotherRapidity && mother.Rapidity() >= 0) { + if (applyRecMotherRapidity && mother.Rapidity() >= 0.5) { continue; } From 9ec6f95fc5b773bc43669d8ec86da8ad9bfc7ed0 Mon Sep 17 00:00:00 2001 From: Sawan Sawan Date: Mon, 28 Jul 2025 22:34:03 +0530 Subject: [PATCH 3/5] rapidity correction --- .../Tasks/Resonances/higherMassResonances.cxx | 111 +++++++++--------- 1 file changed, 55 insertions(+), 56 deletions(-) diff --git a/PWGLF/Tasks/Resonances/higherMassResonances.cxx b/PWGLF/Tasks/Resonances/higherMassResonances.cxx index 61ef930a7c2..535c61e91d1 100644 --- a/PWGLF/Tasks/Resonances/higherMassResonances.cxx +++ b/PWGLF/Tasks/Resonances/higherMassResonances.cxx @@ -43,7 +43,6 @@ #include #include #include -#include #include #include #include @@ -143,7 +142,7 @@ struct HigherMassResonances { Configurable cTVXEvsel{"cTVXEvsel", true, "Triggger selection"}; Configurable avoidsplitrackMC{"avoidsplitrackMC", false, "avoid split track in MC"}; Configurable selectMCparticles{"selectMCparticles", 1, "0: f0(1710), 1: f2(1525), 2: a2(1320), 3: f0(1370), 4: f0(1500)"}; - Configurable apply_rapidityMC{"apply_rapidityMC", true, "Apply rapidity cut on generated and reconstructed particles"}; + Configurable applyRapidityMC{"applyRapidityMC", true, "Apply rapidity cut on generated and reconstructed particles"}; std::vector pdgCodes = {10331, 335, 115, 10221, 9030221}; // output THnSparses @@ -182,13 +181,13 @@ struct HigherMassResonances { ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, motherRot, fourVecDauCM, fourVecDauCM1; ROOT::Math::PxPyPzEVector mother1; ROOT::Math::XYZVector randomVec, beamVec, normalVec; - ROOT::Math::XYZVectorF v1_CM, zaxis_HE, yaxis_HE, xaxis_HE; + ROOT::Math::XYZVectorF v1CM, zaxisHE, yaxisHE, xaxisHE; // ROOT::Math::XYZVector threeVecDauCM, helicityVec, randomVec, beamVec, normalVec; ROOT::Math::XYZVector zBeam; // ẑ: beam direction in lab frame - double BeamMomentum = std::sqrt(13600 * 13600 / 4 - 0.938 * 0.938); // GeV - 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; + double beamMomentum = std::sqrt(13600 * 13600 / 4 - 0.938 * 0.938); // GeV + ROOT::Math::PxPyPzEVector beam1{0., 0., -beamMomentum, 13600. / 2.}; + ROOT::Math::PxPyPzEVector beam2{0., 0., beamMomentum, 13600. / 2.}; + ROOT::Math::XYZVectorF beam1CM, beam2CM; // const double massK0s = o2::constants::physics::MassK0Short; bool isMix = false; @@ -615,8 +614,8 @@ struct HigherMassResonances { fourVecDauCM = boost(daughter1); // boost the frame of daughter to the center of mass frame // threeVecDauCM = fourVecDauCM.Vect(); // get the 3 vector of daughter in the frame of mother - Beam1_CM = ROOT::Math::XYZVectorF((boost(Beam1).Vect()).Unit()); - Beam2_CM = ROOT::Math::XYZVectorF((boost(Beam2).Vect()).Unit()); + beam1CM = ROOT::Math::XYZVectorF((boost(beam1).Vect()).Unit()); + beam2CM = ROOT::Math::XYZVectorF((boost(beam2).Vect()).Unit()); // define y = zBeam x z: Normal to the production plane // ẑ: mother direction in lab, boosted into mother's rest frame @@ -634,22 +633,22 @@ struct HigherMassResonances { // auto p_proj_y = threeVecDauCM.Dot(y_axis); // // Calculate φ in [-π, π] - // auto angle_phi = std::atan2(p_proj_y, p_proj_x); // φ in radians + // auto anglePhi = std::atan2(p_proj_y, p_proj_x); // φ in radians - v1_CM = ROOT::Math::XYZVectorF(boost(daughter1).Vect()).Unit(); + v1CM = ROOT::Math::XYZVectorF(boost(daughter1).Vect()).Unit(); // ROOT::Math::XYZVectorF v2_CM{(boost(daughter1).Vect()).Unit()}; // using positive sign convention for the first track - // 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 + // 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 // Helicity frame - zaxis_HE = ROOT::Math::XYZVectorF(mother.Vect()).Unit(); - yaxis_HE = ROOT::Math::XYZVectorF(Beam1_CM.Cross(Beam2_CM)).Unit(); - xaxis_HE = ROOT::Math::XYZVectorF(yaxis_HE.Cross(zaxis_HE)).Unit(); + zaxisHE = ROOT::Math::XYZVectorF(mother.Vect()).Unit(); + yaxisHE = ROOT::Math::XYZVectorF(beam1CM.Cross(beam2CM)).Unit(); + xaxisHE = ROOT::Math::XYZVectorF(yaxisHE.Cross(zaxisHE)).Unit(); - // CosThetaHE = zaxis_HE.Dot(v_CM); + // CosThetaHE = zaxisHE.Dot(v_CM); - auto angle_phi = TMath::ATan2(yaxis_HE.Dot(v1_CM), xaxis_HE.Dot(v1_CM)); - if (angle_phi < 0) { - angle_phi += 2 * TMath::Pi(); // ensure phi is in [0, 2pi] + auto anglePhi = std::atan2(yaxisHE.Dot(v1CM), xaxisHE.Dot(v1CM)); + if (anglePhi < 0) { + anglePhi += 2 * o2::constants::math::PI; // ensure phi is in [0, 2pi] } // if (std::abs(mother.Rapidity()) < 0.5) { @@ -659,7 +658,7 @@ struct HigherMassResonances { auto cosThetaStarHelicity = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2())); if (!isMix) { if (std::abs(mother.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, angle_phi); + hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi); } for (int i = 0; i < config.cRotations; i++) { @@ -674,11 +673,11 @@ struct HigherMassResonances { auto cosThetaStarHelicityRot = motherRot.Vect().Dot(daughterRotCM.Vect()) / (std::sqrt(daughterRotCM.Vect().Mag2()) * std::sqrt(motherRot.Vect().Mag2())); if (motherRot.Rapidity() < 0.5) - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, angle_phi); + hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, anglePhi); } } else { if (std::abs(mother.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, angle_phi); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi); } } } else if (config.activateTHnSparseCosThStarProduction) { @@ -686,18 +685,18 @@ struct HigherMassResonances { auto cosThetaStarProduction = normalVec.Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(normalVec.Mag2())); if (!isMix) { if (std::abs(mother.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, angle_phi); + hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi); } for (int i = 0; i < config.cRotations; i++) { theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); 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()); if (std::abs(motherRot.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarProduction, angle_phi); + hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarProduction, anglePhi); } } } else { if (std::abs(mother.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, angle_phi); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi); } } } else if (config.activateTHnSparseCosThStarBeam) { @@ -705,18 +704,18 @@ struct HigherMassResonances { auto cosThetaStarBeam = beamVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2()); if (!isMix) { if (std::abs(mother.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, angle_phi); + hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi); } for (int i = 0; i < config.cRotations; i++) { theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); 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()); if (std::abs(motherRot.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam, angle_phi); + hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam, anglePhi); } } } else { if (std::abs(mother.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, angle_phi); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi); } } } else if (config.activateTHnSparseCosThStarRandom) { @@ -727,18 +726,18 @@ struct HigherMassResonances { auto cosThetaStarRandom = randomVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2()); if (!isMix) { if (std::abs(mother.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, angle_phi); + hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi); } for (int i = 0; i < config.cRotations; i++) { theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); 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()); if (std::abs(motherRot.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, angle_phi); + hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, anglePhi); } } } else { if (std::abs(mother.Rapidity()) < 0.5) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, angle_phi); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi); } } } @@ -864,7 +863,7 @@ struct HigherMassResonances { using EventCandidatesDerivedData = soa::Join; using V0CandidatesDerivedData = soa::Join; - using dauTracks = soa::Join; + using DauTracks = soa::Join; void processSEderived(EventCandidatesDerivedData::iterator const& collision, TrackCandidates const& /*tracks*/, aod::V0Datas const& V0s) { @@ -1290,8 +1289,8 @@ struct HigherMassResonances { int counter = 0; float multiplicityGen = 0.0; std::vector passKs; - ROOT::Math::PxPyPzMVector lResonance_gen1; - ROOT::Math::PxPyPzEVector lResonance_gen; + ROOT::Math::PxPyPzMVector lResonanceGen1; + ROOT::Math::PxPyPzEVector lResonanceGen; void processGen(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& collisions) { @@ -1349,7 +1348,7 @@ struct HigherMassResonances { } hMChists.fill(HIST("events_check"), 5.5); - if (config.apply_rapidityMC && std::abs(mcParticle.y()) >= 0.5) { + if (config.applyRapidityMC && std::abs(mcParticle.y()) >= 0.5) { continue; } hMChists.fill(HIST("events_check"), 6.5); @@ -1382,33 +1381,33 @@ struct HigherMassResonances { } } if (passKs.size() == 2) { - lResonance_gen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e()); - lResonance_gen1 = daughter1 + daughter2; + lResonanceGen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e()); + lResonanceGen1 = daughter1 + daughter2; - ROOT::Math::Boost boost{lResonance_gen.BoostToCM()}; - ROOT::Math::Boost boost1{lResonance_gen1.BoostToCM()}; + ROOT::Math::Boost boost{lResonanceGen.BoostToCM()}; + ROOT::Math::Boost boost1{lResonanceGen1.BoostToCM()}; fourVecDauCM = boost(daughter1); fourVecDauCM1 = boost1(daughter1); - auto helicity_gen = lResonance_gen.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonance_gen.Vect().Mag2())); - auto helicity_gen1 = lResonance_gen1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(lResonance_gen1.Vect().Mag2())); + auto helicityGen = lResonanceGen.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonanceGen.Vect().Mag2())); + auto helicityGen1 = lResonanceGen1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(lResonanceGen1.Vect().Mag2())); - hMChists.fill(HIST("Genf1710"), multiplicityGen, lResonance_gen.pt(), helicity_gen); - hMChists.fill(HIST("Genf1710_mass"), lResonance_gen.M()); + hMChists.fill(HIST("Genf1710"), multiplicityGen, lResonanceGen.pt(), helicityGen); + hMChists.fill(HIST("Genf1710_mass"), lResonanceGen.M()); hMChists.fill(HIST("GenRapidity"), mcParticle.y()); hMChists.fill(HIST("GenEta"), mcParticle.eta()); hMChists.fill(HIST("GenPhi"), mcParticle.phi()); - if (config.applyPairRapidityGen && std::abs(lResonance_gen1.Y()) >= 0.5) { + if (config.applyPairRapidityGen && std::abs(lResonanceGen1.Rapidity()) >= 0.5) { continue; } - hMChists.fill(HIST("Genf17102"), multiplicityGen, lResonance_gen1.pt(), helicity_gen1); - hMChists.fill(HIST("Genf1710_mass2"), lResonance_gen1.M()); - hMChists.fill(HIST("GenRapidity2"), lResonance_gen1.Y()); - hMChists.fill(HIST("GenEta2"), lResonance_gen1.Eta()); - hMChists.fill(HIST("GenPhi2"), lResonance_gen1.Phi()); + hMChists.fill(HIST("Genf17102"), multiplicityGen, lResonanceGen1.pt(), helicityGen1); + hMChists.fill(HIST("Genf1710_mass2"), lResonanceGen1.M()); + hMChists.fill(HIST("GenRapidity2"), lResonanceGen1.Rapidity()); + hMChists.fill(HIST("GenEta2"), lResonanceGen1.Eta()); + hMChists.fill(HIST("GenPhi2"), lResonanceGen1.Phi()); } passKs.clear(); // clear the vector for the next iteration } @@ -1555,7 +1554,7 @@ struct HigherMassResonances { } hMChists.fill(HIST("events_checkrec"), 18.5); - if (config.apply_rapidityMC && std::abs(mothertrack1.y()) >= 0.5) { + if (config.applyRapidityMC && std::abs(mothertrack1.y()) >= 0.5) { continue; } hMChists.fill(HIST("events_checkrec"), 19.5); @@ -1578,21 +1577,21 @@ struct HigherMassResonances { 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 helicityRec = 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())); + auto helicityRec2 = mother1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(mother1.Vect().Mag2())); - hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mothertrack1.pt(), mother1.M(), helicity_rec2); + hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mothertrack1.pt(), mother1.M(), helicityRec2); hMChists.fill(HIST("RecRapidity"), mothertrack1.y()); hMChists.fill(HIST("RecPhi"), mothertrack1.phi()); hMChists.fill(HIST("RecEta"), mothertrack1.eta()); - if (config.applyPairRapidityRec && std::abs(mother.Y()) >= 0.5) { + if (config.applyPairRapidityRec && std::abs(mother.Rapidity()) >= 0.5) { continue; } - hMChists.fill(HIST("Recf1710_pt2"), multiplicity, mother.Pt(), mother.M(), helicity_rec); - hMChists.fill(HIST("RecRapidity2"), mother.Y()); + hMChists.fill(HIST("Recf1710_pt2"), multiplicity, mother.Pt(), mother.M(), helicityRec); + hMChists.fill(HIST("RecRapidity2"), mother.Rapidity()); hMChists.fill(HIST("RecPhi2"), mother.Phi()); hMChists.fill(HIST("RecEta2"), mother.Eta()); } From 0d88ba4dae1a1ec98ff507a62fad0127a0103c8e Mon Sep 17 00:00:00 2001 From: Sawan Sawan Date: Tue, 29 Jul 2025 12:34:09 +0530 Subject: [PATCH 4/5] fixed linter error --- .../Tasks/Resonances/higherMassResonances.cxx | 71 ++-- PWGLF/Tasks/Resonances/kstarqa.cxx | 309 ++++++++++-------- 2 files changed, 205 insertions(+), 175 deletions(-) diff --git a/PWGLF/Tasks/Resonances/higherMassResonances.cxx b/PWGLF/Tasks/Resonances/higherMassResonances.cxx index 535c61e91d1..e5ccab67046 100644 --- a/PWGLF/Tasks/Resonances/higherMassResonances.cxx +++ b/PWGLF/Tasks/Resonances/higherMassResonances.cxx @@ -170,6 +170,11 @@ struct HigherMassResonances { // ConfigurableAxis axisdEdx{"axisdEdx", {20000, 0.0f, 200.0f}, "dE/dx (a.u.)"}; // ConfigurableAxis axisPtfordEbydx{"axisPtfordEbydx", {2000, 0, 20}, "pT (GeV/c)"}; // ConfigurableAxis axisMultdist{"axisMultdist", {3500, 0, 70000}, "Multiplicity distribution"}; + + // fixed variables + float rapidityMotherData = 0.5; + std::array numbers = {0, 1, 2, 3, 4, 5}; + double beamMomentum = std::sqrt(13600 * 13600 / 4 - o2::constants::physics::MassProton * o2::constants::physics::MassProton); // GeV } config; // Service PDGdatabase; @@ -183,10 +188,9 @@ struct HigherMassResonances { ROOT::Math::XYZVector randomVec, beamVec, normalVec; ROOT::Math::XYZVectorF v1CM, zaxisHE, yaxisHE, xaxisHE; // ROOT::Math::XYZVector threeVecDauCM, helicityVec, randomVec, beamVec, normalVec; - ROOT::Math::XYZVector zBeam; // ẑ: beam direction in lab frame - double beamMomentum = std::sqrt(13600 * 13600 / 4 - 0.938 * 0.938); // GeV - ROOT::Math::PxPyPzEVector beam1{0., 0., -beamMomentum, 13600. / 2.}; - ROOT::Math::PxPyPzEVector beam2{0., 0., beamMomentum, 13600. / 2.}; + ROOT::Math::XYZVector zBeam; // ẑ: beam direction in lab frame + ROOT::Math::PxPyPzEVector beam1{0., 0., -config.beamMomentum, 13600. / 2.}; + ROOT::Math::PxPyPzEVector beam2{0., 0., config.beamMomentum, 13600. / 2.}; ROOT::Math::XYZVectorF beam1CM, beam2CM; // const double massK0s = o2::constants::physics::MassK0Short; @@ -388,8 +392,8 @@ struct HigherMassResonances { 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; - float highmasscutks0 = 0.497 + config.cWidthKs0 * config.cSigmaMassKs0; + float lowmasscutks0 = o2::constants::physics::MassKPlus - config.cWidthKs0 * config.cSigmaMassKs0; + float highmasscutks0 = o2::constants::physics::MassKPlus + config.cWidthKs0 * config.cSigmaMassKs0; // float decayLength = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * RecoDecay::sqrtSumOfSquares(candidate.px(), candidate.py(), candidate.pz()); if (config.qAv0) { @@ -647,17 +651,18 @@ struct HigherMassResonances { // CosThetaHE = zaxisHE.Dot(v_CM); auto anglePhi = std::atan2(yaxisHE.Dot(v1CM), xaxisHE.Dot(v1CM)); - if (anglePhi < 0) { - anglePhi += 2 * o2::constants::math::PI; // ensure phi is in [0, 2pi] - } + anglePhi = RecoDecay::constrainAngle(anglePhi, 0.0); + // if (anglePhi < 0) { + // anglePhi += o2::constants::math::TwoPI; // ensure phi is in [0, 2pi] + // } - // if (std::abs(mother.Rapidity()) < 0.5) { + // if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { if (config.activateTHnSparseCosThStarHelicity) { // helicityVec = mother.Vect(); // 3 vector of mother in COM frame // auto cosThetaStarHelicity = helicityVec.Dot(threeVecDauCM) / (std::sqrt(threeVecDauCM.Mag2()) * std::sqrt(helicityVec.Mag2())); auto cosThetaStarHelicity = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2())); if (!isMix) { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi); } @@ -672,11 +677,11 @@ struct HigherMassResonances { daughterRotCM = boost2(daughterRot); auto cosThetaStarHelicityRot = motherRot.Vect().Dot(daughterRotCM.Vect()) / (std::sqrt(daughterRotCM.Vect().Mag2()) * std::sqrt(motherRot.Vect().Mag2())); - if (motherRot.Rapidity() < 0.5) + if (motherRot.Rapidity() < config.rapidityMotherData) hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, anglePhi); } } else { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi); } } @@ -684,18 +689,18 @@ struct HigherMassResonances { normalVec = ROOT::Math::XYZVector(mother.Py(), -mother.Px(), 0.f); auto cosThetaStarProduction = normalVec.Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(normalVec.Mag2())); if (!isMix) { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi); } for (int i = 0; i < config.cRotations; i++) { theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); 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()); - if (std::abs(motherRot.Rapidity()) < 0.5) { + if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarProduction, anglePhi); } } } else { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi); } } @@ -703,18 +708,18 @@ struct HigherMassResonances { beamVec = ROOT::Math::XYZVector(0.f, 0.f, 1.f); auto cosThetaStarBeam = beamVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2()); if (!isMix) { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi); } for (int i = 0; i < config.cRotations; i++) { theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); 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()); - if (std::abs(motherRot.Rapidity()) < 0.5) { + if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam, anglePhi); } } } else { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi); } } @@ -725,18 +730,18 @@ struct HigherMassResonances { randomVec = ROOT::Math::XYZVector(std::sin(thetaRandom) * std::cos(phiRandom), std::sin(thetaRandom) * std::sin(phiRandom), std::cos(thetaRandom)); auto cosThetaStarRandom = randomVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2()); if (!isMix) { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi); } for (int i = 0; i < config.cRotations; i++) { theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); 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()); - if (std::abs(motherRot.Rapidity()) < 0.5) { + if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, anglePhi); } } } else { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi); } } @@ -854,7 +859,7 @@ struct HigherMassResonances { } int sizeofv0indexes = v0indexes.size(); rKzeroShort.fill(HIST("NksProduced"), sizeofv0indexes); - if (config.selectTWOKsOnly && sizeofv0indexes == 2 && allConditionsMet) { + if (config.selectTWOKsOnly && sizeofv0indexes == config.numbers[2] && allConditionsMet) { fillInvMass(mother, multiplicity, daughter1, daughter2, false); } v0indexes.clear(); @@ -975,7 +980,7 @@ struct HigherMassResonances { } int sizeofv0indexes = v0indexes.size(); rKzeroShort.fill(HIST("NksProduced"), sizeofv0indexes); - if (config.selectTWOKsOnly && sizeofv0indexes == 2 && allConditionsMet) { + if (config.selectTWOKsOnly && sizeofv0indexes == config.numbers[2] && allConditionsMet) { fillInvMass(mother, multiplicity, daughter1, daughter2, false); } v0indexes.clear(); @@ -1348,7 +1353,7 @@ struct HigherMassResonances { } hMChists.fill(HIST("events_check"), 5.5); - if (config.applyRapidityMC && std::abs(mcParticle.y()) >= 0.5) { + if (config.applyRapidityMC && std::abs(mcParticle.y()) >= config.rapidityMotherData) { continue; } hMChists.fill(HIST("events_check"), 6.5); @@ -1358,7 +1363,7 @@ struct HigherMassResonances { // counter++; auto kDaughters = mcParticle.daughters_as(); - if (kDaughters.size() != 2) { + if (kDaughters.size() != config.numbers[2]) { continue; } hMChists.fill(HIST("events_check"), 7.5); @@ -1370,17 +1375,17 @@ struct HigherMassResonances { continue; } hMChists.fill(HIST("events_check"), 8.5); - if (std::abs(kCurrentDaughter.pdgCode()) == 310) { + if (std::abs(kCurrentDaughter.pdgCode()) == PDG_t::kK0Short) { passKs.push_back(true); hMChists.fill(HIST("events_check"), 9.5); if (passKs.size() == 1) { daughter1 = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassK0Short); - } else if (passKs.size() == 2) { + } else if (passKs.size() == config.numbers[2]) { daughter2 = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassK0Short); } } } - if (passKs.size() == 2) { + if (passKs.size() == config.numbers[2]) { lResonanceGen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e()); lResonanceGen1 = daughter1 + daughter2; @@ -1399,7 +1404,7 @@ struct HigherMassResonances { hMChists.fill(HIST("GenEta"), mcParticle.eta()); hMChists.fill(HIST("GenPhi"), mcParticle.phi()); - if (config.applyPairRapidityGen && std::abs(lResonanceGen1.Rapidity()) >= 0.5) { + if (config.applyPairRapidityGen && std::abs(lResonanceGen1.Rapidity()) >= config.rapidityMotherData) { continue; } @@ -1507,7 +1512,7 @@ struct HigherMassResonances { int trackv0PDG1 = std::abs(mctrackv01.pdgCode()); int trackv0PDG2 = std::abs(mctrackv02.pdgCode()); - if (std::abs(trackv0PDG1) != 310 || std::abs(trackv0PDG2) != 310) { + if (std::abs(trackv0PDG1) != PDG_t::kK0Short || std::abs(trackv0PDG2) != PDG_t::kK0Short) { continue; } hMChists.fill(HIST("events_checkrec"), 12.5); @@ -1554,7 +1559,7 @@ struct HigherMassResonances { } hMChists.fill(HIST("events_checkrec"), 18.5); - if (config.applyRapidityMC && std::abs(mothertrack1.y()) >= 0.5) { + if (config.applyRapidityMC && std::abs(mothertrack1.y()) >= config.rapidityMotherData) { continue; } hMChists.fill(HIST("events_checkrec"), 19.5); @@ -1586,7 +1591,7 @@ struct HigherMassResonances { hMChists.fill(HIST("RecPhi"), mothertrack1.phi()); hMChists.fill(HIST("RecEta"), mothertrack1.eta()); - if (config.applyPairRapidityRec && std::abs(mother.Rapidity()) >= 0.5) { + if (config.applyPairRapidityRec && std::abs(mother.Rapidity()) >= config.rapidityMotherData) { continue; } diff --git a/PWGLF/Tasks/Resonances/kstarqa.cxx b/PWGLF/Tasks/Resonances/kstarqa.cxx index 7c01c7b63ec..c8be78a1e2f 100644 --- a/PWGLF/Tasks/Resonances/kstarqa.cxx +++ b/PWGLF/Tasks/Resonances/kstarqa.cxx @@ -102,6 +102,12 @@ struct Kstarqa { Configurable cfgRCRFC{"cfgRCRFC", 0.8f, "Crossed Rows to Findable Clusters"}; Configurable cfgITSChi2NCl{"cfgITSChi2NCl", 36.0, "ITS Chi2/NCl"}; Configurable cfgTPCChi2NCl{"cfgTPCChi2NCl", 4.0, "TPC Chi2/NCl"}; + + // Other fixed variables + float lowpTcutinpTdepPID = 0.5; + std::array numbers = {0, 1, 2, 3, 4, 5}; + float rapidityMotherData = 0.5; + } selectionConfig; // Histograms are defined with HistogramRegistry @@ -175,8 +181,8 @@ struct Kstarqa { rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}}); rEventSelection.add("hMultiplicity", "Multiplicity percentile", kTH1F, {{110, 0, 110}}); - rEventSelection.add("hEventCutFlow", "No. of event after cuts", kTH1I, {{20, 0, 20}}); - std::shared_ptr hCutFlow = rEventSelection.get(HIST("hEventCutFlow")); + rEventSelection.add("hEventCut", "No. of event after cuts", kTH1I, {{20, 0, 20}}); + std::shared_ptr hCutFlow = rEventSelection.get(HIST("hEventCut")); hCutFlow->GetXaxis()->SetBinLabel(1, "All Events"); hCutFlow->GetXaxis()->SetBinLabel(2, "|Vz| < cut"); hCutFlow->GetXaxis()->SetBinLabel(3, "sel8"); @@ -188,6 +194,7 @@ struct Kstarqa { hCutFlow->GetXaxis()->SetBinLabel(9, "rctChecker"); hCutFlow->GetXaxis()->SetBinLabel(10, "kIsTriggerTVX"); hCutFlow->GetXaxis()->SetBinLabel(11, "kIsGoodZvtxFT0vsPV"); + hCutFlow->GetXaxis()->SetBinLabel(12, "IsINELgt0"); // for primary tracksbinsMultPlot if (cQAplots) { @@ -270,11 +277,38 @@ struct Kstarqa { hInvMass.add("MCcorrections/hImpactParameterRec", "Impact parameter in reconstructed MC", kTH1F, {{impactParAxis}}); hInvMass.add("MCcorrections/hImpactParameterGen", "Impact parameter in generated MC", kTH1F, {{impactParAxis}}); hInvMass.add("MCcorrections/hImpactParametervsMultiplicity", "Impact parameter vs multiplicity in reconstructed MC", kTH1F, {{impactParAxis}, {multiplicityAxis}}); - rEventSelection.add("events_check_data", "No. of events in the data", kTH1I, {{20, 0, 20}}); - rEventSelection.add("events_check", "No. of events in the generated MC", kTH1I, {{20, 0, 20}}); - rEventSelection.add("events_checkrec", "No. of events in the reconstructed MC", kTH1I, {{20, 0, 20}}); + rEventSelection.add("tracksCheckData", "No. of events in the data", kTH1I, {{10, 0, 10}}); + rEventSelection.add("eventsCheckGen", "No. of events in the generated MC", kTH1I, {{10, 0, 10}}); + rEventSelection.add("recMCparticles", "No. of events in the reconstructed MC", kTH1I, {{20, 0, 20}}); rEventSelection.add("hOccupancy", "Occupancy distribution", kTH1F, {{1000, 0, 15000}}); + std::shared_ptr hrecLabel = rEventSelection.get(HIST("hEventCut")); + hrecLabel->GetXaxis()->SetBinLabel(1, "All tracks"); + hrecLabel->GetXaxis()->SetBinLabel(2, "Track selection"); + hrecLabel->GetXaxis()->SetBinLabel(3, "has_MC"); + hrecLabel->GetXaxis()->SetBinLabel(4, "StrictlyUpperIndex"); + hrecLabel->GetXaxis()->SetBinLabel(5, "Unlike Sign"); + hrecLabel->GetXaxis()->SetBinLabel(6, "Physical Primary"); + hrecLabel->GetXaxis()->SetBinLabel(7, "PID Cut"); + hrecLabel->GetXaxis()->SetBinLabel(8, "Same mother"); + hrecLabel->GetXaxis()->SetBinLabel(9, "Generator"); + hrecLabel->GetXaxis()->SetBinLabel(10, "Rapidity"); + hrecLabel->GetXaxis()->SetBinLabel(11, "MotherPID313"); + hrecLabel->GetXaxis()->SetBinLabel(12, "Split track"); + + std::shared_ptr hDataTracks = rEventSelection.get(HIST("tracksCheckData")); + hDataTracks->GetXaxis()->SetBinLabel(1, "All tracks"); + hDataTracks->GetXaxis()->SetBinLabel(2, "Track selection"); + hDataTracks->GetXaxis()->SetBinLabel(3, "PID Cut"); + hDataTracks->GetXaxis()->SetBinLabel(4, "RmFakeTracks"); + hDataTracks->GetXaxis()->SetBinLabel(5, "Global Index"); + + std::shared_ptr hGenTracks = rEventSelection.get(HIST("eventsCheckGen")); + hGenTracks->GetXaxis()->SetBinLabel(1, "All events"); + hGenTracks->GetXaxis()->SetBinLabel(2, "INELgt0+vtz"); + hGenTracks->GetXaxis()->SetBinLabel(3, "INELgt0"); + hGenTracks->GetXaxis()->SetBinLabel(4, "Event Reconstructed"); + // Multplicity distribution if (cQAevents) { rEventSelection.add("multdist_FT0M", "FT0M Multiplicity distribution", kTH1F, {axisMultdist}); @@ -294,37 +328,37 @@ struct Kstarqa { bool selectionEvent(const Coll& collision, bool fillHist = true) { if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 0); + rEventSelection.fill(HIST("hEventCut"), 0); if (std::abs(collision.posZ()) > selectionConfig.cutzvertex) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 1); + rEventSelection.fill(HIST("hEventCut"), 1); if (!collision.sel8()) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 2); + rEventSelection.fill(HIST("hEventCut"), 2); if (selectionConfig.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 3); + rEventSelection.fill(HIST("hEventCut"), 3); if (selectionConfig.isNoITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 4); + rEventSelection.fill(HIST("hEventCut"), 4); if (selectionConfig.isNoSameBunchPileup && (!collision.selection_bit(aod::evsel::kNoSameBunchPileup))) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 5); + rEventSelection.fill(HIST("hEventCut"), 5); if (selectionConfig.isAllLayersGoodITS && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 6); + rEventSelection.fill(HIST("hEventCut"), 6); if (selectionConfig.isApplyOccCut && (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard))) return false; @@ -332,28 +366,28 @@ struct Kstarqa { if (selectionConfig.isApplyOccCut && (std::abs(collision.trackOccupancyInTimeRange()) > selectionConfig.configOccCut)) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 7); + rEventSelection.fill(HIST("hEventCut"), 7); if (rctCut.requireRCTFlagChecker && !rctChecker(collision)) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 8); + rEventSelection.fill(HIST("hEventCut"), 8); if (selectionConfig.isTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 9); + rEventSelection.fill(HIST("hEventCut"), 9); if (selectionConfig.isGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) return false; if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 10); + rEventSelection.fill(HIST("hEventCut"), 10); if (selectionConfig.isINELgt0 && !collision.isInelGt0()) { return false; } if (fillHist) - rEventSelection.fill(HIST("hEventCutFlow"), 11); + rEventSelection.fill(HIST("hEventCut"), 11); return true; } @@ -464,23 +498,23 @@ struct Kstarqa { bool selectionPIDNew(const T& candidate, int PID) { if (PID == 0) { - if (candidate.pt() < 0.5 && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi) { + if (candidate.pt() < selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi) { return true; } - if (candidate.pt() >= 0.5 && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi && candidate.hasTOF() && std::abs(candidate.tofNSigmaPi()) < nsigmaCutTOFPi) { + if (candidate.pt() >= selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi && candidate.hasTOF() && std::abs(candidate.tofNSigmaPi()) < nsigmaCutTOFPi) { return true; } - if (candidate.pt() >= 0.5 && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi && !candidate.hasTOF()) { + if (candidate.pt() >= selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi && !candidate.hasTOF()) { return true; } } else if (PID == 1) { - if (candidate.pt() < 0.5 && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa) { + if (candidate.pt() < selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa) { return true; } - if (candidate.pt() >= 0.5 && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa && candidate.hasTOF() && std::abs(candidate.tofNSigmaKa()) < nsigmaCutTOFKa) { + if (candidate.pt() >= selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa && candidate.hasTOF() && std::abs(candidate.tofNSigmaKa()) < nsigmaCutTOFKa) { return true; } - if (candidate.pt() >= 0.5 && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa && !candidate.hasTOF()) { + if (candidate.pt() >= selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa && !candidate.hasTOF()) { return true; } } @@ -586,6 +620,7 @@ struct Kstarqa { using EventCandidatesMC = soa::Join; using TrackCandidatesMC = soa::Filtered>; + using EventMCGenerated = soa::Join; // aod::CentNGlobals, aod::CentNTPVs, aod::CentMFTs //*********Varibles declaration*************** float multiplicity{-1.0}, theta2; @@ -600,13 +635,13 @@ struct Kstarqa { ROOT::Math::Boost boost{mother.BoostToCM()}; // boost mother to center of mass frame fourVecDauCM = boost(daughterSelected); // boost the frame of daughter same as mother - // if (std::abs(mother.Rapidity()) < 0.5) { + // if (std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { if (activateTHnSparseCosThStarHelicity) { auto cosThetaStarHelicity = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2())); if (track1.sign() * track2.sign() < 0) { if (!isMix) { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { hInvMass.fill(HIST("h3KstarInvMassUnlikeSign"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity); } @@ -622,15 +657,15 @@ struct Kstarqa { auto cosThetaStarHelicityRot = motherRot.Vect().Dot(daughterRotCM.Vect()) / (std::sqrt(daughterRotCM.Vect().Mag2()) * std::sqrt(motherRot.Vect().Mag2())); - if (calcRotational && motherRot.Rapidity() < 0.5) + if (calcRotational && motherRot.Rapidity() < selectionConfig.rapidityMotherData) hInvMass.fill(HIST("h3KstarInvMassRotated"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot); } - } else if (isMix && std::abs(mother.Rapidity()) < 0.5) { + } else if (isMix && std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { hInvMass.fill(HIST("h3KstarInvMassMixed"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity); } } else { if (!isMix) { - if (calcLikeSign && std::abs(mother.Rapidity()) < 0.5) { + if (calcLikeSign && std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { if (track1.sign() > 0 && track2.sign() > 0) { hInvMass.fill(HIST("h3KstarInvMasslikeSignPP"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity); } else if (track1.sign() < 0 && track2.sign() < 0) { @@ -646,7 +681,7 @@ struct Kstarqa { if (track1.sign() * track2.sign() < 0) { if (!isMix) { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { hInvMass.fill(HIST("h3KstarInvMassUnlikeSign"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction); } for (int i = 0; i < cRotations; i++) { @@ -654,15 +689,15 @@ struct Kstarqa { daughterRot = ROOT::Math::PxPyPzMVector(daughter1.Px() * std::cos(theta2) - daughter1.Py() * std::sin(theta2), daughter1.Px() * std::sin(theta2) + daughter1.Py() * std::cos(theta2), daughter1.Pz(), daughter1.M()); motherRot = daughterRot + daughter2; - if (calcRotational && std::abs(motherRot.Rapidity()) < 0.5) + if (calcRotational && std::abs(motherRot.Rapidity()) < selectionConfig.rapidityMotherData) hInvMass.fill(HIST("h3KstarInvMassRotated"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarProduction); } - } else if (isMix && std::abs(mother.Rapidity()) < 0.5) { + } else if (isMix && std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { hInvMass.fill(HIST("h3KstarInvMassMixed"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction); } } else { if (!isMix) { - if (calcLikeSign && std::abs(mother.Rapidity()) < 0.5) { + if (calcLikeSign && std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { if (track1.sign() > 0 && track2.sign() > 0) { hInvMass.fill(HIST("h3KstarInvMasslikeSignPP"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction); } else if (track1.sign() < 0 && track2.sign() < 0) { @@ -677,7 +712,7 @@ struct Kstarqa { if (track1.sign() * track2.sign() < 0) { if (!isMix) { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { hInvMass.fill(HIST("h3KstarInvMassUnlikeSign"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam); } for (int i = 0; i < cRotations; i++) { @@ -685,14 +720,14 @@ struct Kstarqa { daughterRot = ROOT::Math::PxPyPzMVector(daughter1.Px() * std::cos(theta2) - daughter1.Py() * std::sin(theta2), daughter1.Px() * std::sin(theta2) + daughter1.Py() * std::cos(theta2), daughter1.Pz(), daughter1.M()); motherRot = daughterRot + daughter2; - if (calcRotational && std::abs(motherRot.Rapidity()) < 0.5) + if (calcRotational && std::abs(motherRot.Rapidity()) < selectionConfig.rapidityMotherData) hInvMass.fill(HIST("h3KstarInvMassRotated"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam); } - } else if (isMix && std::abs(mother.Rapidity()) < 0.5) { + } else if (isMix && std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { hInvMass.fill(HIST("h3KstarInvMassMixed"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam); } } else { - if (calcLikeSign && std::abs(mother.Rapidity()) < 0.5) { + if (calcLikeSign && std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { if (track1.sign() > 0 && track2.sign() > 0) { hInvMass.fill(HIST("h3KstarInvMasslikeSignPP"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam); } else if (track1.sign() < 0 && track2.sign() < 0) { @@ -709,7 +744,7 @@ struct Kstarqa { if (track1.sign() * track2.sign() < 0) { if (!isMix) { - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { hInvMass.fill(HIST("h3KstarInvMassUnlikeSign"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom); } for (int i = 0; i < cRotations; i++) { @@ -717,15 +752,15 @@ struct Kstarqa { daughterRot = ROOT::Math::PxPyPzMVector(daughter1.Px() * std::cos(theta2) - daughter1.Py() * std::sin(theta2), daughter1.Px() * std::sin(theta2) + daughter1.Py() * std::cos(theta2), daughter1.Pz(), daughter1.M()); motherRot = daughterRot + daughter2; - if (calcRotational && std::abs(motherRot.Rapidity()) < 0.5) + if (calcRotational && std::abs(motherRot.Rapidity()) < selectionConfig.rapidityMotherData) hInvMass.fill(HIST("h3KstarInvMassRotated"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom); } - } else if (isMix && std::abs(mother.Rapidity()) < 0.5) { + } else if (isMix && std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { hInvMass.fill(HIST("h3KstarInvMassMixed"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom); } } else { if (!isMix) { - if (calcLikeSign && std::abs(mother.Rapidity()) < 0.5) { + if (calcLikeSign && std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { if (track1.sign() > 0 && track2.sign() > 0) { hInvMass.fill(HIST("h3KstarInvMasslikeSignPP"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom); } else if (track1.sign() < 0 && track2.sign() < 0) { @@ -739,16 +774,13 @@ struct Kstarqa { } // int counter = 0; - void processSE(EventCandidates::iterator const& collision, TrackCandidates const& tracks, aod::BCs const&) { - rEventSelection.fill(HIST("events_check_data"), 0.5); // if (cTVXEvsel && (!collision.selection_bit(aod::evsel::kIsTriggerTVX))) { // return; // } - // rEventSelection.fill(HIST("events_check_data"), 1.5); // if (timFrameEvsel && (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(aod::evsel::kNoITSROFrameBorder))) { // return; @@ -757,25 +789,23 @@ struct Kstarqa { // if (!collision.sel8()) { // return; // } - // rEventSelection.fill(HIST("events_check_data"), 2.5); if (!selectionEvent(collision, true)) { return; } - rEventSelection.fill(HIST("events_check_data"), 3.5); int occupancy = collision.trackOccupancyInTimeRange(); rEventSelection.fill(HIST("hOccupancy"), occupancy); multiplicity = -1; - if (cSelectMultEstimator == 0) { + if (cSelectMultEstimator == selectionConfig.numbers[0]) { // FT0M multiplicity = collision.centFT0M(); - } else if (cSelectMultEstimator == 1) { + } else if (cSelectMultEstimator == selectionConfig.numbers[1]) { multiplicity = collision.centFT0A(); - } else if (cSelectMultEstimator == 2) { + } else if (cSelectMultEstimator == selectionConfig.numbers[2]) { multiplicity = collision.centFT0C(); - } else if (cSelectMultEstimator == 3) { + } else if (cSelectMultEstimator == selectionConfig.numbers[3]) { multiplicity = collision.centFV0A(); } else { multiplicity = collision.centFT0M(); @@ -801,13 +831,14 @@ struct Kstarqa { } for (const auto& [track1, track2] : combinations(CombinationsFullIndexPolicy(tracks, tracks))) { + rEventSelection.fill(HIST("tracksCheckData"), 0.5); if (!selectionTrack(track1)) { continue; } if (!selectionTrack(track2)) { continue; } - rEventSelection.fill(HIST("events_check_data"), 4.5); + rEventSelection.fill(HIST("tracksCheckData"), 1.5); if (cQAplots) { hPID.fill(HIST("Before/hNsigmaTPC_Ka_before"), track1.pt(), track1.tpcNSigmaKa()); @@ -839,11 +870,6 @@ struct Kstarqa { } } - rEventSelection.fill(HIST("events_check_data"), 5.5); - // if (counter < 1e4) - // std::cout << "TOF beta value is " << track1.beta() << std::endl; - // counter++; - if (cQAevents) { rEventSelection.fill(HIST("hDcaxy"), track1.dcaXY()); rEventSelection.fill(HIST("hDcaz"), track1.dcaZ()); @@ -860,7 +886,7 @@ struct Kstarqa { if (applypTdepPID && !selectionPIDNew(track2, 0)) // Track 2 is checked with Pion continue; - rEventSelection.fill(HIST("events_check_data"), 6.5); + rEventSelection.fill(HIST("tracksCheckData"), 2.5); if (cFakeTrack && isFakeTrack(track1, 1)) // Kaon continue; @@ -876,7 +902,7 @@ struct Kstarqa { // continue; // } - rEventSelection.fill(HIST("events_check_data"), 7.5); + rEventSelection.fill(HIST("tracksCheckData"), 3.5); if (cQAplots) { hPID.fill(HIST("After/hDcaxyPi"), track2.dcaXY()); @@ -901,7 +927,7 @@ struct Kstarqa { if (track1.globalIndex() == track2.globalIndex()) continue; - rEventSelection.fill(HIST("events_check_data"), 8.5); + rEventSelection.fill(HIST("tracksCheckData"), 4.5); daughter1 = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massKa); daughter2 = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massPi); @@ -964,7 +990,7 @@ struct Kstarqa { isMix = true; - if (std::abs(mother.Rapidity()) < 0.5) { + if (std::abs(mother.Rapidity()) < selectionConfig.rapidityMotherData) { fillInvMass(daughter1, daughter2, mother, multiplicity, isMix, t1, t2); } } @@ -972,25 +998,22 @@ struct Kstarqa { }; // Call mixing based on selected estimator - if (cSelectMultEstimator == 0) { + if (cSelectMultEstimator == selectionConfig.numbers[0]) { // FT0M runMixing(pair1, [](const auto& c) { return c.centFT0M(); }); - } else if (cSelectMultEstimator == 1) { + } else if (cSelectMultEstimator == selectionConfig.numbers[1]) { runMixing(pair2, [](const auto& c) { return c.centFT0A(); }); - } else if (cSelectMultEstimator == 2) { + } else if (cSelectMultEstimator == selectionConfig.numbers[2]) { runMixing(pair3, [](const auto& c) { return c.centFT0C(); }); - } else if (cSelectMultEstimator == 3) { + } else if (cSelectMultEstimator == selectionConfig.numbers[3]) { runMixing(pair4, [](const auto& c) { return c.centFV0A(); }); } } PROCESS_SWITCH(Kstarqa, processME, "Process Mixed event", true); - void processGen(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& collisions) + void processGen(EventMCGenerated::iterator const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& collisions) { - rEventSelection.fill(HIST("events_check"), 0.5); - if (std::abs(mcCollision.posZ()) < selectionConfig.cutzvertex) { - rEventSelection.fill(HIST("events_check"), 1.5); - } + rEventSelection.fill(HIST("eventsCheckGen"), 0.5); int nChInel = 0; for (const auto& mcParticle : mcParticles) { @@ -1002,51 +1025,56 @@ struct Kstarqa { } } if (nChInel > 0 && std::abs(mcCollision.posZ()) < selectionConfig.cutzvertex) - rEventSelection.fill(HIST("events_check"), 2.5); + rEventSelection.fill(HIST("eventsCheckGen"), 1.5); std::vector selectedEvents(collisions.size()); int nevts = 0; multiplicity = -1.0; // float impactParameter = mcCollision.impactParameter(); - // if (mcCollision.isInelGt0()) { - // return; - // } + if (selectionConfig.isINELgt0 && !mcCollision.isInelGt0()) { + return; + } + rEventSelection.fill(HIST("eventsCheckGen"), 2.5); for (const auto& collision : collisions) { // if (!collision.sel8() || std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex) { - if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex) { - continue; - } - if (!collision.sel8()) { - continue; - } - if (selectionConfig.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { - continue; - } - if (selectionConfig.isTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) { - continue; - } - if (selectionConfig.isINELgt0 && !collision.isInelGt0()) { - continue; - } - if (selectionConfig.isNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { - continue; - } - if (selectionConfig.isGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + // if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex) { + // continue; + // } + // if (!collision.sel8()) { + // continue; + // } + // if (selectionConfig.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + // continue; + // } + // if (selectionConfig.isTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) { + // continue; + // } + // if (selectionConfig.isINELgt0 && !collision.isInelGt0()) { + // continue; + // } + // if (selectionConfig.isNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { + // continue; + // } + // if (selectionConfig.isGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + // continue; + // } + if (!selectionEvent(collision, true)) { continue; } multiplicity = collision.centFT0M(); hInvMass.fill(HIST("h1GenMult"), multiplicity); - selectedEvents[nevts++] = collision.mcCollision_as().globalIndex(); + int occupancy = collision.trackOccupancyInTimeRange(); rEventSelection.fill(HIST("hOccupancy"), occupancy); + + selectedEvents[nevts++] = collision.mcCollision_as().globalIndex(); } selectedEvents.resize(nevts); - rEventSelection.fill(HIST("events_check"), 3.5); for (const auto& mcParticle : mcParticles) { - if (std::abs(mcParticle.y()) < 0.5 && std::abs(mcParticle.pdgCode()) == o2::constants::physics::kK0Star892) { + if (std::abs(mcParticle.y()) < selectionConfig.rapidityMotherData && std::abs(mcParticle.pdgCode()) == o2::constants::physics::kK0Star892) { hInvMass.fill(HIST("hAllKstarGenCollisisons"), multiplicity, mcParticle.pt()); } } @@ -1057,25 +1085,23 @@ struct Kstarqa { return; } hInvMass.fill(HIST("hAllGenCollisions1Rec"), multiplicity); - rEventSelection.fill(HIST("events_check"), 4.5); + rEventSelection.fill(HIST("eventsCheckGen"), 3.5); for (const auto& mcParticle : mcParticles) { - if (std::abs(mcParticle.y()) >= 0.5) { + + if (std::abs(mcParticle.y()) >= selectionConfig.rapidityMotherData) { continue; } - rEventSelection.fill(HIST("events_check"), 5.5); if (std::abs(mcParticle.pdgCode()) != o2::constants::physics::kK0Star892) { continue; } - rEventSelection.fill(HIST("events_check"), 6.5); hInvMass.fill(HIST("hAllKstarGenCollisisons1Rec"), multiplicity, mcParticle.pt()); auto kDaughters = mcParticle.daughters_as(); - if (kDaughters.size() != 2) { + if (kDaughters.size() != selectionConfig.numbers[2]) { continue; } - rEventSelection.fill(HIST("events_check"), 7.5); auto passkaon = false; auto passpion = false; @@ -1083,7 +1109,6 @@ struct Kstarqa { if (!kCurrentDaughter.isPhysicalPrimary()) { continue; } - rEventSelection.fill(HIST("events_check"), 8.5); if (std::abs(kCurrentDaughter.pdgCode()) == PDG_t::kKPlus) { passkaon = true; @@ -1104,8 +1129,11 @@ struct Kstarqa { } PROCESS_SWITCH(Kstarqa, processGen, "Process Generated", false); - void processEvtLossSigLossMC(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& recCollisions) + void processEvtLossSigLossMC(EventMCGenerated::iterator const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& recCollisions) { + if (selectionConfig.isINELgt0 && !mcCollision.isInelGt0()) { + return; + } auto impactPar = mcCollision.impactParameter(); hInvMass.fill(HIST("MCcorrections/hImpactParameterGen"), impactPar); @@ -1126,7 +1154,7 @@ struct Kstarqa { // Generated MC for (const auto& mcPart : mcParticles) { - if (std::abs(mcPart.y()) >= 0.5 || std::abs(mcPart.pdgCode()) != o2::constants::physics::kK0Star892) + if (std::abs(mcPart.y()) >= selectionConfig.rapidityMotherData || std::abs(mcPart.pdgCode()) != o2::constants::physics::kK0Star892) continue; // signal loss estimation @@ -1141,12 +1169,9 @@ struct Kstarqa { void processRec(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&, aod::McCollisions const& /*mcCollisions*/) { - rEventSelection.fill(HIST("events_checkrec"), 0.5); - if (!collision.has_mcCollision()) { return; } - rEventSelection.fill(HIST("events_checkrec"), 1.5); if (selectionConfig.isINELgt0 && !collision.isInelGt0()) { return; @@ -1154,32 +1179,33 @@ struct Kstarqa { multiplicity = collision.centFT0M(); hInvMass.fill(HIST("hAllRecCollisions"), multiplicity); - // if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex || !collision.sel8()) { - if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex) { + if (!selectionEvent(collision, false)) { return; } - rEventSelection.fill(HIST("events_checkrec"), 2.5); - if (selectionConfig.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { - return; - } - rEventSelection.fill(HIST("events_checkrec"), 3.5); + // // if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex || !collision.sel8()) { + // if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex) { + // return; + // } - if (selectionConfig.isTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) { - return; - } - rEventSelection.fill(HIST("events_checkrec"), 4.5); + // if (selectionConfig.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + // return; + // } - if (!collision.sel8()) { - return; - } + // if (selectionConfig.isTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) { + // return; + // } - if (selectionConfig.isNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { - return; - } - if (selectionConfig.isGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { - return; - } + // if (!collision.sel8()) { + // return; + // } + + // if (selectionConfig.isNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { + // return; + // } + // if (selectionConfig.isGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + // return; + // } multiplicity = collision.centFT0M(); hInvMass.fill(HIST("h1RecMult"), multiplicity); @@ -1189,35 +1215,34 @@ struct Kstarqa { if (!selectionTrack(track1)) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 5.5); if (!track1.has_mcParticle()) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 6.5); auto track1ID = track1.index(); for (const auto& track2 : tracks) { + rEventSelection.fill(HIST("recMCparticles"), 0.5); if (!track2.has_mcParticle()) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 7.5); + rEventSelection.fill(HIST("recMCparticles"), 1.5); if (!selectionTrack(track2)) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 8.5); + rEventSelection.fill(HIST("recMCparticles"), 2.5); auto track2ID = track2.index(); if (track2ID <= track1ID) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 9.5); + rEventSelection.fill(HIST("recMCparticles"), 3.5); if (track1.sign() * track2.sign() >= 0) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 10.5); + rEventSelection.fill(HIST("recMCparticles"), 4.5); const auto mctrack1 = track1.mcParticle(); const auto mctrack2 = track2.mcParticle(); @@ -1256,12 +1281,11 @@ struct Kstarqa { if (!mctrack1.isPhysicalPrimary()) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 11.5); if (!mctrack2.isPhysicalPrimary()) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 12.5); + rEventSelection.fill(HIST("recMCparticles"), 5.5); // if (!(track1PDG == PDG_t::kKPlus && track2PDG == PDG_t::kPiPlus)) { // continue; @@ -1272,34 +1296,33 @@ struct Kstarqa { if ((track2PDG != PDG_t::kPiPlus) && (track2PDG != PDG_t::kKPlus)) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 13.5); - rEventSelection.fill(HIST("events_checkrec"), 14.5); + rEventSelection.fill(HIST("recMCparticles"), 6.5); for (const auto& mothertrack1 : mctrack1.mothers_as()) { for (const auto& mothertrack2 : mctrack2.mothers_as()) { if (mothertrack1.pdgCode() != mothertrack2.pdgCode()) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 15.5); if (mothertrack1.globalIndex() != mothertrack2.globalIndex()) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 16.5); + rEventSelection.fill(HIST("recMCparticles"), 7.5); if (!mothertrack1.producedByGenerator()) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 17.5); + rEventSelection.fill(HIST("recMCparticles"), 8.5); - if (std::abs(mothertrack1.y()) >= 0.5) { + if (std::abs(mothertrack1.y()) >= selectionConfig.rapidityMotherData) { continue; } - rEventSelection.fill(HIST("events_checkrec"), 18.5); + rEventSelection.fill(HIST("recMCparticles"), 9.5); if (std::abs(mothertrack1.pdgCode()) != o2::constants::physics::kK0Star892) { continue; } + rEventSelection.fill(HIST("recMCparticles"), 10.5); if (track1PDG == PDG_t::kPiPlus) { if (!applypTdepPID && !(selectionPID(track1, 0) && selectionPID(track2, 1))) { // pion and kaon @@ -1319,6 +1342,8 @@ struct Kstarqa { hInvMass.fill(HIST("h1KSRecsplit"), mothertrack1.pt()); continue; } + rEventSelection.fill(HIST("recMCparticles"), 11.5); + oldindex = mothertrack1.globalIndex(); if (track1.sign() * track2.sign() < 0) { daughter1 = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massKa); @@ -1327,7 +1352,7 @@ struct Kstarqa { hInvMass.fill(HIST("h2KstarRecpt2"), mothertrack1.pt(), multiplicity, std::sqrt(mothertrack1.e() * mothertrack1.e() - mothertrack1.p() * mothertrack1.p())); - if (applyRecMotherRapidity && mother.Rapidity() >= 0.5) { + if (applyRecMotherRapidity && mother.Rapidity() >= selectionConfig.rapidityMotherData) { continue; } From d2b8df48026272322144bae08853ea8373c24387 Mon Sep 17 00:00:00 2001 From: Sawan Sawan Date: Tue, 29 Jul 2025 14:49:17 +0530 Subject: [PATCH 5/5] optimising code --- .../Tasks/Resonances/higherMassResonances.cxx | 15 +++--- PWGLF/Tasks/Resonances/kstarqa.cxx | 47 ++++++++++++------- 2 files changed, 37 insertions(+), 25 deletions(-) diff --git a/PWGLF/Tasks/Resonances/higherMassResonances.cxx b/PWGLF/Tasks/Resonances/higherMassResonances.cxx index e5ccab67046..3a34c13fad2 100644 --- a/PWGLF/Tasks/Resonances/higherMassResonances.cxx +++ b/PWGLF/Tasks/Resonances/higherMassResonances.cxx @@ -173,8 +173,9 @@ struct HigherMassResonances { // fixed variables float rapidityMotherData = 0.5; - std::array numbers = {0, 1, 2, 3, 4, 5}; - double beamMomentum = std::sqrt(13600 * 13600 / 4 - o2::constants::physics::MassProton * o2::constants::physics::MassProton); // GeV + float beamEnergy = 13600.0; + double beamMomentum = std::sqrt(beamEnergy * beamEnergy / 4 - o2::constants::physics::MassProton * o2::constants::physics::MassProton); // GeV + int noOfDaughters = 2; } config; // Service PDGdatabase; @@ -859,7 +860,7 @@ struct HigherMassResonances { } int sizeofv0indexes = v0indexes.size(); rKzeroShort.fill(HIST("NksProduced"), sizeofv0indexes); - if (config.selectTWOKsOnly && sizeofv0indexes == config.numbers[2] && allConditionsMet) { + if (config.selectTWOKsOnly && sizeofv0indexes == config.noOfDaughters && allConditionsMet) { fillInvMass(mother, multiplicity, daughter1, daughter2, false); } v0indexes.clear(); @@ -980,7 +981,7 @@ struct HigherMassResonances { } int sizeofv0indexes = v0indexes.size(); rKzeroShort.fill(HIST("NksProduced"), sizeofv0indexes); - if (config.selectTWOKsOnly && sizeofv0indexes == config.numbers[2] && allConditionsMet) { + if (config.selectTWOKsOnly && sizeofv0indexes == config.noOfDaughters && allConditionsMet) { fillInvMass(mother, multiplicity, daughter1, daughter2, false); } v0indexes.clear(); @@ -1363,7 +1364,7 @@ struct HigherMassResonances { // counter++; auto kDaughters = mcParticle.daughters_as(); - if (kDaughters.size() != config.numbers[2]) { + if (kDaughters.size() != config.noOfDaughters) { continue; } hMChists.fill(HIST("events_check"), 7.5); @@ -1380,12 +1381,12 @@ struct HigherMassResonances { hMChists.fill(HIST("events_check"), 9.5); if (passKs.size() == 1) { daughter1 = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassK0Short); - } else if (passKs.size() == config.numbers[2]) { + } else if (static_cast(passKs.size()) == config.noOfDaughters) { daughter2 = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), o2::constants::physics::MassK0Short); } } } - if (passKs.size() == config.numbers[2]) { + if (static_cast(passKs.size()) == config.noOfDaughters) { lResonanceGen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e()); lResonanceGen1 = daughter1 + daughter2; diff --git a/PWGLF/Tasks/Resonances/kstarqa.cxx b/PWGLF/Tasks/Resonances/kstarqa.cxx index c8be78a1e2f..21f4db8e9c8 100644 --- a/PWGLF/Tasks/Resonances/kstarqa.cxx +++ b/PWGLF/Tasks/Resonances/kstarqa.cxx @@ -104,12 +104,22 @@ struct Kstarqa { Configurable cfgTPCChi2NCl{"cfgTPCChi2NCl", 4.0, "TPC Chi2/NCl"}; // Other fixed variables - float lowpTcutinpTdepPID = 0.5; - std::array numbers = {0, 1, 2, 3, 4, 5}; + float lowPtCutPID = 0.5; + int noOfDaughters = 2; float rapidityMotherData = 0.5; } selectionConfig; + enum MultEstimator { + kFT0M, + kFT0A, + kFT0C, + kFV0A, + kFV0C, + kFV0M, + kNEstimators // useful if you want to iterate or size things + }; + // Histograms are defined with HistogramRegistry HistogramRegistry rEventSelection{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry hInvMass{"hInvMass", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; @@ -498,23 +508,23 @@ struct Kstarqa { bool selectionPIDNew(const T& candidate, int PID) { if (PID == 0) { - if (candidate.pt() < selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi) { + if (candidate.pt() < selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi) { return true; } - if (candidate.pt() >= selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi && candidate.hasTOF() && std::abs(candidate.tofNSigmaPi()) < nsigmaCutTOFPi) { + if (candidate.pt() >= selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi && candidate.hasTOF() && std::abs(candidate.tofNSigmaPi()) < nsigmaCutTOFPi) { return true; } - if (candidate.pt() >= selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi && !candidate.hasTOF()) { + if (candidate.pt() >= selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaPi()) < nsigmaCutTPCPi && !candidate.hasTOF()) { return true; } } else if (PID == 1) { - if (candidate.pt() < selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa) { + if (candidate.pt() < selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa) { return true; } - if (candidate.pt() >= selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa && candidate.hasTOF() && std::abs(candidate.tofNSigmaKa()) < nsigmaCutTOFKa) { + if (candidate.pt() >= selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa && candidate.hasTOF() && std::abs(candidate.tofNSigmaKa()) < nsigmaCutTOFKa) { return true; } - if (candidate.pt() >= selectionConfig.lowpTcutinpTdepPID && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa && !candidate.hasTOF()) { + if (candidate.pt() >= selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaKa()) < nsigmaCutTPCKa && !candidate.hasTOF()) { return true; } } @@ -799,17 +809,18 @@ struct Kstarqa { multiplicity = -1; - if (cSelectMultEstimator == selectionConfig.numbers[0]) { // FT0M + if (cSelectMultEstimator == kFT0M) { multiplicity = collision.centFT0M(); - } else if (cSelectMultEstimator == selectionConfig.numbers[1]) { + } else if (cSelectMultEstimator == kFT0A) { multiplicity = collision.centFT0A(); - } else if (cSelectMultEstimator == selectionConfig.numbers[2]) { + } else if (cSelectMultEstimator == kFT0C) { multiplicity = collision.centFT0C(); - } else if (cSelectMultEstimator == selectionConfig.numbers[3]) { + } else if (cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); } else { - multiplicity = collision.centFT0M(); + multiplicity = collision.centFT0M(); // default } + /* else if (cSelectMultEstimator == 4) { multiplicity = collision.centMFT(); } */ @@ -998,13 +1009,13 @@ struct Kstarqa { }; // Call mixing based on selected estimator - if (cSelectMultEstimator == selectionConfig.numbers[0]) { // FT0M + if (cSelectMultEstimator == kFT0M) { runMixing(pair1, [](const auto& c) { return c.centFT0M(); }); - } else if (cSelectMultEstimator == selectionConfig.numbers[1]) { + } else if (cSelectMultEstimator == kFT0A) { runMixing(pair2, [](const auto& c) { return c.centFT0A(); }); - } else if (cSelectMultEstimator == selectionConfig.numbers[2]) { + } else if (cSelectMultEstimator == kFT0C) { runMixing(pair3, [](const auto& c) { return c.centFT0C(); }); - } else if (cSelectMultEstimator == selectionConfig.numbers[3]) { + } else if (cSelectMultEstimator == kFV0A) { runMixing(pair4, [](const auto& c) { return c.centFV0A(); }); } } @@ -1099,7 +1110,7 @@ struct Kstarqa { hInvMass.fill(HIST("hAllKstarGenCollisisons1Rec"), multiplicity, mcParticle.pt()); auto kDaughters = mcParticle.daughters_as(); - if (kDaughters.size() != selectionConfig.numbers[2]) { + if (kDaughters.size() != selectionConfig.noOfDaughters) { continue; }