Skip to content

[PWGLF] optimised code #12682

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 2 commits into from
Aug 22, 2025
Merged

[PWGLF] optimised code #12682

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5c5b0cb
optimisation
Aug 22, 2025
19d5489
calibrated multiplicity used
Aug 22, 2025
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
135 changes: 110 additions & 25 deletions PWGLF/Tasks/Resonances/k892analysis.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,19 +15,23 @@
///
/// \author Bong-Hwi Lim <bong-hwi.lim@cern.ch>, Sawan Sawan <sawan.sawan@cern.ch>

#include <TLorentzVector.h>
#include "TF1.h"
#include "TRandom3.h"
#include "PWGLF/DataModel/LFResonanceTables.h"
#include "PWGLF/DataModel/mcCentrality.h"
#include "PWGLF/Utils/inelGt.h"

#include "Common/DataModel/PIDResponse.h"
#include "Common/DataModel/Centrality.h"
#include "Common/DataModel/EventSelection.h"
#include "Framework/AnalysisTask.h"
#include "Common/DataModel/PIDResponse.h"

#include "CommonConstants/PhysicsConstants.h"
#include "DataFormatsParameters/GRPObject.h"
#include "Framework/ASoAHelpers.h"
#include "Framework/AnalysisTask.h"
#include "Framework/runDataProcessing.h"
#include "PWGLF/DataModel/LFResonanceTables.h"
#include "DataFormatsParameters/GRPObject.h"
#include "CommonConstants/PhysicsConstants.h"

#include "TF1.h"
#include "TRandom3.h"
#include <TLorentzVector.h>

Check failure on line 34 in PWGLF/Tasks/Resonances/k892analysis.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.

using namespace o2;
using namespace o2::framework;
Expand Down Expand Up @@ -59,6 +63,9 @@
Configurable<bool> invmass1D{"invmass1D", false, "Invariant mass 1D"};
Configurable<bool> studyAntiparticle{"studyAntiparticle", false, "Study anti-particles separately"};
Configurable<bool> fillPidPlots{"fillPidPlots", false, "Make TPC and TOF PID plots"};
Configurable<bool> cisInelGt0{"cisInelGt0", true, "check if INEL>0"};
Configurable<bool> cMCCent{"cMCCent", true, "Using calibrated MC centrality (for FT0M)"};

// Configurable<bool> applyOccupancyCut{"applyOccupancyCut", false, "Apply occupancy cut"};
// Configurable<int> occupancyCut{"occupancyCut", 1000, "Mimimum Occupancy cut"};

Expand Down Expand Up @@ -426,11 +433,11 @@
return false;
}

template <bool IsMC, bool IsMix, typename CollisionType, typename TracksType>
void fillHistograms(const CollisionType& collision, const TracksType& dTracks1, const TracksType& dTracks2)
template <bool IsMC, bool IsMix, typename CollisionType, typename TracksType, typename Multdatamc>
void fillHistograms(const CollisionType& collision, const TracksType& dTracks1, const TracksType& dTracks2, const Multdatamc& multiplicity)
{
// auto multNTracksPV = collision.multNTracksPV();
auto multiplicity = collision.cent();
// auto multiplicity = collision.cent();
if (additionalEvsel && !eventSelected(collision, multiplicity)) {
return;
}
Expand Down Expand Up @@ -463,7 +470,7 @@
histos.fill(HIST("MultCalib/GloPVpi_after"), dTracks1.size(), collision.multNTracksPV()); // global tracks vs PV tracks after the multiplicity calibration cuts
}

TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance, ldaughterRot, lresonanceRot;

