From 81f3728f0cf697bd2a9277c761b8128406f18d22 Mon Sep 17 00:00:00 2001 From: Swati Saha Date: Mon, 2 Mar 2026 16:13:27 +0530 Subject: [PATCH 1/3] Calculate observables as a func. of Ngen in MCgen and rec --- .../Tasks/MeanptFluctuations.cxx | 86 +++++++++++++++++-- 1 file changed, 80 insertions(+), 6 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx b/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx index 21ed21ad9de..d0548638c66 100644 --- a/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx @@ -27,6 +27,7 @@ #include "Framework/runDataProcessing.h" #include +#include "TDatabasePDG.h" #include #include #include @@ -149,6 +150,8 @@ struct MeanptFluctuationsAnalysis { HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; std::vector>> subSampleMcGen; std::vector>> subSample; + std::vector>> subSampleMcGenV2; + std::vector>> subSampleV2; TRandom3* fRndm = new TRandom3(0); // filtering collisions and tracks*********** @@ -246,6 +249,7 @@ struct MeanptFluctuationsAnalysis { histos.add("MCGenerated/hPtParticleVsTrack", "", kTH2F, {ptAxis, ptAxis}); histos.add("MCGenerated/hEtaParticleVsTrack", "", kTH2F, {{100, -2.01, 2.01}, {100, -2.01, 2.01}}); histos.add("MCGenerated/hPhiParticleVsTrack", "", kTH2F, {{100, 0., o2::constants::math::TwoPI}, {100, 0., o2::constants::math::TwoPI}}); + histos.add("MCGenerated/hNgenVsNrec", "Gen multiplicity vs. Rec multiplicity in MC", kTH2F, {multAxis, multAxis}); // Analysis Profiles for central val histos.add("MCGenerated/AnalysisProfiles/Prof_mean_t1", "", {HistType::kTProfile2D, {centAxis, multAxis}}); @@ -255,16 +259,39 @@ struct MeanptFluctuationsAnalysis { histos.add("MCGenerated/AnalysisProfiles/Hist2D_Nch_centrality", "", {HistType::kTH2D, {centAxis, multAxis}}); histos.add("MCGenerated/AnalysisProfiles/Hist2D_meanpt_centrality", "", {HistType::kTH2D, {centAxis, meanpTAxis}}); + histos.add("MCGenerated/AnalysisProfilesV2/Prof_mean_t1", "", {HistType::kTProfile, {multAxis}}); + histos.add("MCGenerated/AnalysisProfilesV2/Prof_var_t1", "", {HistType::kTProfile, {multAxis}}); + histos.add("MCGenerated/AnalysisProfilesV2/Prof_skew_t1", "", {HistType::kTProfile, {multAxis}}); + histos.add("MCGenerated/AnalysisProfilesV2/Prof_kurt_t1", "", {HistType::kTProfile, {multAxis}}); + histos.add("AnalysisProfilesV2/Prof_mean_t1", "", {HistType::kTProfile2D, {multAxis, multAxis}}); + histos.add("AnalysisProfilesV2/Prof_var_t1", "", {HistType::kTProfile2D, {multAxis, multAxis}}); + histos.add("AnalysisProfilesV2/Prof_skew_t1", "", {HistType::kTProfile2D, {multAxis, multAxis}}); + histos.add("AnalysisProfilesV2/Prof_kurt_t1", "", {HistType::kTProfile2D, {multAxis, multAxis}}); + // Analysis Profiles for error subSampleMcGen.resize(cfgNsubSample); + subSampleMcGenV2.resize(cfgNsubSample); + subSampleV2.resize(cfgNsubSample); for (int i = 0; i < cfgNsubSample; i++) { subSampleMcGen[i].resize(4); + subSampleMcGenV2[i].resize(4); + subSampleV2[i].resize(4); } for (int i = 0; i < cfgNsubSample; i++) { subSampleMcGen[i][0] = std::get>(histos.add(Form("MCGenerated/AnalysisProfiles/subSample_%d/Prof_mean_t1", i), "", {HistType::kTProfile2D, {centAxis, multAxis}})); subSampleMcGen[i][1] = std::get>(histos.add(Form("MCGenerated/AnalysisProfiles/subSample_%d/Prof_var_t1", i), "", {HistType::kTProfile2D, {centAxis, multAxis}})); subSampleMcGen[i][2] = std::get>(histos.add(Form("MCGenerated/AnalysisProfiles/subSample_%d/Prof_skew_t1", i), "", {HistType::kTProfile2D, {centAxis, multAxis}})); subSampleMcGen[i][3] = std::get>(histos.add(Form("MCGenerated/AnalysisProfiles/subSample_%d/Prof_kurt_t1", i), "", {HistType::kTProfile2D, {centAxis, multAxis}})); + + subSampleMcGenV2[i][0] = std::get>(histos.add(Form("MCGenerated/AnalysisProfilesV2/subSample_%d/Prof_mean_t1", i), "", {HistType::kTProfile, {multAxis}})); + subSampleMcGenV2[i][1] = std::get>(histos.add(Form("MCGenerated/AnalysisProfilesV2/subSample_%d/Prof_var_t1", i), "", {HistType::kTProfile, {multAxis}})); + subSampleMcGenV2[i][2] = std::get>(histos.add(Form("MCGenerated/AnalysisProfilesV2/subSample_%d/Prof_skew_t1", i), "", {HistType::kTProfile, {multAxis}})); + subSampleMcGenV2[i][3] = std::get>(histos.add(Form("MCGenerated/AnalysisProfilesV2/subSample_%d/Prof_kurt_t1", i), "", {HistType::kTProfile, {multAxis}})); + + subSampleV2[i][0] = std::get>(histos.add(Form("AnalysisProfilesV2/subSample_%d/Prof_mean_t1", i), "", {HistType::kTProfile2D, {multAxis, multAxis}})); + subSampleV2[i][1] = std::get>(histos.add(Form("AnalysisProfilesV2/subSample_%d/Prof_var_t1", i), "", {HistType::kTProfile2D, {multAxis, multAxis}})); + subSampleV2[i][2] = std::get>(histos.add(Form("AnalysisProfilesV2/subSample_%d/Prof_skew_t1", i), "", {HistType::kTProfile2D, {multAxis, multAxis}})); + subSampleV2[i][3] = std::get>(histos.add(Form("AnalysisProfilesV2/subSample_%d/Prof_kurt_t1", i), "", {HistType::kTProfile2D, {multAxis, multAxis}})); } } @@ -530,9 +557,11 @@ struct MeanptFluctuationsAnalysis { if (!mcParticle.has_mcCollision()) continue; - int pdgCode = std::abs(mcParticle.pdgCode()); - bool extraPDGType = (pdgCode != PDG_t::kK0Short && pdgCode != PDG_t::kLambda0); - if (extraPDGType && pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton) + // charged check + auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcParticle.pdgCode()); + if (!pdgEntry) + continue; + if (pdgEntry->Charge() == 0) continue; if (mcParticle.isPhysicalPrimary()) { @@ -575,6 +604,11 @@ struct MeanptFluctuationsAnalysis { histos.fill(HIST("MCGenerated/AnalysisProfiles/Hist2D_Nch_centrality"), cent, nChgen); histos.fill(HIST("MCGenerated/AnalysisProfiles/Hist2D_meanpt_centrality"), cent, meanTerm1gen); + histos.get(HIST("MCGenerated/AnalysisProfilesV2/Prof_mean_t1"))->Fill(nChgen, meanTerm1gen); + histos.get(HIST("MCGenerated/AnalysisProfilesV2/Prof_var_t1"))->Fill(nChgen, varianceTerm1gen); + histos.get(HIST("MCGenerated/AnalysisProfilesV2/Prof_skew_t1"))->Fill(nChgen, skewnessTerm1gen); + histos.get(HIST("MCGenerated/AnalysisProfilesV2/Prof_kurt_t1"))->Fill(nChgen, kurtosisTerm1gen); + // selecting subsample and filling profiles float lRandomMc = fRndm->Rndm(); int sampleIndexMc = static_cast(cfgNsubSample * lRandomMc); @@ -582,12 +616,17 @@ struct MeanptFluctuationsAnalysis { subSampleMcGen[sampleIndexMc][1]->Fill(cent, nChgen, varianceTerm1gen); subSampleMcGen[sampleIndexMc][2]->Fill(cent, nChgen, skewnessTerm1gen); subSampleMcGen[sampleIndexMc][3]->Fill(cent, nChgen, kurtosisTerm1gen); + + subSampleMcGenV2[sampleIndexMc][0]->Fill(nChgen, meanTerm1gen); + subSampleMcGenV2[sampleIndexMc][1]->Fill(nChgen, varianceTerm1gen); + subSampleMcGenV2[sampleIndexMc][2]->Fill(nChgen, skewnessTerm1gen); + subSampleMcGenV2[sampleIndexMc][3]->Fill(nChgen, kurtosisTerm1gen); } //------------------------------------------------------------------------------------------- } PROCESS_SWITCH(MeanptFluctuationsAnalysis, processMCGen, "Process Generated MC data", true); - void processMCRec(MyMCRecCollisions::iterator const& collision, MyMCTracks const& tracks, aod::McCollisions const&, aod::McParticles const&) + void processMCRec(MyMCRecCollisions::iterator const& collision, MyMCTracks const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles) { histos.fill(HIST("MCGenerated/hMC"), 5.5); @@ -635,6 +674,29 @@ struct MeanptFluctuationsAnalysis { histos.fill(HIST("Hist2D_globalTracks_PVTracks"), collision.multNTracksPV(), tracks.size()); histos.fill(HIST("Hist2D_cent_nch"), tracks.size(), centralityFT0C); + // Calculating generated no of particles for the collision event + double noGen = 0.0; + auto mcColl = collision.mcCollision(); + for (const auto& mcParticle : mcParticles) { + if (!mcParticle.has_mcCollision()) + continue; + + // charged check + auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcParticle.pdgCode()); + if (!pdgEntry) + continue; + if (pdgEntry->Charge() == 0) + continue; + + if (mcParticle.isPhysicalPrimary()) { + if ((mcParticle.pt() > cfgCutPtLower) && (mcParticle.pt() < cfgCutPreSelPt) && (std::abs(mcParticle.eta()) < cfgCutPreSelEta)) { + if (mcParticle.pt() > cfgCutPtLower && mcParticle.pt() < cfgCutPtUpper) { + noGen = noGen + 1.0; + } + } + } + } //! end particle loop + // variables double pTsum = 0.0; double nN = 0.0; @@ -667,7 +729,7 @@ struct MeanptFluctuationsAnalysis { } histos.fill(HIST("his2DdcaXYvsPtBeforePtDepSel"), track.pt(), track.dcaXY()); - histos.fill(HIST("his2DdcaZvsPtBeforePtDepSel"), track.pt(), track.dcaZ()); + histos.fill(HIST("his2DdcaZvsPtPtBeforePtDepSel"), track.pt(), track.dcaZ()); if (cfgUsePtDepDCAxy && !(std::abs(track.dcaXY()) < fPtDepDCAxy->Eval(track.pt()))) { continue; @@ -706,6 +768,8 @@ struct MeanptFluctuationsAnalysis { } } // end track loop + histos.fill(HIST("MCGenerated/hNgenVsNrec"), noGen, nCh); + // MeanPt if (nN > 0.0f) histos.fill(HIST("hMeanPt"), cent, pTsum / nN); @@ -725,6 +789,11 @@ struct MeanptFluctuationsAnalysis { histos.fill(HIST("AnalysisProfiles/Hist2D_Nch_centrality"), cent, nCh); histos.fill(HIST("AnalysisProfiles/Hist2D_meanpt_centrality"), cent, meanTerm1); + histos.get(HIST("AnalysisProfilesV2/Prof_mean_t1"))->Fill(noGen, nCh, meanTerm1); + histos.get(HIST("AnalysisProfilesV2/Prof_var_t1"))->Fill(noGen, nCh, varianceTerm1); + histos.get(HIST("AnalysisProfilesV2/Prof_skew_t1"))->Fill(noGen, nCh, skewnessTerm1); + histos.get(HIST("AnalysisProfilesV2/Prof_kurt_t1"))->Fill(noGen, nCh, kurtosisTerm1); + // selecting subsample and filling profiles float lRandom = fRndm->Rndm(); int sampleIndex = static_cast(cfgNsubSample * lRandom); @@ -732,6 +801,11 @@ struct MeanptFluctuationsAnalysis { subSample[sampleIndex][1]->Fill(cent, nCh, varianceTerm1); subSample[sampleIndex][2]->Fill(cent, nCh, skewnessTerm1); subSample[sampleIndex][3]->Fill(cent, nCh, kurtosisTerm1); + + subSampleV2[sampleIndex][0]->Fill(noGen, nCh, meanTerm1); + subSampleV2[sampleIndex][1]->Fill(noGen, nCh, varianceTerm1); + subSampleV2[sampleIndex][2]->Fill(noGen, nCh, skewnessTerm1); + subSampleV2[sampleIndex][3]->Fill(noGen, nCh, kurtosisTerm1); } //------------------------------------------------------------------------------------------- } @@ -804,7 +878,7 @@ struct MeanptFluctuationsAnalysis { } histos.fill(HIST("his2DdcaXYvsPtBeforePtDepSel"), track.pt(), track.dcaXY()); - histos.fill(HIST("his2DdcaZvsPtBeforePtDepSel"), track.pt(), track.dcaZ()); + histos.fill(HIST("his2DdcaZvsPtPtBeforePtDepSel"), track.pt(), track.dcaZ()); if (cfgUsePtDepDCAxy && !(std::abs(track.dcaXY()) < fPtDepDCAxy->Eval(track.pt()))) { continue; From 3fff5c47f2c56b444fd21fc6d9711bb848a8c793 Mon Sep 17 00:00:00 2001 From: Swati Saha Date: Mon, 2 Mar 2026 17:28:39 +0530 Subject: [PATCH 2/3] corrected the mc-particle loop --- PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx b/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx index d0548638c66..5a13c444592 100644 --- a/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx @@ -164,6 +164,7 @@ struct MeanptFluctuationsAnalysis { using EventCandidatesMC = soa::Join; Preslice perCollision = aod::track::collisionId; + Preslice perMcCollision = aod::mcparticle::mcCollisionId; // Event selection cuts - Alex TF1* fMultPVCutLow = nullptr; @@ -677,7 +678,10 @@ struct MeanptFluctuationsAnalysis { // Calculating generated no of particles for the collision event double noGen = 0.0; auto mcColl = collision.mcCollision(); - for (const auto& mcParticle : mcParticles) { + // Slice particles belonging only to this MC collision + auto particlesThisEvent = mcParticles.sliceBy(perMcCollision, mcColl.globalIndex()); + + for (const auto& mcParticle : particlesThisEvent) { if (!mcParticle.has_mcCollision()) continue; From 7ef9a68ece8747cd95085d35741e1c4c7b9a51f3 Mon Sep 17 00:00:00 2001 From: Swati Saha Date: Mon, 2 Mar 2026 17:37:51 +0530 Subject: [PATCH 3/3] corrected the mc-particle loop --- PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx b/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx index 5a13c444592..c67c9f767f8 100644 --- a/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/MeanptFluctuations.cxx @@ -733,7 +733,7 @@ struct MeanptFluctuationsAnalysis { } histos.fill(HIST("his2DdcaXYvsPtBeforePtDepSel"), track.pt(), track.dcaXY()); - histos.fill(HIST("his2DdcaZvsPtPtBeforePtDepSel"), track.pt(), track.dcaZ()); + histos.fill(HIST("his2DdcaZvsPtBeforePtDepSel"), track.pt(), track.dcaZ()); if (cfgUsePtDepDCAxy && !(std::abs(track.dcaXY()) < fPtDepDCAxy->Eval(track.pt()))) { continue; @@ -882,7 +882,7 @@ struct MeanptFluctuationsAnalysis { } histos.fill(HIST("his2DdcaXYvsPtBeforePtDepSel"), track.pt(), track.dcaXY()); - histos.fill(HIST("his2DdcaZvsPtPtBeforePtDepSel"), track.pt(), track.dcaZ()); + histos.fill(HIST("his2DdcaZvsPtBeforePtDepSel"), track.pt(), track.dcaZ()); if (cfgUsePtDepDCAxy && !(std::abs(track.dcaXY()) < fPtDepDCAxy->Eval(track.pt()))) { continue;