Check failure on line 473 in PWGLF/Tasks/Resonances/k892analysis.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.
for (const auto& [trk1, trk2] : combinations(CombinationsFullIndexPolicy(dTracks1, dTracks2))) {

// Full index policy is needed to consider all possible combinations
Expand Down Expand Up @@ -583,7 +590,7 @@
lDecayDaughter2.SetPtEtaPhiM(trk2.pt(), trk2.eta(), trk2.phi(), massKa);
lResonance = lDecayDaughter1 + lDecayDaughter2;
// Rapidity cut
if (std::abs(lResonance.Rapidity()) >= 0.5)

Check failure on line 593 in PWGLF/Tasks/Resonances/k892analysis.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.
continue;
if (cfgCutsOnMother) {
if (lResonance.Pt() >= cMaxPtMotherCut) // excluding candidates in overflow
Expand Down Expand Up @@ -653,11 +660,11 @@

// MC
if constexpr (IsMC) {
if (std::abs(trk1.pdgCode()) != 211 || std::abs(trk2.pdgCode()) != 321)

Check failure on line 663 in PWGLF/Tasks/Resonances/k892analysis.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.

Check failure on line 663 in PWGLF/Tasks/Resonances/k892analysis.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;
if (trk1.motherId() != trk2.motherId()) // Same mother
continue;
if (std::abs(trk1.motherPDG()) != 313)

Check failure on line 667 in PWGLF/Tasks/Resonances/k892analysis.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[magic-number] 

Avoid magic numbers in expressions. Assign the value to a clearly named variable or constant.

Check failure on line 667 in PWGLF/Tasks/Resonances/k892analysis.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;

// Track selection check.
Expand Down Expand Up @@ -703,59 +710,136 @@
}
}

void processDataLight(aod::ResoCollision const& collision,
aod::ResoTracks const& resotracks)
void processDataLight(aod::ResoCollision const& resocollisions, aod::ResoCollisionColls const& collisionIndex, soa::Join<aod::Collisions, aod::PVMults> const& collisions, aod::ResoTracks const& resotracks)
{
// LOG(info) << "new collision, zvtx: " << collision.posZ();
if (cisInelGt0) {
auto linkRow = collisionIndex.iteratorAt(resocollisions.globalIndex());
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision

auto coll = collisions.iteratorAt(collId); // Take original collision matched with resoCollision

if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision
return;
}
if (additionalQAeventPlots)
histos.fill(HIST("QAevent/hEvtCounterSameE"), 1.0);
fillHistograms<false, false>(collision, resotracks, resotracks);
auto multiplicity = resocollisions.cent();
fillHistograms<false, false>(resocollisions, resotracks, resotracks, multiplicity);
}
PROCESS_SWITCH(K892analysis, processDataLight, "Process Event for data", false);

void processMCLight(ResoMCCols::iterator const& collision,
soa::Join<aod::ResoTracks, aod::ResoMCTracks> const& resotracks)
void processMCLight(ResoMCCols::iterator const& resoCollision,
aod::ResoCollisionColls const& collisionIndex,
soa::Join<aod::ResoCollisionCandidatesMC, aod::PVMults> const& collisionsMC,
soa::Join<aod::ResoTracks, aod::ResoMCTracks> const& resoTracks,
soa::Join<aod::McCollisions, aod::McCentFT0Ms> const&)
{
if (!collision.isInAfterAllCuts() || (std::abs(collision.posZ()) > cZvertCutMC)) // MC event selection, all cuts missing vtx cut
float multiplicity;
if (cMCCent && cisInelGt0) {
auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex());
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision

auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision

if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision
return;

auto mcColl = coll.mcCollision_as<soa::Join<aod::McCollisions, aod::McCentFT0Ms>>();
multiplicity = mcColl.centFT0M();
} else if (!cMCCent && cisInelGt0) {
auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex());
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision

auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision

if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision
return;

multiplicity = resoCollision.cent();
} else if (cMCCent && !cisInelGt0) {
auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex());
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision

auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision

auto mcColl = coll.mcCollision_as<soa::Join<aod::McCollisions, aod::McCentFT0Ms>>();
multiplicity = mcColl.centFT0M();
} else {
multiplicity = resoCollision.cent();
}
if (!resoCollision.isInAfterAllCuts() || (std::abs(resoCollision.posZ()) > cZvertCutMC)) // MC event selection, all cuts missing vtx cut
return;
fillHistograms<true, false>(collision, resotracks, resotracks);
fillHistograms<true, false>(resoCollision, resoTracks, resoTracks, multiplicity);
}
PROCESS_SWITCH(K892analysis, processMCLight, "Process Event for MC (Reconstructed)", false);

void processMCTrue(ResoMCCols::iterator const& collision, aod::ResoMCParents const& resoParents)
void processMCTrue(ResoMCCols::iterator const& resoCollision, aod::ResoCollisionColls const& collisionIndex, aod::ResoMCParents const& resoParents, aod::ResoCollisionCandidatesMC const& collisionsMC, soa::Join<aod::McCollisions, aod::McCentFT0Ms> const&)
{
auto multiplicity = collision.cent();
float multiplicity;
if (cMCCent && cisInelGt0) {
auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex());
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision

auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision

if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision
return;

auto mcColl = coll.mcCollision_as<soa::Join<aod::McCollisions, aod::McCentFT0Ms>>();
multiplicity = mcColl.centFT0M();
} else if (!cMCCent && cisInelGt0) {
auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex());
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision

auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision

if (!coll.isInelGt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision
return;

multiplicity = resoCollision.cent();
} else if (cMCCent && !cisInelGt0) {
auto linkRow = collisionIndex.iteratorAt(resoCollision.globalIndex());
auto collId = linkRow.collisionId(); // Take original collision global index matched with resoCollision

auto coll = collisionsMC.iteratorAt(collId); // Take original collision matched with resoCollision

auto mcColl = coll.mcCollision_as<soa::Join<aod::McCollisions, aod::McCentFT0Ms>>();
multiplicity = mcColl.centFT0M();
} else {
multiplicity = resoCollision.cent();
}
for (const auto& part : resoParents) { // loop over all pre-filtered MC particles
if (std::abs(part.pdgCode()) != 313 || std::abs(part.y()) >= 0.5)

Check failure on line 813 in PWGLF/Tasks/Resonances/k892analysis.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;
bool pass1 = std::abs(part.daughterPDG1()) == 211 || std::abs(part.daughterPDG2()) == 211;

Check failure on line 815 in PWGLF/Tasks/Resonances/k892analysis.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.
bool pass2 = std::abs(part.daughterPDG1()) == 321 || std::abs(part.daughterPDG2()) == 321;

Check failure on line 816 in PWGLF/Tasks/Resonances/k892analysis.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.

if (!pass1 || !pass2)
continue;

if (collision.isVtxIn10()) // INEL10
if (resoCollision.isVtxIn10()) // INEL10
{
if (part.pdgCode() > 0)
histos.fill(HIST("k892Gen"), 0, part.pt(), multiplicity);
else
histos.fill(HIST("k892GenAnti"), 0, part.pt(), multiplicity);
}
if (collision.isVtxIn10() && collision.isInSel8()) // INEL>10, vtx10
if (resoCollision.isVtxIn10() && resoCollision.isInSel8()) // INEL>10, vtx10
{
if (part.pdgCode() > 0)
histos.fill(HIST("k892Gen"), 1, part.pt(), multiplicity);
else
histos.fill(HIST("k892GenAnti"), 1, part.pt(), multiplicity);
}
if (collision.isVtxIn10() && collision.isTriggerTVX()) // vtx10, TriggerTVX
if (resoCollision.isVtxIn10() && resoCollision.isTriggerTVX()) // vtx10, TriggerTVX
{
if (part.pdgCode() > 0)
histos.fill(HIST("k892Gen"), 2, part.pt(), multiplicity);
else
histos.fill(HIST("k892GenAnti"), 2, part.pt(), multiplicity);
}
if (collision.isInAfterAllCuts()) // after all event selection
if (resoCollision.isInAfterAllCuts()) // after all event selection
{
if (part.pdgCode() > 0)
histos.fill(HIST("k892Gen"), 3, part.pt(), multiplicity);
Expand All @@ -775,9 +859,10 @@
SameKindPair<aod::ResoCollisions, aod::ResoTracks, BinningTypeVtxZT0M> pairs{colBinning, nEvtMixing, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip

for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) {
auto multiplicity = collision1.cent();
if (additionalQAeventPlots)
histos.fill(HIST("QAevent/hEvtCounterMixedE"), 1.0);
fillHistograms<false, true>(collision1, tracks1, tracks2);
fillHistograms<false, true>(collision1, tracks1, tracks2, multiplicity);
}
};
PROCESS_SWITCH(K892analysis, processMELight, "Process EventMixing light without partition", false);
Expand Down
56 changes: 11 additions & 45 deletions PWGLF/Tasks/Resonances/kstarqa.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1688,7 +1688,7 @@ struct Kstarqa {
}
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&)
void processRec(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&, EventMCGenerated const&)
{

if (!collision.has_mcCollision()) {
Expand All @@ -1705,8 +1705,6 @@ struct Kstarqa {
}
// multiplicity = collision.centFT0M();

multiplicity = -1.0;

if (cSelectMultEstimator == kFT0M) {
multiplicity = collision.centFT0M();
} else if (cSelectMultEstimator == kFT0A) {
Expand Down Expand Up @@ -1934,20 +1932,22 @@ struct Kstarqa {
}
PROCESS_SWITCH(Kstarqa, processRec, "Process Reconstructed", false);

void processRec2(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&, aod::McCollisions const& /*mcCollisions*/)
void processRec2(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&, EventMCGenerated const&)
{

if (!collision.has_mcCollision()) {
return;
}

double multiplicityRec = -1.0;
const auto& mcCollisionRec = collision.mcCollision_as<EventMCGenerated>();
multiplicityRec = mcCollisionRec.centFT0M();

if (selectionConfig.isINELgt0 && !collision.isInelGt0()) {
return;
}
// multiplicity = collision.centFT0M();

multiplicity = -1.0;

if (cSelectMultEstimator == kFT0M) {
multiplicity = collision.centFT0M();
} else if (cSelectMultEstimator == kFT0A) {
Expand All @@ -1961,50 +1961,14 @@ struct Kstarqa {
}

hInvMass.fill(HIST("hAllRecCollisions"), multiplicity);
hInvMass.fill(HIST("hAllRecCollisionsCalib"), multiplicityRec);

if (!selectionEvent(collision, false)) {
return;
}

// // if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex || !collision.sel8()) {
// if (std::abs(collision.mcCollision().posZ()) > selectionConfig.cutzvertex) {
// return;
// }

// if (selectionConfig.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
// return;
// }

// if (selectionConfig.isTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) {
// 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();

multiplicity = -1.0;

if (cSelectMultEstimator == kFT0M) {
multiplicity = collision.centFT0M();
} else if (cSelectMultEstimator == kFT0A) {
multiplicity = collision.centFT0A();
} else if (cSelectMultEstimator == kFT0C) {
multiplicity = collision.centFT0C();
} else if (cSelectMultEstimator == kFV0A) {
multiplicity = collision.centFV0A();
} else {
multiplicity = collision.centFT0M(); // default
}
hInvMass.fill(HIST("h1RecMult"), multiplicity);
hInvMass.fill(HIST("h1RecMult2"), multiplicityRec);

hInvMass.fill(HIST("h1RecMult"), multiplicity);

Expand Down Expand Up @@ -2169,13 +2133,15 @@ struct Kstarqa {
mother = daughter1 + daughter2; // Kstar meson

hInvMass.fill(HIST("h2KstarRecpt2"), mothertrack1.pt(), multiplicity, std::sqrt(mothertrack1.e() * mothertrack1.e() - mothertrack1.p() * mothertrack1.p()));
hInvMass.fill(HIST("h2KstarRecptCalib2"), mothertrack1.pt(), multiplicityRec, std::sqrt(mothertrack1.e() * mothertrack1.e() - mothertrack1.p() * mothertrack1.p()));

if (applyRecMotherRapidity && mother.Rapidity() >= selectionConfig.rapidityMotherData) {
continue;
}

hInvMass.fill(HIST("h1KstarRecMass"), mother.M());
hInvMass.fill(HIST("h2KstarRecpt1"), mother.Pt(), multiplicity, mother.M());
hInvMass.fill(HIST("h2KstarRecptCalib1"), mother.Pt(), multiplicityRec, mother.M());
}
}
}
Expand Down
Loading