From bae1326e6d8e6f8a6cfecc93890bd33edd69271a Mon Sep 17 00:00:00 2001 From: gianniliveraro Date: Wed, 25 Mar 2026 13:06:48 -0300 Subject: [PATCH 1/4] First EMCal implementation on sigma0builder --- PWGLF/DataModel/LFSigmaTables.h | 36 ++ .../TableProducer/Strangeness/CMakeLists.txt | 2 +- .../Strangeness/sigma0builder.cxx | 562 +++++++++++++++--- 3 files changed, 510 insertions(+), 90 deletions(-) diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index 9d728bab596..0008c8bbf04 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -502,6 +502,42 @@ DECLARE_SOA_TABLE(Sigma0PhotonExtras, "AOD", "SIGMA0PHOTON", sigma0PhotonExtra::PhotonNegTrackCode, sigma0PhotonExtra::PhotonV0Type); + +// For EMCal Photon extra info +namespace sigma0EMPhoton +{ +//______________________________________________________ +// REGULAR COLUMNS FOR SIGMA0EMPHOTON +DECLARE_SOA_COLUMN(PhotonID, photonID, int); +DECLARE_SOA_COLUMN(PhotonEnergy, photonEnergy, float); +DECLARE_SOA_COLUMN(PhotonEta, photonEta, float); +DECLARE_SOA_COLUMN(PhotonPhi, photonPhi, float); +DECLARE_SOA_COLUMN(PhotonM02, photonM02, float); +DECLARE_SOA_COLUMN(PhotonM20, photonM20, float); +DECLARE_SOA_COLUMN(PhotonNCells, photonNCells, int); +DECLARE_SOA_COLUMN(PhotonTime, photonTime, float); +DECLARE_SOA_COLUMN(PhotonIsExotic, photonIsExotic, bool); +DECLARE_SOA_COLUMN(PhotonDistToBad, photonDistToBad, float); +DECLARE_SOA_COLUMN(PhotonNLM, photonNLM, int); +DECLARE_SOA_COLUMN(PhotonDefinition, photonDefinition, int); + +} // namespace sigma0PhotonExtra + +DECLARE_SOA_TABLE(sigma0EMPhotons, "AOD", "SIGMA0EMPHOTON", + sigma0EMPhoton::PhotonID, + sigma0EMPhoton::PhotonEnergy, + sigma0EMPhoton::PhotonEta, + sigma0EMPhoton::PhotonPhi, + sigma0EMPhoton::PhotonM02, + sigma0EMPhoton::PhotonM20, + sigma0EMPhoton::PhotonNCells, + sigma0EMPhoton::PhotonTime, + sigma0EMPhoton::PhotonIsExotic, + sigma0EMPhoton::PhotonDistToBad, + sigma0EMPhoton::PhotonNLM, + sigma0EMPhoton::PhotonDefinition); + + // For Lambda extra info namespace sigma0LambdaExtra { diff --git a/PWGLF/TableProducer/Strangeness/CMakeLists.txt b/PWGLF/TableProducer/Strangeness/CMakeLists.txt index 3f1d8565951..08ec72c6913 100644 --- a/PWGLF/TableProducer/Strangeness/CMakeLists.txt +++ b/PWGLF/TableProducer/Strangeness/CMakeLists.txt @@ -139,7 +139,7 @@ o2physics_add_dpl_workflow(cascademlselection o2physics_add_dpl_workflow(sigma0builder SOURCES sigma0builder.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore O2Physics::AnalysisCCDB + PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2Physics::AnalysisCore O2Physics::MLCore O2Physics::AnalysisCCDB COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(lambdajetpolarizationbuilder diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index cc6dcf403cf..ff11db32d73 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -24,6 +24,7 @@ #include "PWGLF/DataModel/LFStrangenessMLTables.h" #include "PWGLF/DataModel/LFStrangenessPIDTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGJE/DataModel/EMCALClusters.h" #include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" @@ -66,8 +67,15 @@ using V0DerivedMCDatas = soa::Join; using V0TOFDerivedMCDatas = soa::Join; +using EMCalMCClusters = soa::Join; + +// TODO: Should I filter the clusters? +// we need to filter the emcal clusters, since we only want clusters from one specific clusterizer algorithm +//using EMCalClusters = o2::soa::Filtered; + static const std::vector DirList = {"V0BeforeSel", "PhotonSel", "LambdaSel", "KShortSel"}; -; +static const std::vector DirList2 = {"EMCalPhotonBeforeSel", "EMCalPhotonSel"}; + struct sigma0builder { Service ccdb; @@ -86,15 +94,16 @@ struct sigma0builder { //__________________________________________________ // Sigma0 specific - Produces sigmaIndices; // references V0Cores from sigma0Gens - Produces sigma0cores; // sigma0 candidates info for analysis - Produces sigmaPhotonExtras; // photons from sigma0 candidates info - Produces sigmaLambdaExtras; // lambdas from sigma0 candidates info - Produces sigma0CollRefs; // references collisions from Sigma0Cores - Produces sigma0mccores; // Reco sigma0 MC properties - Produces sigma0mclabel; // Link of reco sigma0 to mcparticles - Produces sigma0Gens; // Generated sigma0s - Produces sigma0GenCollRefs; // references collisions from sigma0Gens + Produces sigmaIndices; // references V0Cores from sigma0Gens + Produces sigma0cores; // sigma0 candidates info for analysis + Produces sigmaPhotonExtras; // photons from sigma0 candidates info + Produces sigmaEmCalPhotonExtras; // EMCAL photons from sigma0 candidates info + Produces sigmaLambdaExtras; // lambdas from sigma0 candidates info + Produces sigma0CollRefs; // references collisions from Sigma0Cores + Produces sigma0mccores; // Reco sigma0 MC properties + Produces sigma0mclabel; // Link of reco sigma0 to mcparticles + Produces sigma0Gens; // Generated sigma0s + Produces sigma0GenCollRefs; // references collisions from sigma0Gens //__________________________________________________ // Pi0 specific @@ -161,8 +170,13 @@ struct sigma0builder { Configurable minIR{"minIR", -1, "minimum IR collisions"}; Configurable maxIR{"maxIR", -1, "maximum IR collisions"}; + Configurable fSkipEmptyEMCal{"fSkipEmptyEMCal", true, "Flag to skip events without EMCal clusters"}; + } eventSelections; + // Photon Source + //Configurable fUsePCMPhotons{"fUsePCMPhotons", true, "Use PCM Photons for sigma0/kstar reconstruction. If False, EMCal photons are used instead."}; + // Tables to fill Configurable fillPi0Tables{"fillPi0Tables", false, "fill pi0 tables for QA"}; Configurable fillSigma0Tables{"fillSigma0Tables", true, "fill sigma0 tables for analysis"}; @@ -178,7 +192,7 @@ struct sigma0builder { std::string prefix = "lambdaSelections"; // JSON group name Configurable Lambda_MLThreshold{"Lambda_MLThreshold", 0.1, "Decision Threshold value to select lambdas"}; Configurable AntiLambda_MLThreshold{"AntiLambda_MLThreshold", 0.1, "Decision Threshold value to select antilambdas"}; - Configurable doMCAssociation{"doMCAssociation", false, "if MC, select true lambda/alambdas only"}; + Configurable doMCAssociation{"doMCAssociation", false, "if MC, select true lambda/alambdas only"}; Configurable LambdaMinDCANegToPv{"LambdaMinDCANegToPv", .05, "min DCA Neg To PV (cm)"}; Configurable LambdaMinDCAPosToPv{"LambdaMinDCAPosToPv", .05, "min DCA Pos To PV (cm)"}; Configurable LambdaMaxDCAV0Dau{"LambdaMaxDCAV0Dau", 2.5, "Max DCA V0 Daughters (cm)"}; @@ -207,7 +221,7 @@ struct sigma0builder { struct : ConfigurableGroup { std::string prefix = "photonSelections"; // JSON group name Configurable Gamma_MLThreshold{"Gamma_MLThreshold", 0.1, "Decision Threshold value to select gammas"}; - Configurable doMCAssociation{"doMCAssociation", false, "if MC, select true photons only"}; + Configurable doMCAssociation{"doMCAssociation", false, "if MC, select true photons only"}; Configurable Photonv0TypeSel{"Photonv0TypeSel", 7, "select on a certain V0 type (leave negative if no selection desired)"}; Configurable PhotonMinDCADauToPv{"PhotonMinDCADauToPv", 0.0, "Min DCA daughter To PV (cm)"}; Configurable PhotonMaxDCAV0Dau{"PhotonMaxDCAV0Dau", 3.5, "Max DCA V0 Daughters (cm)"}; @@ -232,11 +246,27 @@ struct sigma0builder { Configurable PhotonPhiMax2{"PhotonPhiMax2", -1, "Phi min value to reject photons, region 2 (leave negative if no selection desired)"}; } photonSelections; + //// Photon criteria: + struct : ConfigurableGroup { + std::string prefix = "EMCalPhotonSelections"; // JSON group name + Configurable definition{"definition", 13, "Cluster definitions to be accepted (e.g. 13 for kV3MostSplitLowSeed)"}; + Configurable MinCells{"MinCells", 1, "Min number of cells in cluster"}; + Configurable MinEnergy{"MinEnergy", 0.0, "Minimum energy of selected clusters (GeV)"}; + Configurable MaxEnergy{"MaxEnergy", 5, "Max energy of selected clusters (GeV)"}; + Configurable MaxEta{"MaxEta", 1.0, "Max absolute cluster Eta"}; + Configurable MinTime{"MinTime", -50, "Minimum time of selected clusters (ns)"}; + Configurable MaxTime{"MaxTime", 50, "Max time of selected clusters (ns)"}; + Configurable RemoveExotic{"RemoveExotic", false, "Flag to enable the removal of exotic clusters"}; + Configurable MinM02{"MinM02", -1., "Minimum shower shape long axis"}; + Configurable MaxM02{"MaxM02", 5., "Max shower shape long axis"}; + + } EMCalPhotonSelections; + // KShort criteria: struct : ConfigurableGroup { std::string prefix = "kshortSelections"; // JSON group name Configurable KShort_MLThreshold{"KShort_MLThreshold", 0.1, "Decision Threshold value to select kshorts"}; - Configurable doMCAssociation{"doMCAssociation", false, "if MC, select true kshorts only"}; + Configurable doMCAssociation{"doMCAssociation", false, "if MC, select true kshorts only"}; Configurable KShortMinDCANegToPv{"KShortMinDCANegToPv", .05, "min DCA Neg To PV (cm)"}; Configurable KShortMinDCAPosToPv{"KShortMinDCAPosToPv", .05, "min DCA Pos To PV (cm)"}; Configurable KShortMaxDCAV0Dau{"KShortMaxDCAV0Dau", 2.5, "Max DCA V0 Daughters (cm)"}; @@ -314,8 +344,16 @@ struct sigma0builder { ConfigurableAxis axisCandSel{"axisCandSel", {15, 0.5f, +15.5f}, "Candidate Selection"}; ConfigurableAxis axisIRBinning{"axisIRBinning", {151, -10, 1500}, "Binning for the interaction rate (kHz)"}; ConfigurableAxis axisLifetime{"axisLifetime", {200, 0, 50}, "Lifetime"}; + + // EMCal-specifc + ConfigurableAxis axisClrDefinition{"axisClrDefinition", {51, -0.5, 50.5}, "Cluster Definition"}; + ConfigurableAxis axisClrNCells{"axisClrNCells", {25, 0.0, 25}, "N cells per cluster"}; + ConfigurableAxis axisClrEnergy{"axisClrEnergy", {400, 0.0, 10}, "Energy per cluster"}; + ConfigurableAxis axisClrTime{"axisClrTime", {300, -30.0, 30.0}, "cluster time (ns)"}; + ConfigurableAxis axisClrShape{"axisClrShape", {100, 0.0, 1.0}, "cluster shape"}; // For manual sliceBy (necessary to calculate the correction factors) + // TODO: remove me! PresliceUnsorted> perMcCollision = aod::v0data::straMCCollisionId; void init(InitContext const&) @@ -323,8 +361,10 @@ struct sigma0builder { LOGF(info, "Initializing now: cross-checking correctness..."); if (doprocessRealData + doprocessRealDataWithTOF + + doprocessRealDataWithEMCal + doprocessMonteCarlo + doprocessMonteCarloWithTOF + + doprocessMonteCarloWithEMCal + doprocessV0QA + doprocessV0MCQA > 1) { @@ -373,9 +413,11 @@ struct sigma0builder { } } + bool fUsePCMPhoton = !doprocessRealDataWithEMCal && !doprocessMonteCarloWithEMCal; + for (const auto& histodir : DirList) { if ((histodir == "V0BeforeSel" && !fFillNoSelV0Histos) || - (histodir == "PhotonSel" && !fFillSelPhotonHistos) || + (histodir == "PhotonSel" && !fFillSelPhotonHistos && !fUsePCMPhoton) || (histodir == "LambdaSel" && !fFillSelLambdaHistos) || (histodir == "KShortSel" && !fFillSelKShortHistos)) { continue; @@ -424,21 +466,46 @@ struct sigma0builder { histos.add(histodir + "/h3dV0XYZ", "h3dV0XYZ", kTH3D, {axisXY, axisXY, axisZ}); } - histos.add("PhotonSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "Mass"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(3, "Y"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(4, "Neg Eta"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(5, "Pos Eta"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(6, "DCAToPV"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(7, "DCADau"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(8, "Radius"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(9, "Z"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(10, "CosPA"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(11, "Phi"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(12, "TPCCR"); - histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(13, "TPC NSigma"); + if (fUsePCMPhoton){ + histos.add("PhotonSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "Mass"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(3, "Y"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(4, "Neg Eta"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(5, "Pos Eta"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(6, "DCAToPV"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(7, "DCADau"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(8, "Radius"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(9, "Z"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(10, "CosPA"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(11, "Phi"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(12, "TPCCR"); + histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(13, "TPC NSigma"); + } + else{ + for (const auto& histodir : DirList2) { + histos.add(histodir + "/hpT", "hpT", kTH1D, {axisPt}); + histos.add(histodir + "/hDefinition", "hDefinition", kTH1D, {axisClrDefinition}); + histos.add(histodir + "/hNCells", "hNCells", kTH1D, {axisClrNCells}); + histos.add(histodir + "/hEnergy", "hEnergy", kTH1D, {axisClrEnergy}); + histos.add(histodir + "/h2dEtaVsPhi", "h2dEtaVsPhi", kTH2D, {axisRapidity, axisPhi}); + histos.add(histodir + "/hTime", "hTime", kTH1D, {axisClrTime}); + histos.add(histodir + "/hExotic", "hExotic", kTH1D, {{2, -0.5f, 1.5f}}); + histos.add(histodir + "/hShape", "hShape", kTH1D, {axisClrShape}); + } + + histos.add("EMCalPhotonSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); + histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); + histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "Definition"); + histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(3, "MinCell"); + histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(4, "Energy"); + histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(5, "Eta"); + histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(6, "Time"); + histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(7, "Exotic"); + histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(8, "Shape"); + } + histos.add("LambdaSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); histos.get(HIST("LambdaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); histos.get(HIST("LambdaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "Mass"); @@ -471,7 +538,7 @@ struct sigma0builder { histos.get(HIST("KShortSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(13, "ITSNCls"); histos.get(HIST("KShortSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(14, "Lifetime"); - if (doprocessRealData || doprocessRealDataWithTOF || doprocessMonteCarlo || doprocessMonteCarloWithTOF) { + if (doprocessRealData || doprocessRealDataWithTOF || doprocessRealDataWithEMCal || doprocessMonteCarlo || doprocessMonteCarloWithTOF || doprocessMonteCarloWithEMCal) { histos.add("SigmaSel/hSigma0DauDeltaIndex", "hSigma0DauDeltaIndex", kTH1F, {{100, -49.5f, 50.5f}}); histos.add("SigmaSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); histos.get(HIST("SigmaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); @@ -501,7 +568,7 @@ struct sigma0builder { } // MC - if (doprocessMonteCarlo || doprocessMonteCarloWithTOF) { + if (doprocessMonteCarlo || doprocessMonteCarloWithTOF || doprocessMonteCarloWithEMCal) { histos.add("MCQA/h2dPhotonNMothersVsPDG", "h2dPhotonNMothersVsPDG", kTHnSparseD, {{10, -0.5f, +9.5f}, {10001, -5000.5f, +5000.5f}}); histos.add("MCQA/h2dPhotonNMothersVsMCProcess", "h2dPhotonNMothersVsMCProcess", kTH2D, {{10, -0.5f, +9.5f}, {50, -0.5f, 49.5f}}); histos.add("MCQA/hPhotonMotherSize", "hPhotonMotherSize", kTH1D, {{10, -0.5f, +9.5f}}); @@ -933,6 +1000,90 @@ struct sigma0builder { return MCinfo; } + template + V0PairMCInfo getClusterV0PairMCInfo(TEMCalCls const& cluster, TV0 const& v0, TCollision const& collision, TMCParticles const& mcparticles) + { + // TODO: you need to improve/recheck this + V0PairMCInfo MCinfo; + + // _________________________ + // V0-specific treatment: + if (!v0.has_v0MCCore()) { + return MCinfo; + } + + auto v0MC = v0.template v0MCCore_as>(); + auto MCParticle_v0 = mcparticles.rawIteratorAt(v0MC.particleIdMC()); + auto const& MCMothersList_v0 = MCParticle_v0.template mothers_as(); + + if (!MCMothersList_v0.empty()) + return MCinfo; + + if (collision.has_straMCCollision()) { + auto MCCollision = collision.template straMCCollision_as>(); + MCinfo.fIsV02CorrectlyAssign = (v0MC.straMCCollisionId() == MCCollision.globalIndex()); + } + + // Getting Lambda's mother + auto const& MCMother_v0 = MCMothersList_v0.front(); // First mother + + // _________________________ + // EMCalClusters-specific treatment: + + // We loop over MC contributors to the cluster + for (size_t i = 0; i < cluster.amplitudeA().size(); i++) { + + // Get MC particle + int clusterInducerId = cluster.mcParticleIds()[i]; + auto mcClusterInducer = mcparticles.iteratorAt(clusterInducerId); + + // Get MC Mother + auto const& MCMothersList_Cluster = mcClusterInducer.template mothers_as(); + + if (MCMothersList_Cluster.empty()) + continue; + + auto const& MCMother_cluster = MCMothersList_Cluster.front(); // First mother + + if (MCMother_cluster.globalIndex() == MCMother_v0.globalIndex()){ + // Is it the same mother? + MCinfo.fV0PairProducedByGenerator = MCMother_v0.producedByGenerator(); + MCinfo.V0PairPDGCode = MCMother_v0.pdgCode(); + MCinfo.V0PairMCProcess = MCMother_v0.getProcess(); + MCinfo.V0PairMCParticleID = MCMother_v0.globalIndex(); + MCinfo.V0PairMCRadius = std::hypot(MCMother_v0.vx(), MCMother_v0.vy()); // production position radius + + auto const& v0pairmothers = MCMother_v0.template mothers_as(); // Get mothers + if (!v0pairmothers.empty()) { + auto& v0PairMother = v0pairmothers.front(); // V0Pair mother, V0s grandmother + MCinfo.V0PairPDGCodeMother = v0PairMother.pdgCode(); + } + + + // Basic kinematic info + MCinfo.V01MCpx = mcClusterInducer.px(); + MCinfo.V01MCpy = mcClusterInducer.py(); + MCinfo.V01MCpz = mcClusterInducer.pz(); + + MCinfo.V02MCpx = v0MC.pxMC(); + MCinfo.V02MCpy = v0MC.pyMC(); + MCinfo.V02MCpz = v0MC.pzMC(); + + // MC association info + MCinfo.fIsV01Primary = mcClusterInducer.isPhysicalPrimary(); + MCinfo.fIsV02Primary = v0MC.isPhysicalPrimary(); + MCinfo.V01PDGCode = mcClusterInducer.pdgCode(); + MCinfo.V02PDGCode = v0MC.pdgCode(); + MCinfo.V01PDGCodeMother = MCMother_cluster.pdgCode(); + MCinfo.V02PDGCodeMother = v0MC.pdgCodeMother(); + + break; // we found what we were looking for! + } + } + + return MCinfo; + } + // ______________________________________________________ // Check whether the collision passes our collision selections // Should work with collisions, mccollisions, stracollisions and stramccollisions tables! @@ -1501,9 +1652,29 @@ struct sigma0builder { } } + // Function to fill QA histograms. mode = 0 (before selections, all clusters), 1 after all selections + template + void fillEMCalHistos(TEMCalClusterObject const& cluster) + { + // Check whether it is before or after selections + static constexpr std::string_view MainDir2[] = {"EMCalPhotonBeforeSel", "EMCalPhotonSel"}; + + // calculate pT for cluster assuming they are photons (so no mass) + float gammapT = sqrt(cluster.energy() * cluster.energy()) / std::cosh(cluster.eta()); + + histos.fill(HIST(MainDir2[mode]) + HIST("/hpT"), gammapT); + histos.fill(HIST(MainDir2[mode]) + HIST("/hDefinition"), cluster.definition()); + histos.fill(HIST(MainDir2[mode]) + HIST("/hNCells"), cluster.nCells()); + histos.fill(HIST(MainDir2[mode]) + HIST("/hEnergy"), cluster.energy()); + histos.fill(HIST(MainDir2[mode]) + HIST("/h2dEtaVsPhi"), cluster.eta(), cluster.phi()); + histos.fill(HIST(MainDir2[mode]) + HIST("/hTime"), cluster.time()); + histos.fill(HIST(MainDir2[mode]) + HIST("/hExotic"), cluster.isExotic()); + histos.fill(HIST(MainDir2[mode]) + HIST("/hShape"), cluster.m02()); + } + // Function to fill QA histograms. mode = 0 (before selections, all v0s), 1 (photon candidates), 2 (lambda/alambda candidates) template - void fillHistos(TV0Object const& v0, TCollision const& collision) + void fillV0Histos(TV0Object const& v0, TCollision const& collision) { // Check whether it is before or after selections static constexpr std::string_view MainDir[] = {"V0BeforeSel", "PhotonSel", "LambdaSel", "KShortSel"}; @@ -1567,7 +1738,56 @@ struct sigma0builder { } //_______________________________________________ - // Process photon candidate + // Process v0 photon candidate + template + bool processEMCalPhotonCandidate(TEMCalClusterObject const& cluster) + { + // Clusterizer + histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 1.); + if (cluster.definition() != EMCalPhotonSelections.definition && EMCalPhotonSelections.definition > -1) + return false; + + // Number of Cells + histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 2.); + if (cluster.nCells() < EMCalPhotonSelections.MinCells) + return false; + + // Energy + histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 3.); + if (cluster.energy() < EMCalPhotonSelections.MinEnergy || cluster.energy() > EMCalPhotonSelections.MaxEnergy) + return false; + + // Eta + histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 4.); + if (TMath::Abs(cluster.eta()) > EMCalPhotonSelections.MaxEta) + return false; + + // Timing + histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 5.); + if (cluster.time() < EMCalPhotonSelections.MinTime || cluster.time() > EMCalPhotonSelections.MaxTime) + return false; + + // Exotic Clusters + histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 6.); + if (cluster.isExotic() && EMCalPhotonSelections.RemoveExotic) { + return false; + } + + // Shower shape long axis + histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 7.); + if (cluster.nCells() > 1) { // Only if we have more than one + if (cluster.m02() < EMCalPhotonSelections.MinM02 || cluster.m02() > EMCalPhotonSelections.MaxM02) { + return false; + } + } + + histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 8.); + + return true; + } + + //_______________________________________________ + // Process v0 photon candidate template bool processPhotonCandidate(TV0Object const& gamma) { @@ -1918,14 +2138,16 @@ struct sigma0builder { } //_______________________________________________ - // Build sigma0 candidate + // Build sigma0 candidate with PCM photons template - bool buildSigma0(TV0Object const& lambda, TV0Object const& gamma, TCollision const& collision, TMCParticles const& mcparticles) + bool buildPCMSigma0(TV0Object const& lambda, TV0Object const& gamma, TCollision const& collision, TMCParticles const& mcparticles) { //_______________________________________________ + // TODO: Add Same index rejection + // Checking if both V0s are made of the very same tracks if (gamma.posTrackExtraId() == lambda.posTrackExtraId() || - gamma.negTrackExtraId() == lambda.negTrackExtraId()) { + gamma.negTrackExtraId() == lambda.negTrackExtraId()) { // TODO: Add extra conditions return false; } @@ -1966,19 +2188,20 @@ struct sigma0builder { if constexpr (requires { gamma.motherMCPartId(); lambda.motherMCPartId(); }) { auto sigma0MCInfo = getV0PairMCInfo(gamma, lambda, collision, mcparticles); + // TODO: Include PDG of V0 Daughters sigma0mccores(sigma0MCInfo.V0PairMCRadius, sigma0MCInfo.V0PairPDGCode, sigma0MCInfo.V0PairPDGCodeMother, sigma0MCInfo.V0PairMCProcess, sigma0MCInfo.fV0PairProducedByGenerator, sigma0MCInfo.V01MCpx, sigma0MCInfo.V01MCpy, sigma0MCInfo.V01MCpz, sigma0MCInfo.fIsV01Primary, sigma0MCInfo.V01PDGCode, sigma0MCInfo.V01PDGCodeMother, sigma0MCInfo.fIsV01CorrectlyAssign, sigma0MCInfo.V02MCpx, sigma0MCInfo.V02MCpy, sigma0MCInfo.V02MCpz, sigma0MCInfo.fIsV02Primary, sigma0MCInfo.V02PDGCode, sigma0MCInfo.V02PDGCodeMother, sigma0MCInfo.fIsV02CorrectlyAssign); - sigma0mclabel(sigma0MCInfo.V0PairMCParticleID); + sigma0mclabel(sigma0MCInfo.V0PairMCParticleID); // TODO: Add this in sigma0mccores } // Sigma0s -> stracollisions link histos.fill(HIST("SigmaSel/hSigma0DauDeltaIndex"), gamma.globalIndex() - lambda.globalIndex()); sigma0CollRefs(collision.globalIndex()); - sigmaIndices(gamma.globalIndex(), lambda.globalIndex()); + sigmaIndices(gamma.globalIndex(), lambda.globalIndex()); // TODO: put these indices in sigma0core table //_______________________________________________ // Photon extra properties @@ -2043,6 +2266,117 @@ struct sigma0builder { return true; } + // Build sigma0 candidate + template + bool buildEMCalSigma0(TV0Object const& lambda, TEMCalClsObject const& gamma, TCollision const& collision, TMCParticles const& mcparticles) + { + // calculate pT for cluster assuming they are photons (so no mass) + float gammapT = sqrt(gamma.energy() * gamma.energy()) / std::cosh(gamma.eta()); + + // Momentum components + float gammapx = gammapT * std::cos(gamma.phi()); + float gammapy = gammapT * std::sin(gamma.phi()); + float gammapz = gammapT * std::sinh(gamma.phi()); + + //_______________________________________________ + // Sigma0 pre-selections + std::array pVecPhotons{gammapx, gammapy, gammapz}; + std::array pVecLambda{lambda.px(), lambda.py(), lambda.pz()}; + + auto arrMom = std::array{pVecPhotons, pVecLambda}; + float sigmaMass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassLambda0}); + float sigmaY = -999.f; + + // TODO: do proper calculation of rapidity if MC + sigmaY = RecoDecay::y(std::array{gammapx + lambda.px(), gammapy + lambda.py(), gammapz + lambda.pz()}, o2::constants::physics::MassSigma0); + + // if constexpr (requires { gamma.pxMC(); lambda.pxMC(); }) // If MC + // sigmaY = RecoDecay::y(std::array{gamma.pxMC() + lambda.pxMC(), gamma.pyMC() + lambda.pyMC(), gamma.pzMC() + lambda.pzMC()}, o2::constants::physics::MassSigma0); + // else // If DATA + // sigmaY = RecoDecay::y(std::array{gamma.px() + lambda.px(), gamma.py() + lambda.py(), gamma.pz() + lambda.pz()}, o2::constants::physics::MassSigma0); + + histos.fill(HIST("SigmaSel/hSelectionStatistics"), 1.); + if (TMath::Abs(sigmaMass - o2::constants::physics::MassSigma0) > Sigma0Window) + return false; + + histos.fill(HIST("SigmaSel/hSelectionStatistics"), 2.); + if (TMath::Abs(sigmaY) > SigmaMaxRap) + return false; + + histos.fill(HIST("SigmaSel/hSigmaMassSelected"), sigmaMass); + histos.fill(HIST("SigmaSel/hSelectionStatistics"), 3.); + + //_______________________________________________ + // Calculate properties & Fill tables + sigma0cores(0.0f, 0.0f, 0.0f, 0.0f, // N.B: Filling with dummy values for now + gammapx, gammapy, gammapz, 0.0f, + lambda.px(), lambda.py(), lambda.pz(), lambda.mLambda(), lambda.mAntiLambda()); + + // MC properties + if constexpr (requires { gamma.mcParticleIds(); lambda.motherMCPartId(); }) { // TODO: Solve MC treatment + + auto sigma0MCInfo = getClusterV0PairMCInfo(gamma, lambda, collision, mcparticles); + + sigma0mccores(sigma0MCInfo.V0PairMCRadius, sigma0MCInfo.V0PairPDGCode, sigma0MCInfo.V0PairPDGCodeMother, sigma0MCInfo.V0PairMCProcess, sigma0MCInfo.fV0PairProducedByGenerator, + sigma0MCInfo.V01MCpx, sigma0MCInfo.V01MCpy, sigma0MCInfo.V01MCpz, + sigma0MCInfo.fIsV01Primary, sigma0MCInfo.V01PDGCode, sigma0MCInfo.V01PDGCodeMother, sigma0MCInfo.fIsV01CorrectlyAssign, + sigma0MCInfo.V02MCpx, sigma0MCInfo.V02MCpy, sigma0MCInfo.V02MCpz, + sigma0MCInfo.fIsV02Primary, sigma0MCInfo.V02PDGCode, sigma0MCInfo.V02PDGCodeMother, sigma0MCInfo.fIsV02CorrectlyAssign); + + sigma0mclabel(sigma0MCInfo.V0PairMCParticleID); + } + + sigma0CollRefs(collision.globalIndex()); + sigmaIndices(gamma.globalIndex(), lambda.globalIndex()); + + //_______________________________________________ + // Photon extra properties + // TODO: include EMCal photons table in header! + sigmaEmCalPhotonExtras(gamma.id(), gamma.energy(), gamma.eta(), gamma.phi(), + gamma.m02(), gamma.m20(), gamma.nCells(), gamma.time(), + gamma.isExotic(), gamma.distanceToBadChannel(), gamma.nlm(), gamma.definition()); + + //_______________________________________________ + // Lambda extra properties + auto posTrackLambda = lambda.template posTrackExtra_as(); + auto negTrackLambda = lambda.template negTrackExtra_as(); + + uint8_t fLambdaPosTrackCode = ((uint8_t(posTrackLambda.hasTPC()) << hasTPC) | + (uint8_t(posTrackLambda.hasITSTracker()) << hasITSTracker) | + (uint8_t(posTrackLambda.hasITSAfterburner()) << hasITSAfterburner) | + (uint8_t(posTrackLambda.hasTRD()) << hasTRD) | + (uint8_t(posTrackLambda.hasTOF()) << hasTOF)); + + uint8_t fLambdaNegTrackCode = ((uint8_t(negTrackLambda.hasTPC()) << hasTPC) | + (uint8_t(negTrackLambda.hasITSTracker()) << hasITSTracker) | + (uint8_t(negTrackLambda.hasITSAfterburner()) << hasITSAfterburner) | + (uint8_t(negTrackLambda.hasTRD()) << hasTRD) | + (uint8_t(negTrackLambda.hasTOF()) << hasTOF)); + + float fLambdaLifeTime = lambda.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassLambda0; + float fLambdaPrTOFNSigma = -999.f; + float fLambdaPiTOFNSigma = -999.f; + float fALambdaPrTOFNSigma = -999.f; + float fALambdaPiTOFNSigma = -999.f; + + if constexpr (requires { lambda.tofNSigmaLaPr(); }) { // If TOF info avaiable + fLambdaPrTOFNSigma = lambda.tofNSigmaLaPr(); + fLambdaPiTOFNSigma = lambda.tofNSigmaLaPi(); + fALambdaPrTOFNSigma = lambda.tofNSigmaALaPr(); + fALambdaPiTOFNSigma = lambda.tofNSigmaALaPi(); + } + + sigmaLambdaExtras(lambda.qtarm(), lambda.alpha(), fLambdaLifeTime, lambda.v0radius(), lambda.v0cosPA(), lambda.dcaV0daughters(), lambda.dcanegtopv(), lambda.dcapostopv(), + posTrackLambda.tpcNSigmaPr(), posTrackLambda.tpcNSigmaPi(), negTrackLambda.tpcNSigmaPr(), negTrackLambda.tpcNSigmaPi(), + fLambdaPrTOFNSigma, fLambdaPiTOFNSigma, fALambdaPrTOFNSigma, fALambdaPiTOFNSigma, + posTrackLambda.tpcCrossedRows(), negTrackLambda.tpcCrossedRows(), + lambda.positiveeta(), lambda.negativeeta(), + posTrackLambda.itsNCls(), negTrackLambda.itsNCls(), posTrackLambda.itsChi2PerNcl(), negTrackLambda.itsChi2PerNcl(), + fLambdaPosTrackCode, fLambdaNegTrackCode, lambda.v0Type()); + + return true; + } + //_______________________________________________ // Build kstar candidate template @@ -2165,25 +2499,32 @@ struct sigma0builder { } // Process photon and lambda candidates to build sigma0 candidates - template - void dataProcess(TCollision const& collisions, TV0s const& fullV0s, TMCParticles const& mcparticles) + template + void dataProcess(TCollision const& collisions, TV0s const& fullV0s, TEMCal const& fullEMCalClusters, TMCParticles const& mcparticles) { //_______________________________________________ // Initial setup + bool fUsePCMPhoton = !doprocessRealDataWithEMCal && !doprocessMonteCarloWithEMCal; + // Auxiliary vectors to store best candidates std::vector bestGammasArray; - std::vector bestLambdasArray; - - std::vector bestKStarGammasArray; + std::vector bestLambdasArray; std::vector bestKShortsArray; // Custom grouping - std::vector> v0grouped(collisions.size()); + std::vector> v0grouped(collisions.size()); + std::vector> emclustersgrouped(collisions.size()); for (const auto& v0 : fullV0s) { v0grouped[v0.straCollisionId()].push_back(v0.globalIndex()); } - + + if constexpr (soa::is_table) { + for (const auto& cluster : fullEMCalClusters) { + emclustersgrouped[cluster.collisionId()].push_back(cluster.globalIndex()); + } + } + //_______________________________________________ // Collisions loop for (const auto& coll : collisions) { @@ -2193,11 +2534,15 @@ struct sigma0builder { continue; } + // TODO: Check! Is this correct? + if constexpr (soa::is_table) { + if (emclustersgrouped[coll.globalIndex()].size()==0 && eventSelections.fSkipEmptyEMCal) + continue; + } + // Clear vectors bestGammasArray.clear(); - bestLambdasArray.clear(); - - bestKStarGammasArray.clear(); + bestLambdasArray.clear(); bestKShortsArray.clear(); float centrality = doPPAnalysis ? coll.centFT0M() : coll.centFT0C(); @@ -2209,41 +2554,55 @@ struct sigma0builder { auto v0 = fullV0s.rawIteratorAt(v0grouped[coll.globalIndex()][i]); if (fFillNoSelV0Histos) - fillHistos<0>(v0, coll); // Filling "all V0s" histograms - - if (processPhotonCandidate(v0)) { // selecting photons - if (fFillSelPhotonHistos) - fillHistos<1>(v0, coll); // QA histos - bestGammasArray.push_back(v0.globalIndex()); // Save indices of best gamma candidates + fillV0Histos<0>(v0, coll); // Filling "all V0s" histograms + + if (fUsePCMPhoton){ + if (processPhotonCandidate(v0)) { // selecting photons + if (fFillSelPhotonHistos) + fillV0Histos<1>(v0, coll); // QA histos + bestGammasArray.push_back(v0.globalIndex()); // Save indices of best gamma candidates + } } - + if (processLambdaCandidate(v0, coll)) { // selecting lambdas if (fFillSelLambdaHistos) - fillHistos<2>(v0, coll); // QA histos + fillV0Histos<2>(v0, coll); // QA histos bestLambdasArray.push_back(v0.globalIndex()); // Save indices of best lambda candidates } if (processKShortCandidate(v0, coll)) { // selecting kshorts if (fFillSelKShortHistos) - fillHistos<3>(v0, coll); // QA histos + fillV0Histos<3>(v0, coll); // QA histos bestKShortsArray.push_back(v0.globalIndex()); // Save indices of best kshort candidates } } - + + // If EMCalClusters is available, we don't use PCM + if constexpr (soa::is_table) { + for (size_t i = 0; i < emclustersgrouped[coll.globalIndex()].size(); i++) { + auto cluster = fullEMCalClusters.rawIteratorAt(emclustersgrouped[coll.globalIndex()][i]); + fillEMCalHistos<0>(cluster); + + if (processEMCalPhotonCandidate(cluster)) { // selecting photons + fillEMCalHistos<1>(cluster); // QA histos + bestGammasArray.push_back(cluster.globalIndex()); // Save indices of best gamma candidates + } + } + } + //_______________________________________________ // Wrongly collision association study (MC-specific) if constexpr (requires { coll.straMCCollisionId(); }) { if (doAssocStudy) { - analyzeV0CollAssoc(coll, fullV0s, bestGammasArray, true); // Photon-analysis + if constexpr (!soa::is_table) analyzeV0CollAssoc(coll, fullV0s, bestGammasArray, true); // Photon-analysis analyzeV0CollAssoc(coll, fullV0s, bestLambdasArray, false); // Lambda-analysis } } //_______________________________________________ - // V0 nested loop + // Photon-V0 nested loop for (size_t i = 0; i < bestGammasArray.size(); ++i) { - auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); - + //_______________________________________________ // Sigma0 loop if (fillSigma0Tables) { @@ -2251,34 +2610,47 @@ struct sigma0builder { auto lambda = fullV0s.rawIteratorAt(bestLambdasArray[j]); // Building sigma0 candidate & filling tables - if (!buildSigma0(lambda, gamma1, coll, mcparticles)) - continue; + if constexpr (soa::is_table) { // using EMCal photons + auto gamma1 = fullEMCalClusters.rawIteratorAt(bestGammasArray[i]); + if (!buildEMCalSigma0(lambda, gamma1, coll, mcparticles)) continue; + } + else { // using PCM photons + auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); + if (!buildPCMSigma0(lambda, gamma1, coll, mcparticles)) continue; + } } } //_______________________________________________ - // KStar loop - if (fillKStarTables) { - for (size_t j = 0; j < bestKShortsArray.size(); ++j) { - auto kshort = fullV0s.rawIteratorAt(bestKShortsArray[j]); - - // Building kstar candidate & filling tables - if (!buildKStar(kshort, gamma1, coll, mcparticles)) - continue; + // KStar loop + if constexpr (!soa::is_table){ // Don't use EMCal clusters here + if (fillKStarTables) { + auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); + for (size_t j = 0; j < bestKShortsArray.size(); ++j) { + auto kshort = fullV0s.rawIteratorAt(bestKShortsArray[j]); + + // Building kstar candidate & filling tables + if (!buildKStar(kshort, gamma1, coll, mcparticles)) + continue; + + } } } - + //_______________________________________________ // pi0 loop - if (fillPi0Tables) { - for (size_t j = i + 1; j < bestGammasArray.size(); ++j) { - auto gamma2 = fullV0s.rawIteratorAt(bestGammasArray[j]); - - // Building pi0 candidate & filling tables - if (!buildPi0(gamma1, gamma2, coll, mcparticles)) - continue; + if constexpr (!soa::is_table){ // Don't use EMCal clusters here + if (fillPi0Tables) { + auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); + for (size_t j = i + 1; j < bestGammasArray.size(); ++j) { + auto gamma2 = fullV0s.rawIteratorAt(bestGammasArray[j]); + + // Building pi0 candidate & filling tables + if (!buildPi0(gamma1, gamma2, coll, mcparticles)) + continue; + } } - } + } } } } @@ -2317,7 +2689,7 @@ struct sigma0builder { auto v0 = fullV0s.rawIteratorAt(v0grouped[coll.globalIndex()][i]); if (fFillNoSelV0Histos) - fillHistos<0>(v0, coll); // Filling "all V0s" histograms + fillV0Histos<0>(v0, coll); // Filling "all V0s" histograms // Selection process fPassPhotonSel = processPhotonCandidate(v0); @@ -2327,7 +2699,7 @@ struct sigma0builder { // Reco part: if (fPassPhotonSel) { if (fFillSelPhotonHistos) - fillHistos<1>(v0, coll); + fillV0Histos<1>(v0, coll); // Fill analysis Histos: float PhotonY = RecoDecay::y(std::array{v0.px(), v0.py(), v0.pz()}, o2::constants::physics::MassGamma); @@ -2337,7 +2709,7 @@ struct sigma0builder { } if (fPassLambdaSel) { if (fFillSelLambdaHistos) - fillHistos<2>(v0, coll); + fillV0Histos<2>(v0, coll); // Fill analysis Histos: histos.fill(HIST("V0QA/h3dLambdaMass"), centrality, v0.pt(), v0.mLambda()); @@ -2352,7 +2724,7 @@ struct sigma0builder { if (fPassKShortSel) { if (fFillSelKShortHistos) - fillHistos<3>(v0, coll); + fillV0Histos<3>(v0, coll); // Fill analysis Histos: histos.fill(HIST("V0QA/h3dKShortMass"), centrality, v0.pt(), v0.mK0Short()); @@ -2406,22 +2778,32 @@ struct sigma0builder { // Sigma0 processing part void processRealData(soa::Join const& collisions, V0StandardDerivedDatas const& fullV0s, dauTracks const&) { - dataProcess(collisions, fullV0s, nullptr); + dataProcess(collisions, fullV0s, nullptr, nullptr); } void processRealDataWithTOF(soa::Join const& collisions, V0TOFStandardDerivedDatas const& fullV0s, dauTracks const&) { - dataProcess(collisions, fullV0s, nullptr); + dataProcess(collisions, fullV0s, nullptr, nullptr); + } + + void processRealDataWithEMCal(soa::Join const& collisions, V0StandardDerivedDatas const& fullV0s, dauTracks const&, aod::EMCALClusters const& fullEMCalClusters) + { + dataProcess(collisions, fullV0s, fullEMCalClusters, nullptr); } void processMonteCarlo(soa::Join const& collisions, V0DerivedMCDatas const& fullV0s, aod::McParticles const& mcParticles, dauTracks const&, aod::MotherMCParts const&, soa::Join const&, soa::Join const&) { - dataProcess(collisions, fullV0s, mcParticles); + dataProcess(collisions, fullV0s, nullptr, mcParticles); } void processMonteCarloWithTOF(soa::Join const& collisions, V0TOFDerivedMCDatas const& fullV0s, aod::McParticles const& mcParticles, dauTracks const&, aod::MotherMCParts const&, soa::Join const&, soa::Join const&) { - dataProcess(collisions, fullV0s, mcParticles); + dataProcess(collisions, fullV0s, nullptr, mcParticles); + } + + void processMonteCarloWithEMCal(soa::Join const& collisions, V0DerivedMCDatas const& fullV0s, aod::McParticles const& mcParticles, dauTracks const&, aod::MotherMCParts const&, soa::Join const&, soa::Join const&, EMCalMCClusters const& fullEMCalMCClusters) + { + dataProcess(collisions, fullV0s, fullEMCalMCClusters, mcParticles); } void processGeneratedRun3(aod::McParticles const& mcParticles) @@ -2446,9 +2828,11 @@ struct sigma0builder { } PROCESS_SWITCH(sigma0builder, processRealData, "process as if real data", true); - PROCESS_SWITCH(sigma0builder, processRealDataWithTOF, "process as if real data", false); + PROCESS_SWITCH(sigma0builder, processRealDataWithTOF, "process as if real data, uses TOF PID info", false); + PROCESS_SWITCH(sigma0builder, processRealDataWithEMCal, "process as if real data, uses EMCal clusters", false); PROCESS_SWITCH(sigma0builder, processMonteCarlo, "process as if MC data", false); PROCESS_SWITCH(sigma0builder, processMonteCarloWithTOF, "process as if MC data, uses TOF PID info", false); + PROCESS_SWITCH(sigma0builder, processMonteCarloWithEMCal, "process as if MC data, uses EMCal clusters", false); PROCESS_SWITCH(sigma0builder, processGeneratedRun3, "process generated MC info", false); PROCESS_SWITCH(sigma0builder, processV0QA, "process QA of lambdas and photons", false); PROCESS_SWITCH(sigma0builder, processV0MCQA, "process QA of lambdas and photons", false); From b444b56478b22f217883e5aa3bfce261c15a22aa Mon Sep 17 00:00:00 2001 From: gianniliveraro Date: Mon, 30 Mar 2026 10:47:50 -0300 Subject: [PATCH 2/4] MC assoc fix + data model changes + analysis histos --- PWGLF/DataModel/LFSigmaTables.h | 50 +- .../Strangeness/sigma0builder.cxx | 673 +++++++++++++----- PWGLF/Tasks/Strangeness/sigmaanalysis.cxx | 397 +++++++---- 3 files changed, 800 insertions(+), 320 deletions(-) diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index 0008c8bbf04..00633578ad4 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -35,16 +35,11 @@ using std::array; namespace o2::aod { -// Indexing -namespace sigma0Core -{ -DECLARE_SOA_INDEX_COLUMN_FULL(PhotonV0, photonV0, int, V0Cores, "_PhotonV0"); //! -DECLARE_SOA_INDEX_COLUMN_FULL(LambdaV0, lambdaV0, int, V0Cores, "_LambdaV0"); //! -} // namespace sigma0Core - // for real data namespace sigma0Core { +DECLARE_SOA_COLUMN(PhotonV0ID, photonV0ID, int); +DECLARE_SOA_COLUMN(LambdaV0ID, lambdaV0ID, int); DECLARE_SOA_COLUMN(X, x, float); DECLARE_SOA_COLUMN(Y, y, float); DECLARE_SOA_COLUMN(Z, z, float); @@ -163,6 +158,9 @@ DECLARE_SOA_DYNAMIC_COLUMN(LambdaPhi, lambdaPhi, //! Phi in the range [0, 2pi) } // namespace sigma0Core DECLARE_SOA_TABLE(Sigma0Cores, "AOD", "SIGMA0CORES", + // Indices for debug + sigma0Core::PhotonV0ID, sigma0Core::LambdaV0ID, + // Basic properties sigma0Core::X, sigma0Core::Y, sigma0Core::Z, sigma0Core::DCADaughters, sigma0Core::PhotonPx, sigma0Core::PhotonPy, sigma0Core::PhotonPz, sigma0Core::PhotonMass, @@ -510,8 +508,8 @@ namespace sigma0EMPhoton // REGULAR COLUMNS FOR SIGMA0EMPHOTON DECLARE_SOA_COLUMN(PhotonID, photonID, int); DECLARE_SOA_COLUMN(PhotonEnergy, photonEnergy, float); -DECLARE_SOA_COLUMN(PhotonEta, photonEta, float); -DECLARE_SOA_COLUMN(PhotonPhi, photonPhi, float); +DECLARE_SOA_COLUMN(PhotonEMCEta, photonEMCEta, float); +DECLARE_SOA_COLUMN(PhotonEMCPhi, photonEMCPhi, float); DECLARE_SOA_COLUMN(PhotonM02, photonM02, float); DECLARE_SOA_COLUMN(PhotonM20, photonM20, float); DECLARE_SOA_COLUMN(PhotonNCells, photonNCells, int); @@ -520,14 +518,15 @@ DECLARE_SOA_COLUMN(PhotonIsExotic, photonIsExotic, bool); DECLARE_SOA_COLUMN(PhotonDistToBad, photonDistToBad, float); DECLARE_SOA_COLUMN(PhotonNLM, photonNLM, int); DECLARE_SOA_COLUMN(PhotonDefinition, photonDefinition, int); +DECLARE_SOA_COLUMN(PhotonHasAssocTrk, photonHasAssocTrk, bool); } // namespace sigma0PhotonExtra -DECLARE_SOA_TABLE(sigma0EMPhotons, "AOD", "SIGMA0EMPHOTON", +DECLARE_SOA_TABLE(Sigma0EMPhotons, "AOD", "SIGMA0EMPHOTON", sigma0EMPhoton::PhotonID, sigma0EMPhoton::PhotonEnergy, - sigma0EMPhoton::PhotonEta, - sigma0EMPhoton::PhotonPhi, + sigma0EMPhoton::PhotonEMCEta, + sigma0EMPhoton::PhotonEMCPhi, sigma0EMPhoton::PhotonM02, sigma0EMPhoton::PhotonM20, sigma0EMPhoton::PhotonNCells, @@ -535,7 +534,8 @@ DECLARE_SOA_TABLE(sigma0EMPhotons, "AOD", "SIGMA0EMPHOTON", sigma0EMPhoton::PhotonIsExotic, sigma0EMPhoton::PhotonDistToBad, sigma0EMPhoton::PhotonNLM, - sigma0EMPhoton::PhotonDefinition); + sigma0EMPhoton::PhotonDefinition, + sigma0EMPhoton::PhotonHasAssocTrk); // For Lambda extra info @@ -603,6 +603,7 @@ DECLARE_SOA_TABLE(Sigma0LambdaExtras, "AOD", "SIGMA0LAMBDA", // for MC namespace sigma0MCCore { +DECLARE_SOA_COLUMN(ParticleIdMC, particleIdMC, int); //! V0 Particle ID DECLARE_SOA_COLUMN(MCradius, mcradius, float); DECLARE_SOA_COLUMN(PDGCode, pdgCode, int); DECLARE_SOA_COLUMN(PDGCodeMother, pdgCodeMother, int); @@ -612,6 +613,9 @@ DECLARE_SOA_COLUMN(IsProducedByGenerator, isProducedByGenerator, bool); DECLARE_SOA_COLUMN(PhotonMCPx, photonmcpx, float); DECLARE_SOA_COLUMN(PhotonMCPy, photonmcpy, float); DECLARE_SOA_COLUMN(PhotonMCPz, photonmcpz, float); +DECLARE_SOA_COLUMN(PhotonAmplitudeA, photonAmplitudeA, float); // Energy fraction deposited by a particle inside this calo cell. +DECLARE_SOA_COLUMN(PhotonPDGCodePos, photonPDGCodePos, int); +DECLARE_SOA_COLUMN(PhotonPDGCodeNeg, photonPDGCodeNeg, int); DECLARE_SOA_COLUMN(IsPhotonPrimary, isPhotonPrimary, bool); DECLARE_SOA_COLUMN(PhotonPDGCode, photonPDGCode, int); DECLARE_SOA_COLUMN(PhotonPDGCodeMother, photonPDGCodeMother, int); @@ -621,6 +625,8 @@ DECLARE_SOA_COLUMN(LambdaMCPx, lambdamcpx, float); DECLARE_SOA_COLUMN(LambdaMCPy, lambdamcpy, float); DECLARE_SOA_COLUMN(LambdaMCPz, lambdamcpz, float); DECLARE_SOA_COLUMN(IsLambdaPrimary, isLambdaPrimary, bool); +DECLARE_SOA_COLUMN(LambdaPDGCodePos, lambdaPDGCodePos, int); +DECLARE_SOA_COLUMN(LambdaPDGCodeNeg, lambdaPDGCodeNeg, int); DECLARE_SOA_COLUMN(LambdaPDGCode, lambdaPDGCode, int); DECLARE_SOA_COLUMN(LambdaPDGCodeMother, lambdaPDGCodeMother, int); DECLARE_SOA_COLUMN(LambdaIsCorrectlyAssoc, lambdaIsCorrectlyAssoc, bool); @@ -727,15 +733,17 @@ DECLARE_SOA_DYNAMIC_COLUMN(LambdaMCPhi, lambdaMCPhi, //! Phi in the range [0, 2p } // namespace sigma0MCCore DECLARE_SOA_TABLE(Sigma0MCCores, "AOD", "SIGMA0MCCORES", + // MC particle index for debug + sigma0MCCore::ParticleIdMC, // Basic properties sigma0MCCore::MCradius, sigma0MCCore::PDGCode, sigma0MCCore::PDGCodeMother, sigma0MCCore::MCprocess, sigma0MCCore::IsProducedByGenerator, - sigma0MCCore::PhotonMCPx, sigma0MCCore::PhotonMCPy, sigma0MCCore::PhotonMCPz, - sigma0MCCore::IsPhotonPrimary, sigma0MCCore::PhotonPDGCode, sigma0MCCore::PhotonPDGCodeMother, sigma0MCCore::PhotonIsCorrectlyAssoc, + sigma0MCCore::PhotonMCPx, sigma0MCCore::PhotonMCPy, sigma0MCCore::PhotonMCPz, sigma0MCCore::PhotonAmplitudeA, + sigma0MCCore::PhotonPDGCodePos, sigma0MCCore::PhotonPDGCodeNeg, sigma0MCCore::IsPhotonPrimary, sigma0MCCore::PhotonPDGCode, sigma0MCCore::PhotonPDGCodeMother, sigma0MCCore::PhotonIsCorrectlyAssoc, sigma0MCCore::LambdaMCPx, sigma0MCCore::LambdaMCPy, sigma0MCCore::LambdaMCPz, - sigma0MCCore::IsLambdaPrimary, sigma0MCCore::LambdaPDGCode, sigma0MCCore::LambdaPDGCodeMother, sigma0MCCore::LambdaIsCorrectlyAssoc, + sigma0MCCore::LambdaPDGCodePos, sigma0MCCore::LambdaPDGCodeNeg, sigma0MCCore::IsLambdaPrimary, sigma0MCCore::LambdaPDGCode, sigma0MCCore::LambdaPDGCodeMother, sigma0MCCore::LambdaIsCorrectlyAssoc, // Dynamic columns sigma0MCCore::IsSigma0, @@ -764,10 +772,6 @@ DECLARE_SOA_TABLE(Sigma0MCCores, "AOD", "SIGMA0MCCORES", sigma0MCCore::LambdaMCY, sigma0MCCore::LambdaMCPhi); -namespace sigma0MCCore -{ -DECLARE_SOA_INDEX_COLUMN(McParticle, mcParticle); //! MC particle for Sigma0 -} // for MC namespace kstarMCCore @@ -959,12 +963,6 @@ DECLARE_SOA_TABLE(KStarGens, "AOD", "KSTARGENS", DECLARE_SOA_TABLE(SigmaCollRef, "AOD", "SIGMACOLLREF", //! optional table to refer back to a collision o2::soa::Index<>, v0data::StraCollisionId); -DECLARE_SOA_TABLE(SigmaIndices, "AOD", "SIGMAINDEX", //! index table when using AO2Ds - o2::soa::Index<>, sigma0Core::PhotonV0Id, sigma0Core::LambdaV0Id, o2::soa::Marker<1>); - -DECLARE_SOA_TABLE(SigmaMCLabels, "AOD", "SIGMAMCLABEL", //! optional table to refer to mcparticles - o2::soa::Index<>, sigma0MCCore::McParticleId); - DECLARE_SOA_TABLE(SigmaGenCollRef, "AOD", "SIGMAGENCOLLREF", //! optional table to refer back to a collision o2::soa::Index<>, v0data::StraMCCollisionId); diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index ff11db32d73..ee1a1db440f 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -25,6 +25,7 @@ #include "PWGLF/DataModel/LFStrangenessPIDTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGJE/DataModel/EMCALClusters.h" +#include "PWGEM/PhotonMeson/Utils/MCUtilities.h" #include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" @@ -69,21 +70,15 @@ using V0TOFDerivedMCDatas = soa::Join; -// TODO: Should I filter the clusters? -// we need to filter the emcal clusters, since we only want clusters from one specific clusterizer algorithm -//using EMCalClusters = o2::soa::Filtered; - static const std::vector DirList = {"V0BeforeSel", "PhotonSel", "LambdaSel", "KShortSel"}; static const std::vector DirList2 = {"EMCalPhotonBeforeSel", "EMCalPhotonSel"}; - struct sigma0builder { Service ccdb; ctpRateFetcher rateFetcher; //___________________________________________________ // KStar Specific - Produces kstarcores; // kstar candidates info for analysis Produces kshortExtras; // lambdas from sigma0 candidates info Produces kstarPhotonExtras; // photons from kstar candidates info @@ -94,24 +89,22 @@ struct sigma0builder { //__________________________________________________ // Sigma0 specific - Produces sigmaIndices; // references V0Cores from sigma0Gens Produces sigma0cores; // sigma0 candidates info for analysis Produces sigmaPhotonExtras; // photons from sigma0 candidates info - Produces sigmaEmCalPhotonExtras; // EMCAL photons from sigma0 candidates info + Produces sigmaEmCalPhotonExtras; // EMCAL photons from sigma0 candidates info Produces sigmaLambdaExtras; // lambdas from sigma0 candidates info Produces sigma0CollRefs; // references collisions from Sigma0Cores - Produces sigma0mccores; // Reco sigma0 MC properties - Produces sigma0mclabel; // Link of reco sigma0 to mcparticles + Produces sigma0mccores; // Reco sigma0 MC properties Produces sigma0Gens; // Generated sigma0s Produces sigma0GenCollRefs; // references collisions from sigma0Gens //__________________________________________________ // Pi0 specific - Produces pi0cores; // pi0 candidates info for analysis - Produces pi0coresRefs; // references collisions from photonpair - Produces pi0coresmc; // Reco pi0 MC properties - Produces pi0Gens; // Generated pi0s - Produces pi0GenCollRefs; // references collisions from pi0Gens + Produces pi0cores; // pi0 candidates info for analysis + Produces pi0coresRefs; // references collisions from photonpair + Produces pi0coresmc; // Reco pi0 MC properties + Produces pi0Gens; // Generated pi0s + Produces pi0GenCollRefs; // references collisions from pi0Gens //__________________________________________________ // pack track quality but separte also afterburner @@ -352,10 +345,6 @@ struct sigma0builder { ConfigurableAxis axisClrTime{"axisClrTime", {300, -30.0, 30.0}, "cluster time (ns)"}; ConfigurableAxis axisClrShape{"axisClrShape", {100, 0.0, 1.0}, "cluster shape"}; - // For manual sliceBy (necessary to calculate the correction factors) - // TODO: remove me! - PresliceUnsorted> perMcCollision = aod::v0data::straMCCollisionId; - void init(InitContext const&) { LOGF(info, "Initializing now: cross-checking correctness..."); @@ -413,7 +402,7 @@ struct sigma0builder { } } - bool fUsePCMPhoton = !doprocessRealDataWithEMCal && !doprocessMonteCarloWithEMCal; + bool fUsePCMPhoton = !doprocessRealDataWithEMCal && !doprocessMonteCarloWithEMCal && !doprocessPCMVsEMCalQA; for (const auto& histodir : DirList) { if ((histodir == "V0BeforeSel" && !fFillNoSelV0Histos) || @@ -466,7 +455,7 @@ struct sigma0builder { histos.add(histodir + "/h3dV0XYZ", "h3dV0XYZ", kTH3D, {axisXY, axisXY, axisZ}); } - if (fUsePCMPhoton){ + if (fUsePCMPhoton || doprocessPCMVsEMCalQA){ histos.add("PhotonSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "Mass"); @@ -481,18 +470,6 @@ struct sigma0builder { histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(11, "Phi"); histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(12, "TPCCR"); histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(13, "TPC NSigma"); - } - else{ - for (const auto& histodir : DirList2) { - histos.add(histodir + "/hpT", "hpT", kTH1D, {axisPt}); - histos.add(histodir + "/hDefinition", "hDefinition", kTH1D, {axisClrDefinition}); - histos.add(histodir + "/hNCells", "hNCells", kTH1D, {axisClrNCells}); - histos.add(histodir + "/hEnergy", "hEnergy", kTH1D, {axisClrEnergy}); - histos.add(histodir + "/h2dEtaVsPhi", "h2dEtaVsPhi", kTH2D, {axisRapidity, axisPhi}); - histos.add(histodir + "/hTime", "hTime", kTH1D, {axisClrTime}); - histos.add(histodir + "/hExotic", "hExotic", kTH1D, {{2, -0.5f, 1.5f}}); - histos.add(histodir + "/hShape", "hShape", kTH1D, {axisClrShape}); - } histos.add("EMCalPhotonSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); @@ -504,6 +481,17 @@ struct sigma0builder { histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(7, "Exotic"); histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(8, "Shape"); + } + else{ + for (const auto& histodir : DirList2) { + histos.add(histodir + "/hDefinition", "hDefinition", kTH1D, {axisClrDefinition}); + histos.add(histodir + "/h2dNCells", "h2dNCells", kTH2D, {axisPt, axisClrNCells}); + histos.add(histodir + "/h2dEnergy", "h2dEnergy", kTH2D, {axisPt, axisClrEnergy}); + histos.add(histodir + "/h2dEtaVsPhi", "h2dEtaVsPhi", kTH2D, {axisRapidity, axisPhi}); + histos.add(histodir + "/h2dTime", "h2dTime", kTH2D, {axisPt, axisClrTime}); + histos.add(histodir + "/hExotic", "hExotic", kTH1D, {{2, -0.5f, 1.5f}}); + histos.add(histodir + "/h2dShape", "h2dShape", kTH2D, {axisPt, axisClrShape}); + } } histos.add("LambdaSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); @@ -544,7 +532,8 @@ struct sigma0builder { histos.get(HIST("SigmaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); histos.get(HIST("SigmaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "Sigma Mass Window"); histos.get(HIST("SigmaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(3, "Sigma Y Window"); - + + histos.add("SigmaSel/hSigmaMassBeforeSel", "hSigmaMassBeforeSel", kTH1F, {axisSigmaMass}); histos.add("SigmaSel/hSigmaMassSelected", "hSigmaMassSelected", kTH1F, {axisSigmaMass}); histos.add("KStarSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); @@ -596,6 +585,31 @@ struct sigma0builder { histos.add("MCQA/hNoV0MCCores", "hNoV0MCCores", kTH1D, {{4, -0.5f, +3.5f}}); } + if (doprocessPCMVsEMCalQA){ + histos.add("PhotonMCQA/hPCMPhotonMCpT", "hPCMPhotonMCpT", kTH1D, {axisPt}); + histos.add("PhotonMCQA/h2dPCMPhotonMCpTResolution", "h2dPCMPhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/hPCMSigma0PhotonMCpT", "hPCMSigma0PhotonMCpT", kTH1D, {axisPt}); + histos.add("PhotonMCQA/h2dPCMSigma0PhotonMCpTResolution", "h2dPCMSigma0PhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); + + histos.add("PhotonMCQA/hEMCalPhotonMCpT", "hEMCalPhotonMCpT", kTH1D, {axisPt}); + histos.add("PhotonMCQA/h2dEMCalPhotonMCpTResolution", "h2dEMCalPhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dEMCalPhotonMCEnergyResolution", "h2dEMCalPhotonMCEnergyResolution", kTH2D, {axisClrEnergy, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dEMCalPhotonMCEtaResolution", "h2dEMCalPhotonMCEtaResolution", kTH2D, {axisRapidity, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dEMCalPhotonMCPhiResolution", "h2dEMCalPhotonMCPhiResolution", kTH2D, {axisPhi, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dEMCalPhotonMCFractionEnergy", "h2dEMCalPhotonMCFractionEnergy", kTH2D, {axisPt, {100, -1.0f, 1.0f}}); + + histos.add("PhotonMCQA/hEMCalSigma0PhotonMCpT", "hEMCalSigma0PhotonMCpT", kTH1D, {axisPt}); + histos.add("PhotonMCQA/h2dEMCalSigma0PhotonMCpTResolution", "h2dEMCalSigma0PhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dEMCalSigma0PhotonMCEnergyResolution", "h2dEMCalSigma0PhotonMCEnergyResolution", kTH2D, {axisClrEnergy, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dEMCalSigma0PhotonMCEtaResolution", "h2dEMCalSigma0PhotonMCEtaResolution", kTH2D, {axisRapidity, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dEMCalSigma0PhotonMCPhiResolution", "h2dEMCalSigma0PhotonMCPhiResolution", kTH2D, {axisPhi, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dEMCalSigma0PhotonMCFractionEnergy", "h2dEMCalSigma0PhotonMCFractionEnergy", kTH2D, {axisPt, {100, -1.0f, 1.0f}}); + + histos.add("PhotonMCQA/hGenPhoton", "hGenPhoton", kTH1D, {axisPt}); + histos.add("PhotonMCQA/hGenSigma0Photon", "hGenSigma0Photon", kTH1D, {axisPt}); + } + + if (doprocessGeneratedRun3 && genSelections.doQA) { // Pi0s @@ -772,6 +786,10 @@ struct sigma0builder { bool fIsV01Primary = false; bool fIsV02Primary = false; bool fV0PairProducedByGenerator = false; + int V01PDGCodePos = 0; + int V02PDGCodePos = 0; + int V01PDGCodeNeg = 0; + int V02PDGCodeNeg = 0; int V01PDGCode = 0; int V02PDGCode = 0; int V01PDGCodeMother = 0; @@ -787,6 +805,7 @@ struct sigma0builder { float V02MCpy = -999.f; float V02MCpz = -999.f; float V0PairMCRadius = -999.f; + float EMCalClusterAmplitude = -999.f; }; // ______________________________________________________ @@ -898,6 +917,10 @@ struct sigma0builder { // MC association info MCinfo.fIsV01Primary = v01MC.isPhysicalPrimary(); MCinfo.fIsV02Primary = v02MC.isPhysicalPrimary(); + MCinfo.V01PDGCodePos = v01MC.pdgCodePositive(); + MCinfo.V01PDGCodeNeg = v01MC.pdgCodeNegative(); + MCinfo.V02PDGCodePos = v02MC.pdgCodePositive(); + MCinfo.V02PDGCodeNeg = v02MC.pdgCodeNegative(); MCinfo.V01PDGCode = v01MC.pdgCode(); MCinfo.V02PDGCode = v02MC.pdgCode(); MCinfo.V01PDGCodeMother = v01MC.pdgCodeMother(); @@ -999,89 +1022,168 @@ struct sigma0builder { return MCinfo; } - + template - V0PairMCInfo getClusterV0PairMCInfo(TEMCalCls const& cluster, TV0 const& v0, TCollision const& collision, TMCParticles const& mcparticles) + V0PairMCInfo getClusterV0PairMCInfo(TEMCalCls const& cluster, + TV0 const& v0, + TCollision const& collision, + TMCParticles const& mcparticles) { - // TODO: you need to improve/recheck this + // Output container V0PairMCInfo MCinfo; - // _________________________ - // V0-specific treatment: - if (!v0.has_v0MCCore()) { + // ============================================================ + // 1) --- V0 (Lambda/AntiLambda) MC information --- + // ============================================================ + + // Check if V0 has MC information + if (!v0.has_v0MCCore()) { return MCinfo; } - - auto v0MC = v0.template v0MCCore_as>(); - auto MCParticle_v0 = mcparticles.rawIteratorAt(v0MC.particleIdMC()); - auto const& MCMothersList_v0 = MCParticle_v0.template mothers_as(); - if (!MCMothersList_v0.empty()) - return MCinfo; - + auto v0MC = v0.template v0MCCore_as>(); + auto mcLambda = mcparticles.rawIteratorAt(v0MC.particleIdMC()); + + // Save basic Lambda MC info (always saved, independent of matching) + MCinfo.V02MCpx = v0MC.pxMC(); + MCinfo.V02MCpy = v0MC.pyMC(); + MCinfo.V02MCpz = v0MC.pzMC(); + MCinfo.fIsV02Primary = v0MC.isPhysicalPrimary(); + MCinfo.V02PDGCode = v0MC.pdgCode(); + MCinfo.V02PDGCodeMother = v0MC.pdgCodeMother(); + MCinfo.V02PDGCodePos = v0MC.pdgCodePositive(); + MCinfo.V02PDGCodeNeg = v0MC.pdgCodeNegative(); + + // Check correct MC collision assignment (if available) if (collision.has_straMCCollision()) { - auto MCCollision = collision.template straMCCollision_as>(); - MCinfo.fIsV02CorrectlyAssign = (v0MC.straMCCollisionId() == MCCollision.globalIndex()); + auto MCCollision = collision.template straMCCollision_as< + soa::Join>(); + MCinfo.fIsV02CorrectlyAssign = (v0MC.straMCCollisionId() == MCCollision.globalIndex()); } - // Getting Lambda's mother - auto const& MCMother_v0 = MCMothersList_v0.front(); // First mother + // Retrieve Lambda mothers + auto const& lambdaMothers = mcLambda.template mothers_as(); + if (lambdaMothers.empty()) { + // No ancestry -> cannot match to any parent + return MCinfo; + } - // _________________________ - // EMCalClusters-specific treatment: - - // We loop over MC contributors to the cluster - for (size_t i = 0; i < cluster.amplitudeA().size(); i++) { - - // Get MC particle - int clusterInducerId = cluster.mcParticleIds()[i]; - auto mcClusterInducer = mcparticles.iteratorAt(clusterInducerId); - - // Get MC Mother - auto const& MCMothersList_Cluster = mcClusterInducer.template mothers_as(); - - if (MCMothersList_Cluster.empty()) + // Assumption: first mother is the physical one + auto const& lambdaMother = lambdaMothers.front(); + int lambdaMotherIndex = lambdaMother.globalIndex(); + + // ============================================================ + // 2) --- EMCal cluster: loop over MC contributors --- + // ============================================================ + + int matchedPhotonId = -1; // MC photon candidate + int matchedMotherIndex = -1; // Common Sigma0 candidate + + // Fallback: sum of all contributor momenta (useful for resolution studies, perhaps?) + float sumPx = 0.f, sumPy = 0.f, sumPz = 0.f; + + // Loop over all MC contributors to the cluster + for (size_t i = 0; i < cluster.mcParticleIds().size(); i++) { + + int mcId = cluster.mcParticleIds()[i]; + auto mcPart = mcparticles.iteratorAt(mcId); + + // Accumulate total momentum (fallback strategy) + sumPx += mcPart.px(); + sumPy += mcPart.py(); + sumPz += mcPart.pz(); + + // ------------------------------------------------------------ + // Check 1: + // Does this contributor come from a Sigma0/AntiSigma0 somewhere in its ancestry? + // ------------------------------------------------------------ + int daughterId = aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcPart, mcparticles, std::vector{PDG_t::kSigma0, PDG_t::kSigma0Bar}); + + if (daughterId < 0) + continue; // Not from Sigma0 -> try next contributor + + auto mcPhoton = mcparticles.iteratorAt(daughterId); + + // Sanity check: are we getting the correct particles? + auto dummy = mcparticles.rawIteratorAt(daughterId); + if (mcPhoton.globalIndex() != dummy.globalIndex()) + LOGF(fatal, "The behave of rawIteratorAt != iteratorAt. Index %i != %i. Please check. Aborting.", mcPhoton.globalIndex(), dummy.globalIndex()); + + // Require true photon, please + if (mcPhoton.pdgCode() != PDG_t::kGamma) continue; - auto const& MCMother_cluster = MCMothersList_Cluster.front(); // First mother + // Get Sigma0 index from photon mother + auto mothers = mcPhoton.mothersIds(); + if (mothers.empty()) // No mothers? Weird + continue; - if (MCMother_cluster.globalIndex() == MCMother_v0.globalIndex()){ - // Is it the same mother? - MCinfo.fV0PairProducedByGenerator = MCMother_v0.producedByGenerator(); - MCinfo.V0PairPDGCode = MCMother_v0.pdgCode(); - MCinfo.V0PairMCProcess = MCMother_v0.getProcess(); - MCinfo.V0PairMCParticleID = MCMother_v0.globalIndex(); - MCinfo.V0PairMCRadius = std::hypot(MCMother_v0.vx(), MCMother_v0.vy()); // production position radius + int sigmaIndex = mothers[0]; + + // ------------------------------------------------------------ + // Check 2: + // Does this photon share the same mother as the Lambda? + // ------------------------------------------------------------ + if (sigmaIndex == lambdaMotherIndex) { + matchedPhotonId = daughterId; + matchedMotherIndex = sigmaIndex; + MCinfo.EMCalClusterAmplitude = cluster.amplitudeA()[i]; + break; // SUCCESS -> stop loop + } + } - auto const& v0pairmothers = MCMother_v0.template mothers_as(); // Get mothers - if (!v0pairmothers.empty()) { - auto& v0PairMother = v0pairmothers.front(); // V0Pair mother, V0s grandmother - MCinfo.V0PairPDGCodeMother = v0PairMother.pdgCode(); - } + // ============================================================ + // 3) --- SUCCESS: true Lambda–photon pair from same Sigma0 --- + // ============================================================ + if (matchedPhotonId >= 0 && matchedMotherIndex >= 0) { - // Basic kinematic info - MCinfo.V01MCpx = mcClusterInducer.px(); - MCinfo.V01MCpy = mcClusterInducer.py(); - MCinfo.V01MCpz = mcClusterInducer.pz(); + auto mcPhoton = mcparticles.iteratorAt(matchedPhotonId); + auto mcSigma = mcparticles.iteratorAt(matchedMotherIndex); + + // --- Pair (Sigma0) information + MCinfo.fV0PairProducedByGenerator = mcSigma.producedByGenerator(); + MCinfo.V0PairPDGCode = mcSigma.pdgCode(); + MCinfo.V0PairMCProcess = mcSigma.getProcess(); + MCinfo.V0PairMCParticleID = mcSigma.globalIndex(); + MCinfo.V0PairMCRadius = std::hypot(mcSigma.vx(), mcSigma.vy()); + + // Sigma0 mother (optional) + auto const& sigmaMothers = mcSigma.template mothers_as(); + if (!sigmaMothers.empty()) { + MCinfo.V0PairPDGCodeMother = sigmaMothers.front().pdgCode(); + } - MCinfo.V02MCpx = v0MC.pxMC(); - MCinfo.V02MCpy = v0MC.pyMC(); - MCinfo.V02MCpz = v0MC.pzMC(); + // --- Photon MC info + MCinfo.V01MCpx = mcPhoton.px(); + MCinfo.V01MCpy = mcPhoton.py(); + MCinfo.V01MCpz = mcPhoton.pz(); + MCinfo.fIsV01Primary = mcPhoton.isPhysicalPrimary(); + MCinfo.V01PDGCode = mcPhoton.pdgCode(); - // MC association info - MCinfo.fIsV01Primary = mcClusterInducer.isPhysicalPrimary(); - MCinfo.fIsV02Primary = v0MC.isPhysicalPrimary(); - MCinfo.V01PDGCode = mcClusterInducer.pdgCode(); - MCinfo.V02PDGCode = v0MC.pdgCode(); - MCinfo.V01PDGCodeMother = MCMother_cluster.pdgCode(); - MCinfo.V02PDGCodeMother = v0MC.pdgCodeMother(); + if (!mcPhoton.mothersIds().empty()) { + auto mcMother = mcparticles.iteratorAt(mcPhoton.mothersIds()[0]); + MCinfo.V01PDGCodeMother = mcMother.pdgCode(); + } - break; // we found what we were looking for! - } + return MCinfo; } - - return MCinfo; + + // ============================================================ + // 4) --- FAILURE: no true matching photon found --- + // ============================================================ + + // Strategy: + // - Keep Lambda MC info (already filled) + // - For cluster: + // use summed momentum of contributors (proxy for cluster truth) + // - Leave PDG / primary flags as default (dummy) + + MCinfo.V01MCpx = sumPx; + MCinfo.V01MCpy = sumPy; + MCinfo.V01MCpz = sumPz; + + return MCinfo; } // ______________________________________________________ @@ -1230,12 +1332,21 @@ struct sigma0builder { std::vector getListOfRecoCollIndices(TMCollisions const& mcCollisions, TCollisions const& collisions) { std::vector listBestCollisionIdx(mcCollisions.size()); - for (auto const& mcCollision : mcCollisions) { - auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex()); + + // Custom grouping + std::vector> groupedCollisions(mcCollisions.size()); + + for (const auto& coll : collisions) { + groupedCollisions[coll.straMCCollisionId()].push_back(coll.globalIndex()); + } + + for (auto const& mcCollision : mcCollisions) { int biggestNContribs = -1; int bestCollisionIndex = -1; - for (auto const& collision : groupedCollisions) { + for (size_t i = 0; i < groupedCollisions[mcCollision.globalIndex()].size(); i++) { // consider event selections in the recoed <-> gen collision association, for the denominator (or numerator) of the efficiency (or signal loss)? + auto collision = collisions.rawIteratorAt(groupedCollisions[mcCollision.globalIndex()][i]); + if (eventSelections.useEvtSelInDenomEff && eventSelections.fUseEventSelection) { if (!IsEventAccepted(collision, false)) continue; @@ -1259,6 +1370,14 @@ struct sigma0builder { void fillGeneratedEventProperties(TMCCollisions const& mcCollisions, TCollisions const& collisions) { std::vector listBestCollisionIdx(mcCollisions.size()); + + // Custom grouping + std::vector> groupedCollisions(mcCollisions.size()); + + for (const auto& coll : collisions) { + groupedCollisions[coll.straMCCollisionId()].push_back(coll.globalIndex()); + } + for (auto const& mcCollision : mcCollisions) { // Apply selections on MC collisions if (eventSelections.applyZVtxSelOnMCPV && std::abs(mcCollision.posZ()) > eventSelections.maxZVtxPosition) { @@ -1275,15 +1394,16 @@ struct sigma0builder { } histos.fill(HIST("V0QA/hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events*/); - - auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex()); + // Check if there is at least one of the reconstructed collisions associated to this MC collision // If so, we consider it bool atLeastOne = false; int biggestNContribs = -1; float centrality = 100.5f; int nCollisions = 0; - for (auto const& collision : groupedCollisions) { + for (size_t i = 0; i < groupedCollisions[mcCollision.globalIndex()].size(); i++) { + auto collision = collisions.rawIteratorAt(groupedCollisions[mcCollision.globalIndex()][i]); + if (eventSelections.fUseEventSelection) { if (!IsEventAccepted(collision, false)) continue; @@ -1661,15 +1781,15 @@ struct sigma0builder { // calculate pT for cluster assuming they are photons (so no mass) float gammapT = sqrt(cluster.energy() * cluster.energy()) / std::cosh(cluster.eta()); - - histos.fill(HIST(MainDir2[mode]) + HIST("/hpT"), gammapT); + histos.fill(HIST(MainDir2[mode]) + HIST("/hDefinition"), cluster.definition()); - histos.fill(HIST(MainDir2[mode]) + HIST("/hNCells"), cluster.nCells()); - histos.fill(HIST(MainDir2[mode]) + HIST("/hEnergy"), cluster.energy()); + histos.fill(HIST(MainDir2[mode]) + HIST("/h2dNCells"), gammapT, cluster.nCells()); + histos.fill(HIST(MainDir2[mode]) + HIST("/h2dEnergy"), gammapT, cluster.energy()); + histos.fill(HIST(MainDir2[mode]) + HIST("/h2dShape"), gammapT, cluster.m02()); histos.fill(HIST(MainDir2[mode]) + HIST("/h2dEtaVsPhi"), cluster.eta(), cluster.phi()); - histos.fill(HIST(MainDir2[mode]) + HIST("/hTime"), cluster.time()); + histos.fill(HIST(MainDir2[mode]) + HIST("/h2dTime"), gammapT, cluster.time()); histos.fill(HIST(MainDir2[mode]) + HIST("/hExotic"), cluster.isExotic()); - histos.fill(HIST(MainDir2[mode]) + HIST("/hShape"), cluster.m02()); + } // Function to fill QA histograms. mode = 0 (before selections, all v0s), 1 (photon candidates), 2 (lambda/alambda candidates) @@ -2143,11 +2263,15 @@ struct sigma0builder { bool buildPCMSigma0(TV0Object const& lambda, TV0Object const& gamma, TCollision const& collision, TMCParticles const& mcparticles) { //_______________________________________________ - // TODO: Add Same index rejection + // Same index rejection + if (gamma.globalIndex() == lambda.globalIndex()) + return false; // Checking if both V0s are made of the very same tracks if (gamma.posTrackExtraId() == lambda.posTrackExtraId() || - gamma.negTrackExtraId() == lambda.negTrackExtraId()) { // TODO: Add extra conditions + gamma.negTrackExtraId() == lambda.negTrackExtraId() || + gamma.posTrackExtraId() == lambda.negTrackExtraId() || + gamma.negTrackExtraId() == lambda.posTrackExtraId()) { return false; } @@ -2181,27 +2305,24 @@ struct sigma0builder { // Sigma0 topological info auto sigma0TopoInfo = propagateV0PairToDCA(gamma, lambda); - sigma0cores(sigma0TopoInfo.X, sigma0TopoInfo.Y, sigma0TopoInfo.Z, sigma0TopoInfo.DCADau, + sigma0cores(gamma.globalIndex(), lambda.globalIndex(), sigma0TopoInfo.X, sigma0TopoInfo.Y, sigma0TopoInfo.Z, sigma0TopoInfo.DCADau, gamma.px(), gamma.py(), gamma.pz(), gamma.mGamma(), lambda.px(), lambda.py(), lambda.pz(), lambda.mLambda(), lambda.mAntiLambda()); // MC properties if constexpr (requires { gamma.motherMCPartId(); lambda.motherMCPartId(); }) { auto sigma0MCInfo = getV0PairMCInfo(gamma, lambda, collision, mcparticles); - - // TODO: Include PDG of V0 Daughters - sigma0mccores(sigma0MCInfo.V0PairMCRadius, sigma0MCInfo.V0PairPDGCode, sigma0MCInfo.V0PairPDGCodeMother, sigma0MCInfo.V0PairMCProcess, sigma0MCInfo.fV0PairProducedByGenerator, - sigma0MCInfo.V01MCpx, sigma0MCInfo.V01MCpy, sigma0MCInfo.V01MCpz, - sigma0MCInfo.fIsV01Primary, sigma0MCInfo.V01PDGCode, sigma0MCInfo.V01PDGCodeMother, sigma0MCInfo.fIsV01CorrectlyAssign, + + sigma0mccores(sigma0MCInfo.V0PairMCParticleID, sigma0MCInfo.V0PairMCRadius, sigma0MCInfo.V0PairPDGCode, sigma0MCInfo.V0PairPDGCodeMother, sigma0MCInfo.V0PairMCProcess, sigma0MCInfo.fV0PairProducedByGenerator, + sigma0MCInfo.V01MCpx, sigma0MCInfo.V01MCpy, sigma0MCInfo.V01MCpz, sigma0MCInfo.EMCalClusterAmplitude, + sigma0MCInfo.V01PDGCodePos, sigma0MCInfo.V01PDGCodeNeg, sigma0MCInfo.fIsV01Primary, sigma0MCInfo.V01PDGCode, sigma0MCInfo.V01PDGCodeMother, sigma0MCInfo.fIsV01CorrectlyAssign, sigma0MCInfo.V02MCpx, sigma0MCInfo.V02MCpy, sigma0MCInfo.V02MCpz, - sigma0MCInfo.fIsV02Primary, sigma0MCInfo.V02PDGCode, sigma0MCInfo.V02PDGCodeMother, sigma0MCInfo.fIsV02CorrectlyAssign); - - sigma0mclabel(sigma0MCInfo.V0PairMCParticleID); // TODO: Add this in sigma0mccores + sigma0MCInfo.V02PDGCodePos, sigma0MCInfo.V02PDGCodeNeg, sigma0MCInfo.fIsV02Primary, sigma0MCInfo.V02PDGCode, sigma0MCInfo.V02PDGCodeMother, sigma0MCInfo.fIsV02CorrectlyAssign); + } // Sigma0s -> stracollisions link histos.fill(HIST("SigmaSel/hSigma0DauDeltaIndex"), gamma.globalIndex() - lambda.globalIndex()); - sigma0CollRefs(collision.globalIndex()); - sigmaIndices(gamma.globalIndex(), lambda.globalIndex()); // TODO: put these indices in sigma0core table + sigma0CollRefs(collision.globalIndex()); //_______________________________________________ // Photon extra properties @@ -2268,7 +2389,7 @@ struct sigma0builder { // Build sigma0 candidate template - bool buildEMCalSigma0(TV0Object const& lambda, TEMCalClsObject const& gamma, TCollision const& collision, TMCParticles const& mcparticles) + bool buildEMCalSigma0(TV0Object const& lambda, TEMCalClsObject const& gamma, TCollision const& collision, TMCParticles const& mcparticles, std::vector const& emcaltracksmatched) { // calculate pT for cluster assuming they are photons (so no mass) float gammapT = sqrt(gamma.energy() * gamma.energy()) / std::cosh(gamma.eta()); @@ -2284,18 +2405,14 @@ struct sigma0builder { std::array pVecLambda{lambda.px(), lambda.py(), lambda.pz()}; auto arrMom = std::array{pVecPhotons, pVecLambda}; - float sigmaMass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassLambda0}); - float sigmaY = -999.f; + float sigmaMass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassLambda0}); - // TODO: do proper calculation of rapidity if MC - sigmaY = RecoDecay::y(std::array{gammapx + lambda.px(), gammapy + lambda.py(), gammapz + lambda.pz()}, o2::constants::physics::MassSigma0); - - // if constexpr (requires { gamma.pxMC(); lambda.pxMC(); }) // If MC - // sigmaY = RecoDecay::y(std::array{gamma.pxMC() + lambda.pxMC(), gamma.pyMC() + lambda.pyMC(), gamma.pzMC() + lambda.pzMC()}, o2::constants::physics::MassSigma0); - // else // If DATA - // sigmaY = RecoDecay::y(std::array{gamma.px() + lambda.px(), gamma.py() + lambda.py(), gamma.pz() + lambda.pz()}, o2::constants::physics::MassSigma0); + // N.B. At this stage, we are only using the reconstructed rapidity (ideally with a very loose cut) + // A proper selection should be done in the sigmaanalysis + float sigmaY = RecoDecay::y(std::array{gammapx + lambda.px(), gammapy + lambda.py(), gammapz + lambda.pz()}, o2::constants::physics::MassSigma0); histos.fill(HIST("SigmaSel/hSelectionStatistics"), 1.); + histos.fill(HIST("SigmaSel/hSigmaMassBeforeSel"), sigmaMass); if (TMath::Abs(sigmaMass - o2::constants::physics::MassSigma0) > Sigma0Window) return false; @@ -2308,33 +2425,32 @@ struct sigma0builder { //_______________________________________________ // Calculate properties & Fill tables - sigma0cores(0.0f, 0.0f, 0.0f, 0.0f, // N.B: Filling with dummy values for now + sigma0cores(gamma.globalIndex(), lambda.globalIndex(), + 0.0f, 0.0f, 0.0f, 0.0f, // N.B: Filling with dummy values for now gammapx, gammapy, gammapz, 0.0f, lambda.px(), lambda.py(), lambda.pz(), lambda.mLambda(), lambda.mAntiLambda()); // MC properties - if constexpr (requires { gamma.mcParticleIds(); lambda.motherMCPartId(); }) { // TODO: Solve MC treatment + if constexpr (requires { gamma.mcParticleIds(); lambda.motherMCPartId(); }) { auto sigma0MCInfo = getClusterV0PairMCInfo(gamma, lambda, collision, mcparticles); - - sigma0mccores(sigma0MCInfo.V0PairMCRadius, sigma0MCInfo.V0PairPDGCode, sigma0MCInfo.V0PairPDGCodeMother, sigma0MCInfo.V0PairMCProcess, sigma0MCInfo.fV0PairProducedByGenerator, - sigma0MCInfo.V01MCpx, sigma0MCInfo.V01MCpy, sigma0MCInfo.V01MCpz, - sigma0MCInfo.fIsV01Primary, sigma0MCInfo.V01PDGCode, sigma0MCInfo.V01PDGCodeMother, sigma0MCInfo.fIsV01CorrectlyAssign, + + sigma0mccores(sigma0MCInfo.V0PairMCParticleID, sigma0MCInfo.V0PairMCRadius, sigma0MCInfo.V0PairPDGCode, sigma0MCInfo.V0PairPDGCodeMother, sigma0MCInfo.V0PairMCProcess, sigma0MCInfo.fV0PairProducedByGenerator, + sigma0MCInfo.V01MCpx, sigma0MCInfo.V01MCpy, sigma0MCInfo.V01MCpz, sigma0MCInfo.EMCalClusterAmplitude, + sigma0MCInfo.V01PDGCodePos, sigma0MCInfo.V01PDGCodeNeg, sigma0MCInfo.fIsV01Primary, sigma0MCInfo.V01PDGCode, sigma0MCInfo.V01PDGCodeMother, sigma0MCInfo.fIsV01CorrectlyAssign, sigma0MCInfo.V02MCpx, sigma0MCInfo.V02MCpy, sigma0MCInfo.V02MCpz, - sigma0MCInfo.fIsV02Primary, sigma0MCInfo.V02PDGCode, sigma0MCInfo.V02PDGCodeMother, sigma0MCInfo.fIsV02CorrectlyAssign); - - sigma0mclabel(sigma0MCInfo.V0PairMCParticleID); + sigma0MCInfo.V02PDGCodePos, sigma0MCInfo.V02PDGCodeNeg, sigma0MCInfo.fIsV02Primary, sigma0MCInfo.V02PDGCode, sigma0MCInfo.V02PDGCodeMother, sigma0MCInfo.fIsV02CorrectlyAssign); + } - sigma0CollRefs(collision.globalIndex()); - sigmaIndices(gamma.globalIndex(), lambda.globalIndex()); + sigma0CollRefs(collision.globalIndex()); //_______________________________________________ - // Photon extra properties - // TODO: include EMCal photons table in header! + // Photon extra properties + bool hasAssociatedTrack = emcaltracksmatched[gamma.globalIndex()]; sigmaEmCalPhotonExtras(gamma.id(), gamma.energy(), gamma.eta(), gamma.phi(), gamma.m02(), gamma.m20(), gamma.nCells(), gamma.time(), - gamma.isExotic(), gamma.distanceToBadChannel(), gamma.nlm(), gamma.definition()); + gamma.isExotic(), gamma.distanceToBadChannel(), gamma.nlm(), gamma.definition(), hasAssociatedTrack); //_______________________________________________ // Lambda extra properties @@ -2499,8 +2615,8 @@ struct sigma0builder { } // Process photon and lambda candidates to build sigma0 candidates - template - void dataProcess(TCollision const& collisions, TV0s const& fullV0s, TEMCal const& fullEMCalClusters, TMCParticles const& mcparticles) + template + void dataProcess(TCollision const& collisions, TV0s const& fullV0s, TEMCal const& fullEMCalClusters, TEMCalTracks const& emcaltracks, TMCParticles const& mcparticles) { //_______________________________________________ // Initial setup @@ -2513,16 +2629,24 @@ struct sigma0builder { // Custom grouping std::vector> v0grouped(collisions.size()); - std::vector> emclustersgrouped(collisions.size()); + std::vector> emclustersgrouped(collisions.size()); + std::vector emcaltracksgrouped; + // Grouping step: for (const auto& v0 : fullV0s) { v0grouped[v0.straCollisionId()].push_back(v0.globalIndex()); } if constexpr (soa::is_table) { + emcaltracksgrouped.resize(fullEMCalClusters.size(), false); for (const auto& cluster : fullEMCalClusters) { emclustersgrouped[cluster.collisionId()].push_back(cluster.globalIndex()); } + + // Mapping emccluster to matched tracks + for (const auto& calotrack : emcaltracks) { + emcaltracksgrouped[calotrack.emcalclusterId()] = true; + } } //_______________________________________________ @@ -2533,8 +2657,7 @@ struct sigma0builder { if (!IsEventAccepted(coll, true)) continue; } - - // TODO: Check! Is this correct? + if constexpr (soa::is_table) { if (emclustersgrouped[coll.globalIndex()].size()==0 && eventSelections.fSkipEmptyEMCal) continue; @@ -2598,9 +2721,10 @@ struct sigma0builder { analyzeV0CollAssoc(coll, fullV0s, bestLambdasArray, false); // Lambda-analysis } } - + //_______________________________________________ // Photon-V0 nested loop + int nSigma0Candidates = 0; for (size_t i = 0; i < bestGammasArray.size(); ++i) { //_______________________________________________ @@ -2612,12 +2736,14 @@ struct sigma0builder { // Building sigma0 candidate & filling tables if constexpr (soa::is_table) { // using EMCal photons auto gamma1 = fullEMCalClusters.rawIteratorAt(bestGammasArray[i]); - if (!buildEMCalSigma0(lambda, gamma1, coll, mcparticles)) continue; + if (!buildEMCalSigma0(lambda, gamma1, coll, mcparticles, emcaltracksgrouped)) continue; } else { // using PCM photons auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); if (!buildPCMSigma0(lambda, gamma1, coll, mcparticles)) continue; - } + } + + nSigma0Candidates++; } } @@ -2652,6 +2778,8 @@ struct sigma0builder { } } } + + LOGF(info, "N. photons: %i, N. lambdas: %i, expected pairs: %i, got: %i", bestGammasArray.size(), bestLambdasArray.size(), bestGammasArray.size()*bestLambdasArray.size(), nSigma0Candidates); } } @@ -2775,35 +2903,252 @@ struct sigma0builder { } } } + + template + void runPCMVsEMCalQA(TCollision const& collisions, TV0s const& fullV0s, TEMCal const& fullEMCalClusters, TEMCalTracks const& emcaltracks, TMCParticles const& mcparticles) + { + + // Custom grouping + std::vector> v0grouped(collisions.size()); + std::vector> emclustersgrouped(collisions.size()); + std::vector emcaltracksgrouped(collisions.size(), false); + + // Grouping step: + for (const auto& v0 : fullV0s) { + v0grouped[v0.straCollisionId()].push_back(v0.globalIndex()); + } + + for (const auto& cluster : fullEMCalClusters) { + emclustersgrouped[cluster.collisionId()].push_back(cluster.globalIndex()); + } + + // Mapping emccluster to matched tracks + for (const auto& calotrack : emcaltracks) { + emcaltracksgrouped[calotrack.emcalclusterId()] = true; + } + + //_______________________________________________ + // Collisions loop + for (const auto& coll : collisions) { + // Event selection + if (eventSelections.fUseEventSelection) { + if (!IsEventAccepted(coll, true)) + continue; + } + + float centrality = doPPAnalysis ? coll.centFT0M() : coll.centFT0C(); + histos.fill(HIST("hEventCentrality"), centrality); + + //_______________________________________________ + // V0s loop + for (size_t i = 0; i < v0grouped[coll.globalIndex()].size(); i++) { + auto v0 = fullV0s.rawIteratorAt(v0grouped[coll.globalIndex()][i]); + + if (processPhotonCandidate(v0)) { // selecting PCM photons + + // Check if V0 has MC information + if (!v0.has_v0MCCore()) { + continue; + } + + auto v0MC = v0.template v0MCCore_as>(); + auto mcv0Photon = mcparticles.rawIteratorAt(v0MC.particleIdMC()); + + if (mcv0Photon.pdgCode() != PDG_t::kGamma || !mcv0Photon.isPhysicalPrimary() || TMath::Abs(mcv0Photon.y()) > 0.5) + continue; + + // MC pT histo + histos.fill(HIST("PhotonMCQA/hPCMPhotonMCpT"), mcv0Photon.pt()); + + // pT resolution + histos.fill(HIST("PhotonMCQA/h2dPCMPhotonMCpTResolution"), mcv0Photon.pt(), 1 - v0.pt()/mcv0Photon.pt()); + + // photon from sigma0/asigma0 + auto const& v0photonMothers = mcv0Photon.template mothers_as(); + if (!v0photonMothers.empty()) { + // Assumption: first mother is the physical one + auto const& v0photonMother = v0photonMothers.front(); + if (TMath::Abs(v0photonMother.pdgCode()) == PDG_t::kSigma0) { // Sigma0 or ASigma0 + // For efficiency + histos.fill(HIST("PhotonMCQA/hPCMSigma0PhotonMCpT"), mcv0Photon.pt()); + + // pT resolution + histos.fill(HIST("PhotonMCQA/h2dPCMSigma0PhotonMCpTResolution"), mcv0Photon.pt(), 1 - v0.pt()/mcv0Photon.pt()); + + } + } + } + } + + // EMCal clusters loop + for (size_t icl = 0; icl < emclustersgrouped[coll.globalIndex()].size(); icl++) { + auto cluster = fullEMCalClusters.rawIteratorAt(emclustersgrouped[coll.globalIndex()][icl]); + + if (!processEMCalPhotonCandidate(cluster)) + continue; + + // ============================================================ + // --- Reco kinematics (assume photon mass = 0) + // ============================================================ + + float E = cluster.energy(); + float eta = cluster.eta(); + + float pt = E / std::cosh(eta); + + // ============================================================ + // --- MC matching + // ============================================================ + // Map to avoid double counting of photons in case of conversions (one cluster can be matched to both e+ and e-) + std::unordered_set uniquePhotonIds; + + // Loop over cluster contributors + for (size_t i = 0; i < cluster.mcParticleIds().size(); i++) { + + int mcId = cluster.mcParticleIds()[i]; + auto mcPart = mcparticles.iteratorAt(mcId); + + // ============================================================ + // Find photon ancestor (including conversions) + // ============================================================ + + int photonId = -1; + + // Case 1: particle itself is photon + if (mcPart.pdgCode() == PDG_t::kGamma) { + photonId = mcPart.globalIndex(); + } + else { + // Case 2: climb ancestry to find photon + int candidateId = + aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain( + mcPart, mcparticles, std::vector{22}); + + if (candidateId >= 0) + photonId = candidateId; + } + + if (photonId < 0) + continue; + + // ============================================================ + // Avoid double counting (conversion case!) + // ============================================================ + + if (uniquePhotonIds.find(photonId) != uniquePhotonIds.end()) + continue; + + uniquePhotonIds.insert(photonId); + + auto mcPhoton = mcparticles.iteratorAt(photonId); + + // ============================================================ + // Select TRUE + PRIMARY photons + // ============================================================ + + if (mcPhoton.pdgCode() != PDG_t::kGamma || !mcPhoton.isPhysicalPrimary() || TMath::Abs(mcPhoton.y()) > 0.5) + continue; + + // ============================================================ + // Fill QA histos + // ============================================================ + + // MC pT histo + histos.fill(HIST("PhotonMCQA/hEMCalPhotonMCpT"), mcPhoton.pt()); + + // pT resolution + histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCpTResolution"), mcPhoton.pt(), 1 - pt/mcPhoton.pt()); + + // Energy resolution + histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCEnergyResolution"), mcPhoton.e(), 1 - cluster.energy()/mcPhoton.e()); + + // Eta resolution + histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCEtaResolution"), mcPhoton.eta(), cluster.eta() - mcPhoton.eta()); + + // Phi resolution + histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCPhiResolution"), mcPhoton.phi(), cluster.phi() - mcPhoton.phi()); + + // Fraction of energy Vs MC pT + histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCFractionEnergy"), mcPhoton.pt(), cluster.amplitudeA()[i]); + + // photon from sigma0/asigma0 + auto const& photonMothers = mcPhoton.template mothers_as(); + if (!photonMothers.empty()) { + // Assumption: first mother is the physical one + auto const& photonMother = photonMothers.front(); + if (TMath::Abs(photonMother.pdgCode()) == PDG_t::kSigma0) { // Sigma0 or ASigma0 + // For efficiency + histos.fill(HIST("PhotonMCQA/hEMCalSigma0PhotonMCpT"), mcPhoton.pt()); + + // pT resolution + histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCpTResolution"), mcPhoton.pt(), 1 - pt/mcPhoton.pt()); + + // Energy resolution + histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCEnergyResolution"), mcPhoton.e(), 1 - cluster.energy()/mcPhoton.e()); + + // Eta resolution + histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCEtaResolution"), mcPhoton.eta(), cluster.eta() - mcPhoton.eta()); + + // Phi resolution + histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCPhiResolution"), mcPhoton.phi(), cluster.phi() - mcPhoton.phi()); + + // Fraction of energy Vs MC pT + histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCFractionEnergy"), mcPhoton.pt(), cluster.amplitudeA()[i]); + } + } + } + } + + } // end of collisions loop + + // Process MC generated photons + for (const auto& mcpart : mcparticles) { + if (mcpart.pdgCode() != PDG_t::kGamma || !mcpart.isPhysicalPrimary() || TMath::Abs(mcpart.y()) > 0.5) + continue; + + histos.fill(HIST("PhotonMCQA/hGenPhoton"), mcpart.pt()); + + // photon from sigma0/asigma0 + auto const& photonMothers = mcpart.template mothers_as(); + if (!photonMothers.empty()) { + // Assumption: first mother is the physical one + auto const& photonMother = photonMothers.front(); + if (TMath::Abs(photonMother.pdgCode()) == PDG_t::kSigma0) { // Sigma0 or ASigma0 + histos.fill(HIST("PhotonMCQA/hGenSigma0Photon"), mcpart.pt()); + } + } + } + } + // Sigma0 processing part void processRealData(soa::Join const& collisions, V0StandardDerivedDatas const& fullV0s, dauTracks const&) { - dataProcess(collisions, fullV0s, nullptr, nullptr); + dataProcess(collisions, fullV0s, nullptr, nullptr, nullptr); } void processRealDataWithTOF(soa::Join const& collisions, V0TOFStandardDerivedDatas const& fullV0s, dauTracks const&) { - dataProcess(collisions, fullV0s, nullptr, nullptr); + dataProcess(collisions, fullV0s, nullptr, nullptr, nullptr); } - void processRealDataWithEMCal(soa::Join const& collisions, V0StandardDerivedDatas const& fullV0s, dauTracks const&, aod::EMCALClusters const& fullEMCalClusters) + void processRealDataWithEMCal(soa::Join const& collisions, V0StandardDerivedDatas const& fullV0s, dauTracks const&, aod::EMCALClusters const& fullEMCalClusters, aod::EMCALMatchedTracks const& emcmatchedtracks) { - dataProcess(collisions, fullV0s, fullEMCalClusters, nullptr); + dataProcess(collisions, fullV0s, fullEMCalClusters, emcmatchedtracks, nullptr); } void processMonteCarlo(soa::Join const& collisions, V0DerivedMCDatas const& fullV0s, aod::McParticles const& mcParticles, dauTracks const&, aod::MotherMCParts const&, soa::Join const&, soa::Join const&) { - dataProcess(collisions, fullV0s, nullptr, mcParticles); + dataProcess(collisions, fullV0s, nullptr, nullptr, mcParticles); } void processMonteCarloWithTOF(soa::Join const& collisions, V0TOFDerivedMCDatas const& fullV0s, aod::McParticles const& mcParticles, dauTracks const&, aod::MotherMCParts const&, soa::Join const&, soa::Join const&) { - dataProcess(collisions, fullV0s, nullptr, mcParticles); + dataProcess(collisions, fullV0s, nullptr, nullptr, mcParticles); } - void processMonteCarloWithEMCal(soa::Join const& collisions, V0DerivedMCDatas const& fullV0s, aod::McParticles const& mcParticles, dauTracks const&, aod::MotherMCParts const&, soa::Join const&, soa::Join const&, EMCalMCClusters const& fullEMCalMCClusters) + void processMonteCarloWithEMCal(soa::Join const& collisions, V0DerivedMCDatas const& fullV0s, aod::McParticles const& mcParticles, dauTracks const&, aod::MotherMCParts const&, soa::Join const&, soa::Join const&, EMCalMCClusters const& fullEMCalMCClusters, aod::EMCALMatchedTracks const& emcmatchedtracks) { - dataProcess(collisions, fullV0s, fullEMCalMCClusters, mcParticles); + dataProcess(collisions, fullV0s, fullEMCalMCClusters, emcmatchedtracks, mcParticles); } void processGeneratedRun3(aod::McParticles const& mcParticles) @@ -2827,6 +3172,11 @@ struct sigma0builder { runGenV0QA(mcCollisions, V0MCCores, collisions); } + void processPCMVsEMCalQA(soa::Join const& collisions, V0DerivedMCDatas const& fullV0s, aod::McParticles const& mcParticles, dauTracks const&, aod::MotherMCParts const&, soa::Join const&, soa::Join const&, EMCalMCClusters const& fullEMCalMCClusters, aod::EMCALMatchedTracks const& emcmatchedtracks) + { + runPCMVsEMCalQA(collisions, fullV0s, fullEMCalMCClusters, emcmatchedtracks, mcParticles); + } + PROCESS_SWITCH(sigma0builder, processRealData, "process as if real data", true); PROCESS_SWITCH(sigma0builder, processRealDataWithTOF, "process as if real data, uses TOF PID info", false); PROCESS_SWITCH(sigma0builder, processRealDataWithEMCal, "process as if real data, uses EMCal clusters", false); @@ -2837,6 +3187,7 @@ struct sigma0builder { PROCESS_SWITCH(sigma0builder, processV0QA, "process QA of lambdas and photons", false); PROCESS_SWITCH(sigma0builder, processV0MCQA, "process QA of lambdas and photons", false); PROCESS_SWITCH(sigma0builder, processV0Generated, "process QA of gen lambdas and photons", false); + PROCESS_SWITCH(sigma0builder, processPCMVsEMCalQA, "process QA of PCM and EMCal photons", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx b/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx index 9c37b4a8b1b..89cb5982da1 100644 --- a/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx @@ -60,8 +60,11 @@ using namespace o2::framework; using namespace o2::framework::expressions; using std::array; -using MCSigma0s = soa::Join; using Sigma0s = soa::Join; +using MCSigma0s = soa::Join; +using Sigma0sWithEMCal = soa::Join; +using MCSigma0sWithEMCal = soa::Join; + static const std::vector PhotonSels = {"NoSel", "V0Type", "DCADauToPV", "DCADau", "DauTPCCR", "TPCNSigmaEl", "V0pT", @@ -88,10 +91,6 @@ struct sigmaanalysis { ctpRateFetcher rateFetcher; //__________________________________________________ - // For manual sliceBy - // SliceCache cache; - PresliceUnsorted> perMcCollision = aod::v0data::straMCCollisionId; - HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; // Event level @@ -220,6 +219,23 @@ struct sigmaanalysis { Configurable PhotonPhiMax2{"PhotonPhiMax2", -1, "Phi min value to reject photons, region 2 (leave negative if no selection desired)"}; } photonSelections; + //// Photon criteria: + struct : ConfigurableGroup { + std::string prefix = "EMCalPhotonSelections"; // JSON group name + Configurable definition{"definition", 13, "Cluster definitions to be accepted (e.g. 13 for kV3MostSplitLowSeed)"}; + Configurable MinCells{"MinCells", 1, "Min number of cells in cluster"}; + Configurable MinEnergy{"MinEnergy", 0.0, "Minimum energy of selected clusters (GeV)"}; + Configurable MaxEnergy{"MaxEnergy", 5, "Max energy of selected clusters (GeV)"}; + Configurable MaxEta{"MaxEta", 1.0, "Max absolute cluster Eta"}; + Configurable MinTime{"MinTime", -50, "Minimum time of selected clusters (ns)"}; + Configurable MaxTime{"MaxTime", 50, "Max time of selected clusters (ns)"}; + Configurable RemoveExotic{"RemoveExotic", false, "Flag to enable the removal of exotic clusters"}; + Configurable MinM02{"MinM02", -1., "Minimum shower shape long axis"}; + Configurable MaxM02{"MaxM02", 5., "Max shower shape long axis"}; + Configurable RemoveMatchedTrack{"RemoveMatchedTrack", true, "Flag to enable the removal of clusters matched to tracks"}; + + } EMCalPhotonSelections; + struct : ConfigurableGroup { std::string prefix = "sigma0Selections"; // JSON group name Configurable Sigma0MinRapidity{"Sigma0MinRapidity", -0.5, "sigma0 min rapidity"}; @@ -277,6 +293,13 @@ struct sigmaanalysis { ConfigurableAxis axisPhi{"axisPhi", {200, 0, 2 * o2::constants::math::PI}, "Phi for photons"}; ConfigurableAxis axisZ{"axisZ", {120, -120.0f, 120.0f}, "V0 Z position (cm)"}; + // EMCal-specifc + ConfigurableAxis axisClrDefinition{"axisClrDefinition", {51, -0.5, 50.5}, "Cluster Definition"}; + ConfigurableAxis axisClrNCells{"axisClrNCells", {25, 0.0, 25}, "N cells per cluster"}; + ConfigurableAxis axisClrEnergy{"axisClrEnergy", {400, 0.0, 10}, "Energy per cluster"}; + ConfigurableAxis axisClrTime{"axisClrTime", {300, -30.0, 30.0}, "cluster time (ns)"}; + ConfigurableAxis axisClrShape{"axisClrShape", {100, 0.0, 1.0}, "cluster shape"}; + ConfigurableAxis axisCandSel{"axisCandSel", {20, 0.5f, +20.5f}, "Candidate Selection"}; // ML @@ -286,7 +309,7 @@ struct sigmaanalysis { void init(InitContext const&) { LOGF(info, "Initializing now: cross-checking correctness..."); - if ((doprocessRealData + doprocessMonteCarlo + doprocessPi0RealData + doprocessPi0MonteCarlo > 1) || + if ((doprocessRealData + doprocessRealDataWithEMCal + doprocessMonteCarlo + doprocessMonteCarloWithEMCal+ doprocessPi0RealData + doprocessPi0MonteCarlo > 1) || (doprocessGeneratedRun3 + doprocessPi0GeneratedRun3 > 1)) { LOGF(fatal, "You have enabled more than one process function. Please check your configuration! Aborting now."); } @@ -333,33 +356,53 @@ struct sigmaanalysis { histos.add("GeneralQA/hCentralityVsInteractionRate", "hCentralityVsInteractionRate", kTH2D, {axisCentrality, axisIRBinning}); } - if (doprocessRealData || doprocessMonteCarlo) { + if (doprocessRealData || doprocessRealDataWithEMCal || doprocessMonteCarlo || doprocessMonteCarloWithEMCal) { for (const auto& histodir : DirList) { if (fillSelhistos) { - histos.add(histodir + "/Photon/hTrackCode", "hTrackCode", kTH1D, {{11, 0.5f, 11.5f}}); - histos.add(histodir + "/Photon/hV0Type", "hV0Type", kTH1D, {{8, 0.5f, 8.5f}}); - histos.add(histodir + "/Photon/hDCANegToPV", "hDCANegToPV", kTH1D, {axisDCAtoPV}); - histos.add(histodir + "/Photon/hDCAPosToPV", "hDCAPosToPV", kTH1D, {axisDCAtoPV}); - histos.add(histodir + "/Photon/hDCADau", "hDCADau", kTH1D, {axisDCAdau}); - histos.add(histodir + "/Photon/hPosTPCCR", "hPosTPCCR", kTH1D, {axisTPCrows}); - histos.add(histodir + "/Photon/hNegTPCCR", "hNegTPCCR", kTH1D, {axisTPCrows}); - histos.add(histodir + "/Photon/hPosTPCNSigmaEl", "hPosTPCNSigmaEl", kTH1D, {axisTPCNSigma}); - histos.add(histodir + "/Photon/hNegTPCNSigmaEl", "hNegTPCNSigmaEl", kTH1D, {axisTPCNSigma}); - - histos.add(histodir + "/Photon/hpT", "hpT", kTH1D, {axisPt}); - histos.add(histodir + "/Photon/hY", "hY", kTH1D, {axisRapidity}); - histos.add(histodir + "/Photon/hPosEta", "hPosEta", kTH1D, {axisRapidity}); - histos.add(histodir + "/Photon/hNegEta", "hNegEta", kTH1D, {axisRapidity}); - histos.add(histodir + "/Photon/hRadius", "hRadius", kTH1D, {axisV0Radius}); - histos.add(histodir + "/Photon/hZ", "hZ", kTH1D, {axisZ}); - histos.add(histodir + "/Photon/h2dRZCut", "h2dRZCut", kTH2D, {axisZ, axisV0Radius}); - histos.add(histodir + "/Photon/h2dRZPlane", "h2dRZPlane", kTH2D, {axisZ, axisV0Radius}); - histos.add(histodir + "/Photon/hCosPA", "hCosPA", kTH1D, {axisCosPA}); - histos.add(histodir + "/Photon/hPsiPair", "hPsiPair", kTH1D, {axisPsiPair}); - histos.add(histodir + "/Photon/hPhi", "hPhi", kTH1D, {axisPhi}); - histos.add(histodir + "/Photon/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisPhotonMass}); - histos.add(histodir + "/Photon/hMass", "hMass", kTH1D, {axisPhotonMass}); + + if (doprocessRealData || doprocessMonteCarlo) { + + histos.add(histodir + "/Photon/hTrackCode", "hTrackCode", kTH1D, {{11, 0.5f, 11.5f}}); + histos.add(histodir + "/Photon/hV0Type", "hV0Type", kTH1D, {{8, 0.5f, 8.5f}}); + histos.add(histodir + "/Photon/hDCANegToPV", "hDCANegToPV", kTH1D, {axisDCAtoPV}); + histos.add(histodir + "/Photon/hDCAPosToPV", "hDCAPosToPV", kTH1D, {axisDCAtoPV}); + histos.add(histodir + "/Photon/hDCADau", "hDCADau", kTH1D, {axisDCAdau}); + histos.add(histodir + "/Photon/hPosTPCCR", "hPosTPCCR", kTH1D, {axisTPCrows}); + histos.add(histodir + "/Photon/hNegTPCCR", "hNegTPCCR", kTH1D, {axisTPCrows}); + histos.add(histodir + "/Photon/hPosTPCNSigmaEl", "hPosTPCNSigmaEl", kTH1D, {axisTPCNSigma}); + histos.add(histodir + "/Photon/hNegTPCNSigmaEl", "hNegTPCNSigmaEl", kTH1D, {axisTPCNSigma}); + + // PCM photons + histos.add(histodir + "/Photon/hpT", "hpT", kTH1D, {axisPt}); + histos.add(histodir + "/Photon/hY", "hY", kTH1D, {axisRapidity}); + histos.add(histodir + "/Photon/hPosEta", "hPosEta", kTH1D, {axisRapidity}); + histos.add(histodir + "/Photon/hNegEta", "hNegEta", kTH1D, {axisRapidity}); + histos.add(histodir + "/Photon/hRadius", "hRadius", kTH1D, {axisV0Radius}); + histos.add(histodir + "/Photon/hZ", "hZ", kTH1D, {axisZ}); + histos.add(histodir + "/Photon/h2dRZCut", "h2dRZCut", kTH2D, {axisZ, axisV0Radius}); + histos.add(histodir + "/Photon/h2dRZPlane", "h2dRZPlane", kTH2D, {axisZ, axisV0Radius}); + histos.add(histodir + "/Photon/hCosPA", "hCosPA", kTH1D, {axisCosPA}); + histos.add(histodir + "/Photon/hPsiPair", "hPsiPair", kTH1D, {axisPsiPair}); + histos.add(histodir + "/Photon/hPhi", "hPhi", kTH1D, {axisPhi}); + histos.add(histodir + "/Photon/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisPhotonMass}); + histos.add(histodir + "/Photon/hMass", "hMass", kTH1D, {axisPhotonMass}); + histos.add(histodir + "/Photon/h3dPhotonYSigma0Mass", "h3dPhotonYSigma0Mass", kTH3D, {axisRapidity, axisPt, axisSigmaMass}); + histos.add(histodir + "/Photon/h3dPhotonRadiusSigma0Mass", "h3dPhotonRadiusSigma0Mass", kTH3D, {axisV0Radius, axisPt, axisSigmaMass}); + + } + + if (doprocessRealDataWithEMCal || doprocessMonteCarloWithEMCal) { + // EMCal photons + histos.add(histodir + "/EMCalPhotonSel/hpT", "hpT", kTH1D, {axisPt}); + histos.add(histodir + "/EMCalPhotonSel/hDefinition", "hDefinition", kTH1D, {axisClrDefinition}); + histos.add(histodir + "/EMCalPhotonSel/hNCells", "hNCells", kTH1D, {axisClrNCells}); + histos.add(histodir + "/EMCalPhotonSel/hEnergy", "hEnergy", kTH1D, {axisClrEnergy}); + histos.add(histodir + "/EMCalPhotonSel/h2dEtaVsPhi", "h2dEtaVsPhi", kTH2D, {axisRapidity, axisPhi}); + histos.add(histodir + "/EMCalPhotonSel/hTime", "hTime", kTH1D, {axisClrTime}); + histos.add(histodir + "/EMCalPhotonSel/hExotic", "hExotic", kTH1D, {{2, -0.5f, 1.5f}}); + histos.add(histodir + "/EMCalPhotonSel/hShape", "hShape", kTH1D, {axisClrShape}); + } histos.add(histodir + "/Lambda/hTrackCode", "hTrackCode", kTH1D, {{11, 0.5f, 11.5f}}); histos.add(histodir + "/Lambda/hRadius", "hRadius", kTH1D, {axisV0Radius}); @@ -399,9 +442,7 @@ struct sigmaanalysis { histos.add(histodir + "/Sigma0/hRadius", "hRadius", kTH1D, {axisV0PairRadius}); histos.add(histodir + "/Sigma0/h2dRadiusVspT", "h2dRadiusVspT", kTH2D, {axisV0PairRadius, axisPt}); histos.add(histodir + "/Sigma0/hDCAPairDau", "hDCAPairDau", kTH1D, {axisDCAdau}); - histos.add(histodir + "/Sigma0/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisSigmaMass}); - histos.add(histodir + "/Sigma0/h3dPhotonYMass", "h3dPhotonYMass", kTH3D, {axisRapidity, axisPt, axisSigmaMass}); - histos.add(histodir + "/Sigma0/h3dPhotonRadiusMass", "h3dPhotonRadiusMass", kTH3D, {axisV0Radius, axisPt, axisSigmaMass}); + histos.add(histodir + "/Sigma0/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisSigmaMass}); histos.add(histodir + "/Sigma0/h3dOPAngleVsMass", "h3dOPAngleVsMass", kTH3D, {{140, 0.0f, +7.0f}, axisPt, axisSigmaMass}); histos.add(histodir + "/ASigma0/hMass", "hMass", kTH1D, {axisSigmaMass}); @@ -410,24 +451,23 @@ struct sigmaanalysis { histos.add(histodir + "/ASigma0/hRadius", "hRadius", kTH1D, {axisV0PairRadius}); histos.add(histodir + "/ASigma0/h2dRadiusVspT", "h2dRadiusVspT", kTH2D, {axisV0PairRadius, axisPt}); histos.add(histodir + "/ASigma0/hDCAPairDau", "hDCAPairDau", kTH1D, {axisDCAdau}); - histos.add(histodir + "/ASigma0/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisSigmaMass}); - histos.add(histodir + "/ASigma0/h3dPhotonYMass", "h3dPhotonYMass", kTH3D, {axisRapidity, axisPt, axisSigmaMass}); - histos.add(histodir + "/ASigma0/h3dPhotonRadiusMass", "h3dPhotonRadiusMass", kTH3D, {axisV0Radius, axisPt, axisSigmaMass}); + histos.add(histodir + "/ASigma0/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisSigmaMass}); histos.add(histodir + "/ASigma0/h3dOPAngleVsMass", "h3dOPAngleVsMass", kTH3D, {{140, 0.0f, +7.0f}, axisPt, axisSigmaMass}); // Process MC - if (doprocessMonteCarlo) { + if (doprocessMonteCarlo || doprocessMonteCarloWithEMCal) { if (fillSelhistos) { histos.add(histodir + "/MC/Photon/hV0ToCollAssoc", "hV0ToCollAssoc", kTH1D, {{2, 0.0f, 2.0f}}); histos.add(histodir + "/MC/Photon/hPt", "hPt", kTH1D, {axisPt}); - histos.add(histodir + "/MC/Photon/hMCPt", "hMCPt", kTH1D, {axisPt}); - histos.add(histodir + "/MC/Photon/hPosTPCNSigmaEl", "hPosTPCNSigmaEl", kTH1D, {axisTPCNSigma}); - histos.add(histodir + "/MC/Photon/hNegTPCNSigmaEl", "hNegTPCNSigmaEl", kTH1D, {axisTPCNSigma}); - histos.add(histodir + "/MC/Photon/h2dPAVsPt", "h2dPAVsPt", kTH2D, {axisPA, axisPt}); - histos.add(histodir + "/MC/Photon/hPt_BadCollAssig", "hPt_BadCollAssig", kTH1D, {axisPt}); - histos.add(histodir + "/MC/Photon/h2dPAVsPt_BadCollAssig", "h2dPAVsPt_BadCollAssig", kTH2D, {axisPA, axisPt}); + histos.add(histodir + "/MC/Photon/hMCPt", "hMCPt", kTH1D, {axisPt}); + if (doprocessMonteCarlo){ + histos.add(histodir + "/MC/Photon/h2dPAVsPt", "h2dPAVsPt", kTH2D, {axisPA, axisPt}); + histos.add(histodir + "/MC/Photon/hPt_BadCollAssig", "hPt_BadCollAssig", kTH1D, {axisPt}); + histos.add(histodir + "/MC/Photon/h2dPAVsPt_BadCollAssig", "h2dPAVsPt_BadCollAssig", kTH2D, {axisPA, axisPt}); + } + histos.add(histodir + "/MC/Lambda/hV0ToCollAssoc", "hV0ToCollAssoc", kTH1D, {{2, 0.0f, 2.0f}}); histos.add(histodir + "/MC/Lambda/hPt", "hPt", kTH1D, {axisPt}); histos.add(histodir + "/MC/Lambda/hMCPt", "hMCPt", kTH1D, {axisPt}); @@ -440,9 +480,7 @@ struct sigmaanalysis { histos.add(histodir + "/MC/ALambda/h3dTPCvsTOFNSigma_Pr", "h3dTPCvsTOFNSigma_Pr", kTH3D, {axisTPCNSigma, axisTOFNSigma, axisPt}); histos.add(histodir + "/MC/ALambda/h3dTPCvsTOFNSigma_Pi", "h3dTPCvsTOFNSigma_Pi", kTH3D, {axisTPCNSigma, axisTOFNSigma, axisPt}); } - - histos.add(histodir + "/MC/h2dArmenteros", "h2dArmenteros", kTH2D, {axisAPAlpha, axisAPQt}); - + histos.add(histodir + "/MC/Sigma0/hPt", "hPt", kTH1D, {axisPt}); histos.add(histodir + "/MC/Sigma0/hMCPt", "hMCPt", kTH1D, {axisPt}); histos.add(histodir + "/MC/Sigma0/hMass", "hMass", kTH1D, {axisSigmaMass}); @@ -465,11 +503,10 @@ struct sigmaanalysis { histos.add(histodir + "/MC/ASigma0/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisSigmaMass}); histos.add(histodir + "/MC/ASigma0/h3dMCProcess", "h3dMCProcess", kTH3D, {{50, -0.5f, 49.5f}, axisPt, axisSigmaMass}); - // 1/pT Resolution: + // pT Resolution: if (fillResoQAhistos && histodir == "BeforeSel") { - histos.add(histodir + "/MC/Reso/h3dGammaPtResoVsTPCCR", "h3dGammaPtResoVsTPCCR", kTH3D, {axisInvPt, axisDeltaPt, axisTPCrows}); - histos.add(histodir + "/MC/Reso/h3dGammaPtResoVsTPCCR", "h3dGammaPtResoVsTPCCR", kTH3D, {axisInvPt, axisDeltaPt, axisTPCrows}); - histos.add(histodir + "/MC/Reso/h2dGammaPtResolution", "h2dGammaPtResolution", kTH2D, {axisInvPt, axisDeltaPt}); + histos.add(histodir + "/MC/Reso/h2dGammaInvPtResolution", "h2dGammaPtResolution", kTH2D, {axisInvPt, axisDeltaPt}); + histos.add(histodir + "/MC/Reso/h2dGammaInvPtResolution", "h2dGammaInvPtResolution", kTH2D, {axisInvPt, axisDeltaPt}); histos.add(histodir + "/MC/Reso/h2dLambdaPtResolution", "h2dLambdaPtResolution", kTH2D, {axisInvPt, axisDeltaPt}); histos.add(histodir + "/MC/Reso/h3dLambdaPtResoVsTPCCR", "h3dLambdaPtResoVsTPCCR", kTH3D, {axisInvPt, axisDeltaPt, axisTPCrows}); histos.add(histodir + "/MC/Reso/h3dLambdaPtResoVsTPCCR", "h3dLambdaPtResoVsTPCCR", kTH3D, {axisInvPt, axisDeltaPt, axisTPCrows}); @@ -477,6 +514,7 @@ struct sigmaanalysis { histos.add(histodir + "/MC/Reso/h3dAntiLambdaPtResoVsTPCCR", "h3dAntiLambdaPtResoVsTPCCR", kTH3D, {axisInvPt, axisDeltaPt, axisTPCrows}); histos.add(histodir + "/MC/Reso/h3dAntiLambdaPtResoVsTPCCR", "h3dAntiLambdaPtResoVsTPCCR", kTH3D, {axisInvPt, axisDeltaPt, axisTPCrows}); histos.add(histodir + "/MC/Reso/h2dSigma0PtResolution", "h2dSigma0PtResolution", kTH2D, {axisInvPt, axisDeltaPt}); + histos.add(histodir + "/MC/Reso/h2dSigma0InvPtResolution", "h2dSigma0InvPtResolution", kTH2D, {axisInvPt, axisDeltaPt}); histos.add(histodir + "/MC/Reso/h2dAntiSigma0PtResolution", "h2dAntiSigma0PtResolution", kTH2D, {axisInvPt, axisDeltaPt}); histos.add(histodir + "/MC/Reso/h2dSigma0RadiusResolution", "h2dSigma0RadiusResolution", kTH2D, {axisPt, axisDeltaPt}); histos.add(histodir + "/MC/Reso/h2dASigma0RadiusResolution", "h2dASigma0RadiusResolution", kTH2D, {axisPt, axisDeltaPt}); @@ -751,12 +789,21 @@ struct sigmaanalysis { std::vector getListOfRecoCollIndices(TMCollisions const& mcCollisions, TCollisions const& collisions) { std::vector listBestCollisionIdx(mcCollisions.size()); - for (auto const& mcCollision : mcCollisions) { - auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex()); + + // Custom grouping + std::vector> groupedCollisions(mcCollisions.size()); + + for (const auto& coll : collisions) { + groupedCollisions[coll.straMCCollisionId()].push_back(coll.globalIndex()); + } + + for (auto const& mcCollision : mcCollisions) { int biggestNContribs = -1; int bestCollisionIndex = -1; - for (auto const& collision : groupedCollisions) { + for (size_t i = 0; i < groupedCollisions[mcCollision.globalIndex()].size(); i++) { // consider event selections in the recoed <-> gen collision association, for the denominator (or numerator) of the efficiency (or signal loss)? + auto collision = collisions.rawIteratorAt(groupedCollisions[mcCollision.globalIndex()][i]); + if (eventSelections.useEvtSelInDenomEff) { if (!IsEventAccepted(collision, false)) { continue; @@ -781,6 +828,14 @@ struct sigmaanalysis { void fillGeneratedEventProperties(TMCCollisions const& mcCollisions, TCollisions const& collisions) { std::vector listBestCollisionIdx(mcCollisions.size()); + + // Custom grouping + std::vector> groupedCollisions(mcCollisions.size()); + + for (const auto& coll : collisions) { + groupedCollisions[coll.straMCCollisionId()].push_back(coll.globalIndex()); + } + for (auto const& mcCollision : mcCollisions) { // Apply selections on MC collisions if (eventSelections.applyZVtxSelOnMCPV && std::abs(mcCollision.posZ()) > eventSelections.maxZVtxPosition) { @@ -797,15 +852,15 @@ struct sigmaanalysis { } histos.fill(HIST("Gen/hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events*/); - - auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex()); + // Check if there is at least one of the reconstructed collisions associated to this MC collision // If so, we consider it bool atLeastOne = false; int biggestNContribs = -1; float centrality = 100.5f; - int nCollisions = 0; - for (auto const& collision : groupedCollisions) { + int nCollisions = 0; + for (size_t i = 0; i < groupedCollisions[mcCollision.globalIndex()].size(); i++) { + auto collision = collisions.rawIteratorAt(groupedCollisions[mcCollision.globalIndex()][i]); if (!IsEventAccepted(collision, false)) { continue; @@ -926,14 +981,16 @@ struct sigmaanalysis { int TrkCode = 10; // 1: TPC-only, 2: TPC+Something, 3: ITS-Only, 4: ITS+TPC + Something, 10: anything else if (isGamma) { - if (sigma.photonPosTrackCode() == 1 && sigma.photonNegTrackCode() == 1) - TrkCode = 1; - if ((sigma.photonPosTrackCode() != 1 && sigma.photonNegTrackCode() == 1) || (sigma.photonPosTrackCode() == 1 && sigma.photonNegTrackCode() != 1)) - TrkCode = 2; - if (sigma.photonPosTrackCode() == 3 && sigma.photonNegTrackCode() == 3) - TrkCode = 3; - if (sigma.photonPosTrackCode() == 2 || sigma.photonNegTrackCode() == 2) - TrkCode = 4; + if constexpr (requires { sigma.photonV0Type();}) { + if (sigma.photonPosTrackCode() == 1 && sigma.photonNegTrackCode() == 1) + TrkCode = 1; + if ((sigma.photonPosTrackCode() != 1 && sigma.photonNegTrackCode() == 1) || (sigma.photonPosTrackCode() == 1 && sigma.photonNegTrackCode() != 1)) + TrkCode = 2; + if (sigma.photonPosTrackCode() == 3 && sigma.photonNegTrackCode() == 3) + TrkCode = 3; + if (sigma.photonPosTrackCode() == 2 || sigma.photonNegTrackCode() == 2) + TrkCode = 4; + } } else { if (sigma.lambdaPosTrackCode() == 1 && sigma.lambdaNegTrackCode() == 1) TrkCode = 1; @@ -955,10 +1012,9 @@ struct sigmaanalysis { //_______________________________________ // Gamma MC association if (sigma.photonPDGCode() == 22) { - if (sigma.photonmcpt() > 0) { - histos.fill(HIST("BeforeSel/MC/Reso/h3dGammaPtResoVsTPCCR"), 1.f / sigma.lambdamcpt(), 1.f / sigma.lambdaPt() - 1.f / sigma.lambdamcpt(), -1 * sigma.photonNegTPCCrossedRows()); // 1/pT resolution - histos.fill(HIST("BeforeSel/MC/Reso/h3dGammaPtResoVsTPCCR"), 1.f / sigma.lambdamcpt(), 1.f / sigma.lambdaPt() - 1.f / sigma.lambdamcpt(), sigma.photonPosTPCCrossedRows()); // 1/pT resolution - histos.fill(HIST("BeforeSel/MC/Reso/h2dGammaPtResolution"), 1.f / sigma.photonmcpt(), 1.f / sigma.photonPt() - 1.f / sigma.photonmcpt()); // pT resolution + if (sigma.photonmcpt() > 0) { + histos.fill(HIST("BeforeSel/MC/Reso/h2dGammaPtResolution"), sigma.photonmcpt(), sigma.photonPt() - sigma.photonmcpt()); // pT resolution + histos.fill(HIST("BeforeSel/MC/Reso/h2dGammaInvPtResolution"), 1.f / sigma.photonmcpt(), 1.f / sigma.photonPt() - 1.f / sigma.photonmcpt()); // pT resolution } } @@ -986,8 +1042,10 @@ struct sigmaanalysis { // Sigma and AntiSigma MC association if (sigma.isSigma0()) { histos.fill(HIST("BeforeSel/MC/Reso/h2dSigma0RadiusResolution"), sigma.mcpt(), sigma.radius() - sigma.mcradius()); // pT resolution - if (sigma.mcpt() > 0) - histos.fill(HIST("BeforeSel/MC/Reso/h2dSigma0PtResolution"), 1.f / sigma.mcpt(), 1.f / sigma.pt() - 1.f / sigma.mcpt()); // pT resolution + if (sigma.mcpt() > 0){ + histos.fill(HIST("BeforeSel/MC/Reso/h2dSigma0PtResolution"), sigma.mcpt(), sigma.pt() - sigma.mcpt()); // pT resolution + histos.fill(HIST("BeforeSel/MC/Reso/h2dSigma0InvPtResolution"), 1.f / sigma.mcpt(), 1.f / sigma.pt() - 1.f / sigma.mcpt()); // pT resolution + } } if (sigma.isAntiSigma0()) { histos.fill(HIST("BeforeSel/MC/Reso/h2dASigma0RadiusResolution"), sigma.mcpt(), sigma.radius() - sigma.mcradius()); // pT resolution @@ -1043,40 +1101,63 @@ struct sigmaanalysis { // Check whether it is before or after selections static constexpr std::string_view MainDir[] = {"BeforeSel", "AfterSel"}; - - // Get V0trackCode - int GammaTrkCode = retrieveV0TrackCode(sigma); - int LambdaTrkCode = retrieveV0TrackCode(sigma); - - float photonRZLineCut = TMath::Abs(sigma.photonZconv()) * TMath::Tan(2 * TMath::ATan(TMath::Exp(-photonSelections.PhotonMaxDauEta))) - photonSelections.PhotonLineCutZ0; + float centrality = getCentralityRun3(collision); if (fillSelhistos) { - //_______________________________________ - // Photon - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hTrackCode"), GammaTrkCode); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hV0Type"), sigma.photonV0Type()); - - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hDCANegToPV"), sigma.photonDCANegPV()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hDCAPosToPV"), sigma.photonDCAPosPV()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hDCADau"), sigma.photonDCADau()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPosTPCCR"), sigma.photonPosTPCCrossedRows()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hNegTPCCR"), sigma.photonNegTPCCrossedRows()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPosTPCNSigmaEl"), sigma.photonPosTPCNSigmaEl()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hNegTPCNSigmaEl"), sigma.photonNegTPCNSigmaEl()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hpT"), sigma.photonPt()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hY"), sigma.photonY()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPosEta"), sigma.photonPosEta()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hNegEta"), sigma.photonNegEta()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hRadius"), sigma.photonRadius()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hZ"), sigma.photonZconv()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h2dRZCut"), sigma.photonRadius(), photonRZLineCut); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h2dRZPlane"), sigma.photonZconv(), sigma.photonRadius()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hCosPA"), sigma.photonCosPA()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPsiPair"), sigma.photonPsiPair()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPhi"), sigma.photonPhi()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h3dMass"), centrality, sigma.photonPt(), sigma.photonMass()); - histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hMass"), sigma.photonMass()); + + // Get V0trackCode + int LambdaTrkCode = retrieveV0TrackCode(sigma); + + if constexpr (requires { sigma.photonV0Type();}) { // Processing PCM photon + + // Base properties + int GammaTrkCode = retrieveV0TrackCode(sigma); + float photonRZLineCut = TMath::Abs(sigma.photonZconv()) * TMath::Tan(2 * TMath::ATan(TMath::Exp(-photonSelections.PhotonMaxDauEta))) - photonSelections.PhotonLineCutZ0; + + //_______________________________________ + // Photon + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hTrackCode"), GammaTrkCode); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hV0Type"), sigma.photonV0Type()); + + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hDCANegToPV"), sigma.photonDCANegPV()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hDCAPosToPV"), sigma.photonDCAPosPV()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hDCADau"), sigma.photonDCADau()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPosTPCCR"), sigma.photonPosTPCCrossedRows()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hNegTPCCR"), sigma.photonNegTPCCrossedRows()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPosTPCNSigmaEl"), sigma.photonPosTPCNSigmaEl()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hNegTPCNSigmaEl"), sigma.photonNegTPCNSigmaEl()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hpT"), sigma.photonPt()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hY"), sigma.photonY()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPosEta"), sigma.photonPosEta()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hNegEta"), sigma.photonNegEta()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hRadius"), sigma.photonRadius()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hZ"), sigma.photonZconv()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h2dRZCut"), sigma.photonRadius(), photonRZLineCut); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h2dRZPlane"), sigma.photonZconv(), sigma.photonRadius()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hCosPA"), sigma.photonCosPA()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPsiPair"), sigma.photonPsiPair()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hPhi"), sigma.photonPhi()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h3dMass"), centrality, sigma.photonPt(), sigma.photonMass()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/hMass"), sigma.photonMass()); + + histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenteros"), sigma.photonAlpha(), sigma.photonQt()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h3dPhotonYSigma0Mass"), sigma.photonY(), sigma.pt(), sigma.sigma0Mass()); + histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h3dPhotonRadiusSigma0Mass"), sigma.photonRadius(), sigma.pt(), sigma.sigma0Mass()); + + } + + if constexpr (requires { sigma.photonDefinition();}) { // Processing EMCal photon + + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hpT"), sigma.photonPt()); + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hDefinition"), sigma.photonDefinition()); + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hNCells"), sigma.photonNCells()); + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hEnergy"), sigma.photonEnergy()); + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/h2dEtaVsPhi"), sigma.photonEMCEta(), sigma.photonEMCPhi()); + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hTime"), sigma.photonTime()); + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hExotic"), sigma.photonIsExotic()); + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hShape"), sigma.photonM02()); + } //_______________________________________ // Lambdas @@ -1094,13 +1175,12 @@ struct sigmaanalysis { histos.fill(HIST(MainDir[mode]) + HIST("/Lambda/hPosChi2PerNc"), sigma.lambdaPosChi2PerNcl()); histos.fill(HIST(MainDir[mode]) + HIST("/Lambda/hNegChi2PerNc"), sigma.lambdaNegChi2PerNcl()); histos.fill(HIST(MainDir[mode]) + HIST("/Lambda/hLifeTime"), sigma.lambdaLifeTime()); + + histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenteros"), sigma.lambdaAlpha(), sigma.lambdaQt()); } //_______________________________________ - // Sigmas and Lambdas - histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenteros"), sigma.photonAlpha(), sigma.photonQt()); - histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenteros"), sigma.lambdaAlpha(), sigma.lambdaQt()); - + // Sigmas and Lambdas if (sigma.lambdaAlpha() > 0) { if (fillSelhistos) { histos.fill(HIST(MainDir[mode]) + HIST("/Lambda/h2dTPCvsTOFNSigma_LambdaPr"), sigma.lambdaPosPrTPCNSigma(), sigma.lambdaPrTOFNSigma()); @@ -1118,9 +1198,7 @@ struct sigmaanalysis { histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/hRadius"), sigma.radius()); histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h2dRadiusVspT"), sigma.radius(), sigma.pt()); histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/hDCAPairDau"), sigma.dcadaughters()); - histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h3dMass"), centrality, sigma.pt(), sigma.sigma0Mass()); - histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h3dPhotonYMass"), sigma.photonY(), sigma.pt(), sigma.sigma0Mass()); - histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h3dPhotonRadiusMass"), sigma.photonRadius(), sigma.pt(), sigma.sigma0Mass()); + histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h3dMass"), centrality, sigma.pt(), sigma.sigma0Mass()); histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h3dOPAngleVsMass"), sigma.opAngle(), sigma.pt(), sigma.sigma0Mass()); } else { if (fillSelhistos) { @@ -1139,9 +1217,7 @@ struct sigmaanalysis { histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/hRadius"), sigma.radius()); histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h2dRadiusVspT"), sigma.radius(), sigma.pt()); histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/hDCAPairDau"), sigma.dcadaughters()); - histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h3dMass"), centrality, sigma.pt(), sigma.sigma0Mass()); - histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h3dPhotonYMass"), sigma.photonY(), sigma.pt(), sigma.sigma0Mass()); - histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h3dPhotonRadiusMass"), sigma.photonRadius(), sigma.pt(), sigma.sigma0Mass()); + histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h3dMass"), centrality, sigma.pt(), sigma.sigma0Mass()); histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h3dOPAngleVsMass"), sigma.opAngle(), sigma.pt(), sigma.sigma0Mass()); } @@ -1154,17 +1230,18 @@ struct sigmaanalysis { //_______________________________________ // Gamma MC association if (sigma.photonPDGCode() == 22) { + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hV0ToCollAssoc"), sigma.photonIsCorrectlyAssoc()); histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hPt"), sigma.photonPt()); - histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hMCPt"), sigma.photonmcpt()); - histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hPosTPCNSigmaEl"), sigma.photonPosTPCNSigmaEl()); - histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hNegTPCNSigmaEl"), sigma.photonNegTPCNSigmaEl()); - - histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/h2dPAVsPt"), TMath::ACos(sigma.photonCosPA()), sigma.photonmcpt()); - - if (!sigma.photonIsCorrectlyAssoc()) { - histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hPt_BadCollAssig"), sigma.photonmcpt()); - histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/h2dPAVsPt_BadCollAssig"), TMath::ACos(sigma.photonCosPA()), sigma.photonmcpt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hMCPt"), sigma.photonmcpt()); + + if constexpr (requires { sigma.photonV0Type();}) { // Processing PCM photon + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/h2dPAVsPt"), TMath::ACos(sigma.photonCosPA()), sigma.photonmcpt()); + + if (!sigma.photonIsCorrectlyAssoc()) { + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hPt_BadCollAssig"), sigma.photonmcpt()); + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/h2dPAVsPt_BadCollAssig"), TMath::ACos(sigma.photonCosPA()), sigma.photonmcpt()); + } } } @@ -1190,10 +1267,7 @@ struct sigmaanalysis { } //_______________________________________ // Sigma0 MC association - if (sigma.isSigma0()) { - histos.fill(HIST(MainDir[mode]) + HIST("/MC/h2dArmenteros"), sigma.photonAlpha(), sigma.photonQt()); - histos.fill(HIST(MainDir[mode]) + HIST("/MC/h2dArmenteros"), sigma.lambdaAlpha(), sigma.lambdaQt()); - + if (sigma.isSigma0()) { histos.fill(HIST(MainDir[mode]) + HIST("/MC/Sigma0/hPt"), sigma.pt()); histos.fill(HIST(MainDir[mode]) + HIST("/MC/Sigma0/hMCPt"), sigma.mcpt()); @@ -1211,9 +1285,6 @@ struct sigmaanalysis { //_______________________________________ // AntiSigma0 MC association if (sigma.isAntiSigma0()) { - histos.fill(HIST(MainDir[mode]) + HIST("/MC/h2dArmenteros"), sigma.photonAlpha(), sigma.photonQt()); - histos.fill(HIST(MainDir[mode]) + HIST("/MC/h2dArmenteros"), sigma.lambdaAlpha(), sigma.lambdaQt()); - histos.fill(HIST(MainDir[mode]) + HIST("/MC/ASigma0/hPt"), sigma.pt()); histos.fill(HIST(MainDir[mode]) + HIST("/MC/ASigma0/hMCPt"), sigma.mcpt()); @@ -1382,6 +1453,47 @@ struct sigmaanalysis { return true; } + // Apply specific selections for photons + template + bool selectEMCalPhoton(TClusterObject const& cand) + { + // Clusterizer + if (cand.photonDefinition() != EMCalPhotonSelections.definition && EMCalPhotonSelections.definition > -1) + return false; + + // Number of Cells + if (cand.photonNCells() < EMCalPhotonSelections.MinCells) + return false; + + // Energy + if (cand.photonEnergy() < EMCalPhotonSelections.MinEnergy || cand.photonEnergy() > EMCalPhotonSelections.MaxEnergy) + return false; + + // Eta + if (TMath::Abs(cand.photonEMCEta()) > EMCalPhotonSelections.MaxEta) + return false; + + // Timing + if (cand.photonTime() < EMCalPhotonSelections.MinTime || cand.photonTime() > EMCalPhotonSelections.MaxTime) + return false; + + // Exotic Clusters + if (cand.photonIsExotic() && EMCalPhotonSelections.RemoveExotic) { + return false; + } + + // Shower shape long axis + if (cand.photonNCells() > 1) { // Only if we have more than one + if (cand.photonM02() < EMCalPhotonSelections.MinM02 || cand.photonM02() > EMCalPhotonSelections.MaxM02) { + return false; + } + } + + // Has matched track? + if (cand.photonHasAssocTrk() && EMCalPhotonSelections.RemoveMatchedTrack) + return false; + } + // Apply specific selections for lambdas template bool selectLambda(TV0Object const& cand) @@ -1491,9 +1603,16 @@ struct sigmaanalysis { template bool processSigma0Candidate(TSigma0Object const& cand) { - // Photon specific selections - if (!selectPhoton(cand)) - return false; + // Photon specific selections + if constexpr (requires { cand.photonV0Type();}) { // Processing PCM photon + if (!selectPhoton(cand)) + return false; + } + + if constexpr (requires { cand.photonDefinition();}) { // Processing EMCal photon + if (!selectEMCalPhoton(cand)) + return false; + } // Lambda specific selections if (!selectLambda(cand)) @@ -1539,7 +1658,7 @@ struct sigmaanalysis { for (const auto& coll : collisions) { // Event selection - if (!IsEventAccepted(coll, true)) + if (!IsEventAccepted(coll, true)) // TODO: Should I Add event selection for events without EMCal? continue; // Sigma0s loop @@ -1559,7 +1678,7 @@ struct sigmaanalysis { fillHistos<0>(sigma0, coll); // Select sigma0 candidates - if (!processSigma0Candidate(sigma0)) + if (!processSigma0Candidate(sigma0)) continue; // Fill histos after all selections @@ -1721,11 +1840,21 @@ struct sigmaanalysis { analyzeRecoeSigma0s(collisions, fullSigma0s); } + void processRealDataWithEMCal(soa::Join const& collisions, Sigma0sWithEMCal const& fullSigma0s) + { + analyzeRecoeSigma0s(collisions, fullSigma0s); + } + void processMonteCarlo(soa::Join const& collisions, MCSigma0s const& fullSigma0s) { analyzeRecoeSigma0s(collisions, fullSigma0s); } + void processMonteCarloWithEMCal(soa::Join const& collisions, MCSigma0sWithEMCal const& fullSigma0s) + { + analyzeRecoeSigma0s(collisions, fullSigma0s); + } + // Simulated processing in Run 3 void processGeneratedRun3(soa::Join const& mcCollisions, soa::Join const& collisions, soa::Join const& Sigma0Gens) { @@ -1751,7 +1880,9 @@ struct sigmaanalysis { // _____________________________________________________ PROCESS_SWITCH(sigmaanalysis, processRealData, "Do real data analysis", true); + PROCESS_SWITCH(sigmaanalysis, processRealDataWithEMCal, "Do real data analysis with EMCal photons", false); PROCESS_SWITCH(sigmaanalysis, processMonteCarlo, "Do Monte-Carlo-based analysis", false); + PROCESS_SWITCH(sigmaanalysis, processMonteCarloWithEMCal, "Do Monte-Carlo-based analysis with EMCal photons", false); PROCESS_SWITCH(sigmaanalysis, processGeneratedRun3, "process MC generated Run 3", false); PROCESS_SWITCH(sigmaanalysis, processPi0RealData, "Do real data analysis for pi0 QA", false); PROCESS_SWITCH(sigmaanalysis, processPi0MonteCarlo, "Do Monte-Carlo-based analysis for pi0 QA", false); From 4ee76c9714badd726a060cd0f7296c3ee8f4fcb8 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Mon, 30 Mar 2026 14:02:53 +0000 Subject: [PATCH 3/4] Please consider the following formatting changes --- PWGLF/DataModel/LFSigmaTables.h | 11 +- .../Strangeness/sigma0builder.cxx | 300 +++++++++--------- PWGLF/Tasks/Strangeness/sigmaanalysis.cxx | 127 ++++---- 3 files changed, 213 insertions(+), 225 deletions(-) diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index 00633578ad4..6c8c9c3a232 100644 --- a/PWGLF/DataModel/LFSigmaTables.h +++ b/PWGLF/DataModel/LFSigmaTables.h @@ -500,7 +500,6 @@ DECLARE_SOA_TABLE(Sigma0PhotonExtras, "AOD", "SIGMA0PHOTON", sigma0PhotonExtra::PhotonNegTrackCode, sigma0PhotonExtra::PhotonV0Type); - // For EMCal Photon extra info namespace sigma0EMPhoton { @@ -520,7 +519,7 @@ DECLARE_SOA_COLUMN(PhotonNLM, photonNLM, int); DECLARE_SOA_COLUMN(PhotonDefinition, photonDefinition, int); DECLARE_SOA_COLUMN(PhotonHasAssocTrk, photonHasAssocTrk, bool); -} // namespace sigma0PhotonExtra +} // namespace sigma0EMPhoton DECLARE_SOA_TABLE(Sigma0EMPhotons, "AOD", "SIGMA0EMPHOTON", sigma0EMPhoton::PhotonID, @@ -537,7 +536,6 @@ DECLARE_SOA_TABLE(Sigma0EMPhotons, "AOD", "SIGMA0EMPHOTON", sigma0EMPhoton::PhotonDefinition, sigma0EMPhoton::PhotonHasAssocTrk); - // For Lambda extra info namespace sigma0LambdaExtra { @@ -603,7 +601,7 @@ DECLARE_SOA_TABLE(Sigma0LambdaExtras, "AOD", "SIGMA0LAMBDA", // for MC namespace sigma0MCCore { -DECLARE_SOA_COLUMN(ParticleIdMC, particleIdMC, int); //! V0 Particle ID +DECLARE_SOA_COLUMN(ParticleIdMC, particleIdMC, int); //! V0 Particle ID DECLARE_SOA_COLUMN(MCradius, mcradius, float); DECLARE_SOA_COLUMN(PDGCode, pdgCode, int); DECLARE_SOA_COLUMN(PDGCodeMother, pdgCodeMother, int); @@ -613,7 +611,7 @@ DECLARE_SOA_COLUMN(IsProducedByGenerator, isProducedByGenerator, bool); DECLARE_SOA_COLUMN(PhotonMCPx, photonmcpx, float); DECLARE_SOA_COLUMN(PhotonMCPy, photonmcpy, float); DECLARE_SOA_COLUMN(PhotonMCPz, photonmcpz, float); -DECLARE_SOA_COLUMN(PhotonAmplitudeA, photonAmplitudeA, float); // Energy fraction deposited by a particle inside this calo cell. +DECLARE_SOA_COLUMN(PhotonAmplitudeA, photonAmplitudeA, float); // Energy fraction deposited by a particle inside this calo cell. DECLARE_SOA_COLUMN(PhotonPDGCodePos, photonPDGCodePos, int); DECLARE_SOA_COLUMN(PhotonPDGCodeNeg, photonPDGCodeNeg, int); DECLARE_SOA_COLUMN(IsPhotonPrimary, isPhotonPrimary, bool); @@ -739,7 +737,7 @@ DECLARE_SOA_TABLE(Sigma0MCCores, "AOD", "SIGMA0MCCORES", // Basic properties sigma0MCCore::MCradius, sigma0MCCore::PDGCode, sigma0MCCore::PDGCodeMother, sigma0MCCore::MCprocess, sigma0MCCore::IsProducedByGenerator, - sigma0MCCore::PhotonMCPx, sigma0MCCore::PhotonMCPy, sigma0MCCore::PhotonMCPz, sigma0MCCore::PhotonAmplitudeA, + sigma0MCCore::PhotonMCPx, sigma0MCCore::PhotonMCPy, sigma0MCCore::PhotonMCPz, sigma0MCCore::PhotonAmplitudeA, sigma0MCCore::PhotonPDGCodePos, sigma0MCCore::PhotonPDGCodeNeg, sigma0MCCore::IsPhotonPrimary, sigma0MCCore::PhotonPDGCode, sigma0MCCore::PhotonPDGCodeMother, sigma0MCCore::PhotonIsCorrectlyAssoc, sigma0MCCore::LambdaMCPx, sigma0MCCore::LambdaMCPy, sigma0MCCore::LambdaMCPz, @@ -772,7 +770,6 @@ DECLARE_SOA_TABLE(Sigma0MCCores, "AOD", "SIGMA0MCCORES", sigma0MCCore::LambdaMCY, sigma0MCCore::LambdaMCPhi); - // for MC namespace kstarMCCore { diff --git a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx index ee1a1db440f..bad2b681d13 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -20,12 +20,12 @@ // gianni.shigeru.setoue.liveraro@cern.ch // +#include "PWGEM/PhotonMeson/Utils/MCUtilities.h" +#include "PWGJE/DataModel/EMCALClusters.h" #include "PWGLF/DataModel/LFSigmaTables.h" #include "PWGLF/DataModel/LFStrangenessMLTables.h" #include "PWGLF/DataModel/LFStrangenessPIDTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" -#include "PWGJE/DataModel/EMCALClusters.h" -#include "PWGEM/PhotonMeson/Utils/MCUtilities.h" #include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" @@ -68,7 +68,7 @@ using V0DerivedMCDatas = soa::Join; using V0TOFDerivedMCDatas = soa::Join; -using EMCalMCClusters = soa::Join; +using EMCalMCClusters = soa::Join; static const std::vector DirList = {"V0BeforeSel", "PhotonSel", "LambdaSel", "KShortSel"}; static const std::vector DirList2 = {"EMCalPhotonBeforeSel", "EMCalPhotonSel"}; @@ -94,17 +94,17 @@ struct sigma0builder { Produces sigmaEmCalPhotonExtras; // EMCAL photons from sigma0 candidates info Produces sigmaLambdaExtras; // lambdas from sigma0 candidates info Produces sigma0CollRefs; // references collisions from Sigma0Cores - Produces sigma0mccores; // Reco sigma0 MC properties + Produces sigma0mccores; // Reco sigma0 MC properties Produces sigma0Gens; // Generated sigma0s Produces sigma0GenCollRefs; // references collisions from sigma0Gens //__________________________________________________ // Pi0 specific - Produces pi0cores; // pi0 candidates info for analysis - Produces pi0coresRefs; // references collisions from photonpair - Produces pi0coresmc; // Reco pi0 MC properties - Produces pi0Gens; // Generated pi0s - Produces pi0GenCollRefs; // references collisions from pi0Gens + Produces pi0cores; // pi0 candidates info for analysis + Produces pi0coresRefs; // references collisions from photonpair + Produces pi0coresmc; // Reco pi0 MC properties + Produces pi0Gens; // Generated pi0s + Produces pi0GenCollRefs; // references collisions from pi0Gens //__________________________________________________ // pack track quality but separte also afterburner @@ -163,12 +163,12 @@ struct sigma0builder { Configurable minIR{"minIR", -1, "minimum IR collisions"}; Configurable maxIR{"maxIR", -1, "maximum IR collisions"}; - Configurable fSkipEmptyEMCal{"fSkipEmptyEMCal", true, "Flag to skip events without EMCal clusters"}; + Configurable fSkipEmptyEMCal{"fSkipEmptyEMCal", true, "Flag to skip events without EMCal clusters"}; } eventSelections; // Photon Source - //Configurable fUsePCMPhotons{"fUsePCMPhotons", true, "Use PCM Photons for sigma0/kstar reconstruction. If False, EMCal photons are used instead."}; + // Configurable fUsePCMPhotons{"fUsePCMPhotons", true, "Use PCM Photons for sigma0/kstar reconstruction. If False, EMCal photons are used instead."}; // Tables to fill Configurable fillPi0Tables{"fillPi0Tables", false, "fill pi0 tables for QA"}; @@ -243,15 +243,15 @@ struct sigma0builder { struct : ConfigurableGroup { std::string prefix = "EMCalPhotonSelections"; // JSON group name Configurable definition{"definition", 13, "Cluster definitions to be accepted (e.g. 13 for kV3MostSplitLowSeed)"}; - Configurable MinCells{"MinCells", 1, "Min number of cells in cluster"}; - Configurable MinEnergy{"MinEnergy", 0.0, "Minimum energy of selected clusters (GeV)"}; - Configurable MaxEnergy{"MaxEnergy", 5, "Max energy of selected clusters (GeV)"}; - Configurable MaxEta{"MaxEta", 1.0, "Max absolute cluster Eta"}; - Configurable MinTime{"MinTime", -50, "Minimum time of selected clusters (ns)"}; - Configurable MaxTime{"MaxTime", 50, "Max time of selected clusters (ns)"}; - Configurable RemoveExotic{"RemoveExotic", false, "Flag to enable the removal of exotic clusters"}; - Configurable MinM02{"MinM02", -1., "Minimum shower shape long axis"}; - Configurable MaxM02{"MaxM02", 5., "Max shower shape long axis"}; + Configurable MinCells{"MinCells", 1, "Min number of cells in cluster"}; + Configurable MinEnergy{"MinEnergy", 0.0, "Minimum energy of selected clusters (GeV)"}; + Configurable MaxEnergy{"MaxEnergy", 5, "Max energy of selected clusters (GeV)"}; + Configurable MaxEta{"MaxEta", 1.0, "Max absolute cluster Eta"}; + Configurable MinTime{"MinTime", -50, "Minimum time of selected clusters (ns)"}; + Configurable MaxTime{"MaxTime", 50, "Max time of selected clusters (ns)"}; + Configurable RemoveExotic{"RemoveExotic", false, "Flag to enable the removal of exotic clusters"}; + Configurable MinM02{"MinM02", -1., "Minimum shower shape long axis"}; + Configurable MaxM02{"MaxM02", 5., "Max shower shape long axis"}; } EMCalPhotonSelections; @@ -337,7 +337,7 @@ struct sigma0builder { ConfigurableAxis axisCandSel{"axisCandSel", {15, 0.5f, +15.5f}, "Candidate Selection"}; ConfigurableAxis axisIRBinning{"axisIRBinning", {151, -10, 1500}, "Binning for the interaction rate (kHz)"}; ConfigurableAxis axisLifetime{"axisLifetime", {200, 0, 50}, "Lifetime"}; - + // EMCal-specifc ConfigurableAxis axisClrDefinition{"axisClrDefinition", {51, -0.5, 50.5}, "Cluster Definition"}; ConfigurableAxis axisClrNCells{"axisClrNCells", {25, 0.0, 25}, "N cells per cluster"}; @@ -455,7 +455,7 @@ struct sigma0builder { histos.add(histodir + "/h3dV0XYZ", "h3dV0XYZ", kTH3D, {axisXY, axisXY, axisZ}); } - if (fUsePCMPhoton || doprocessPCMVsEMCalQA){ + if (fUsePCMPhoton || doprocessPCMVsEMCalQA) { histos.add("PhotonSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); histos.get(HIST("PhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "Mass"); @@ -481,19 +481,18 @@ struct sigma0builder { histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(7, "Exotic"); histos.get(HIST("EMCalPhotonSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(8, "Shape"); - } - else{ - for (const auto& histodir : DirList2) { + } else { + for (const auto& histodir : DirList2) { histos.add(histodir + "/hDefinition", "hDefinition", kTH1D, {axisClrDefinition}); histos.add(histodir + "/h2dNCells", "h2dNCells", kTH2D, {axisPt, axisClrNCells}); histos.add(histodir + "/h2dEnergy", "h2dEnergy", kTH2D, {axisPt, axisClrEnergy}); histos.add(histodir + "/h2dEtaVsPhi", "h2dEtaVsPhi", kTH2D, {axisRapidity, axisPhi}); histos.add(histodir + "/h2dTime", "h2dTime", kTH2D, {axisPt, axisClrTime}); histos.add(histodir + "/hExotic", "hExotic", kTH1D, {{2, -0.5f, 1.5f}}); - histos.add(histodir + "/h2dShape", "h2dShape", kTH2D, {axisPt, axisClrShape}); + histos.add(histodir + "/h2dShape", "h2dShape", kTH2D, {axisPt, axisClrShape}); } } - + histos.add("LambdaSel/hSelectionStatistics", "hSelectionStatistics", kTH1D, {axisCandSel}); histos.get(HIST("LambdaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); histos.get(HIST("LambdaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "Mass"); @@ -532,7 +531,7 @@ struct sigma0builder { histos.get(HIST("SigmaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); histos.get(HIST("SigmaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(2, "Sigma Mass Window"); histos.get(HIST("SigmaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(3, "Sigma Y Window"); - + histos.add("SigmaSel/hSigmaMassBeforeSel", "hSigmaMassBeforeSel", kTH1F, {axisSigmaMass}); histos.add("SigmaSel/hSigmaMassSelected", "hSigmaMassSelected", kTH1F, {axisSigmaMass}); @@ -585,14 +584,14 @@ struct sigma0builder { histos.add("MCQA/hNoV0MCCores", "hNoV0MCCores", kTH1D, {{4, -0.5f, +3.5f}}); } - if (doprocessPCMVsEMCalQA){ + if (doprocessPCMVsEMCalQA) { histos.add("PhotonMCQA/hPCMPhotonMCpT", "hPCMPhotonMCpT", kTH1D, {axisPt}); - histos.add("PhotonMCQA/h2dPCMPhotonMCpTResolution", "h2dPCMPhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dPCMPhotonMCpTResolution", "h2dPCMPhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); histos.add("PhotonMCQA/hPCMSigma0PhotonMCpT", "hPCMSigma0PhotonMCpT", kTH1D, {axisPt}); - histos.add("PhotonMCQA/h2dPCMSigma0PhotonMCpTResolution", "h2dPCMSigma0PhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dPCMSigma0PhotonMCpTResolution", "h2dPCMSigma0PhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); histos.add("PhotonMCQA/hEMCalPhotonMCpT", "hEMCalPhotonMCpT", kTH1D, {axisPt}); - histos.add("PhotonMCQA/h2dEMCalPhotonMCpTResolution", "h2dEMCalPhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); + histos.add("PhotonMCQA/h2dEMCalPhotonMCpTResolution", "h2dEMCalPhotonMCpTResolution", kTH2D, {axisPt, {100, -2.0f, 2.0f}}); histos.add("PhotonMCQA/h2dEMCalPhotonMCEnergyResolution", "h2dEMCalPhotonMCEnergyResolution", kTH2D, {axisClrEnergy, {100, -2.0f, 2.0f}}); histos.add("PhotonMCQA/h2dEMCalPhotonMCEtaResolution", "h2dEMCalPhotonMCEtaResolution", kTH2D, {axisRapidity, {100, -2.0f, 2.0f}}); histos.add("PhotonMCQA/h2dEMCalPhotonMCPhiResolution", "h2dEMCalPhotonMCPhiResolution", kTH2D, {axisPhi, {100, -2.0f, 2.0f}}); @@ -608,7 +607,6 @@ struct sigma0builder { histos.add("PhotonMCQA/hGenPhoton", "hGenPhoton", kTH1D, {axisPt}); histos.add("PhotonMCQA/hGenSigma0Photon", "hGenSigma0Photon", kTH1D, {axisPt}); } - if (doprocessGeneratedRun3 && genSelections.doQA) { @@ -786,7 +784,7 @@ struct sigma0builder { bool fIsV01Primary = false; bool fIsV02Primary = false; bool fV0PairProducedByGenerator = false; - int V01PDGCodePos = 0; + int V01PDGCodePos = 0; int V02PDGCodePos = 0; int V01PDGCodeNeg = 0; int V02PDGCodeNeg = 0; @@ -1022,7 +1020,7 @@ struct sigma0builder { return MCinfo; } - + template V0PairMCInfo getClusterV0PairMCInfo(TEMCalCls const& cluster, TV0 const& v0, @@ -1076,8 +1074,8 @@ struct sigma0builder { // 2) --- EMCal cluster: loop over MC contributors --- // ============================================================ - int matchedPhotonId = -1; // MC photon candidate - int matchedMotherIndex = -1; // Common Sigma0 candidate + int matchedPhotonId = -1; // MC photon candidate + int matchedMotherIndex = -1; // Common Sigma0 candidate // Fallback: sum of all contributor momenta (useful for resolution studies, perhaps?) float sumPx = 0.f, sumPy = 0.f, sumPz = 0.f; @@ -1097,16 +1095,16 @@ struct sigma0builder { // Check 1: // Does this contributor come from a Sigma0/AntiSigma0 somewhere in its ancestry? // ------------------------------------------------------------ - int daughterId = aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcPart, mcparticles, std::vector{PDG_t::kSigma0, PDG_t::kSigma0Bar}); + int daughterId = aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcPart, mcparticles, std::vector{PDG_t::kSigma0, PDG_t::kSigma0Bar}); if (daughterId < 0) continue; // Not from Sigma0 -> try next contributor - auto mcPhoton = mcparticles.iteratorAt(daughterId); + auto mcPhoton = mcparticles.iteratorAt(daughterId); // Sanity check: are we getting the correct particles? - auto dummy = mcparticles.rawIteratorAt(daughterId); - if (mcPhoton.globalIndex() != dummy.globalIndex()) + auto dummy = mcparticles.rawIteratorAt(daughterId); + if (mcPhoton.globalIndex() != dummy.globalIndex()) LOGF(fatal, "The behave of rawIteratorAt != iteratorAt. Index %i != %i. Please check. Aborting.", mcPhoton.globalIndex(), dummy.globalIndex()); // Require true photon, please @@ -1139,7 +1137,7 @@ struct sigma0builder { if (matchedPhotonId >= 0 && matchedMotherIndex >= 0) { auto mcPhoton = mcparticles.iteratorAt(matchedPhotonId); - auto mcSigma = mcparticles.iteratorAt(matchedMotherIndex); + auto mcSigma = mcparticles.iteratorAt(matchedMotherIndex); // --- Pair (Sigma0) information MCinfo.fV0PairProducedByGenerator = mcSigma.producedByGenerator(); @@ -1159,7 +1157,7 @@ struct sigma0builder { MCinfo.V01MCpy = mcPhoton.py(); MCinfo.V01MCpz = mcPhoton.pz(); MCinfo.fIsV01Primary = mcPhoton.isPhysicalPrimary(); - MCinfo.V01PDGCode = mcPhoton.pdgCode(); + MCinfo.V01PDGCode = mcPhoton.pdgCode(); if (!mcPhoton.mothersIds().empty()) { auto mcMother = mcparticles.iteratorAt(mcPhoton.mothersIds()[0]); @@ -1182,7 +1180,7 @@ struct sigma0builder { MCinfo.V01MCpx = sumPx; MCinfo.V01MCpy = sumPy; MCinfo.V01MCpz = sumPz; - + return MCinfo; } @@ -1340,7 +1338,7 @@ struct sigma0builder { groupedCollisions[coll.straMCCollisionId()].push_back(coll.globalIndex()); } - for (auto const& mcCollision : mcCollisions) { + for (auto const& mcCollision : mcCollisions) { int biggestNContribs = -1; int bestCollisionIndex = -1; for (size_t i = 0; i < groupedCollisions[mcCollision.globalIndex()].size(); i++) { @@ -1394,7 +1392,7 @@ struct sigma0builder { } histos.fill(HIST("V0QA/hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events*/); - + // Check if there is at least one of the reconstructed collisions associated to this MC collision // If so, we consider it bool atLeastOne = false; @@ -1772,7 +1770,7 @@ struct sigma0builder { } } - // Function to fill QA histograms. mode = 0 (before selections, all clusters), 1 after all selections + // Function to fill QA histograms. mode = 0 (before selections, all clusters), 1 after all selections template void fillEMCalHistos(TEMCalClusterObject const& cluster) { @@ -1781,15 +1779,14 @@ struct sigma0builder { // calculate pT for cluster assuming they are photons (so no mass) float gammapT = sqrt(cluster.energy() * cluster.energy()) / std::cosh(cluster.eta()); - + histos.fill(HIST(MainDir2[mode]) + HIST("/hDefinition"), cluster.definition()); histos.fill(HIST(MainDir2[mode]) + HIST("/h2dNCells"), gammapT, cluster.nCells()); histos.fill(HIST(MainDir2[mode]) + HIST("/h2dEnergy"), gammapT, cluster.energy()); histos.fill(HIST(MainDir2[mode]) + HIST("/h2dShape"), gammapT, cluster.m02()); - histos.fill(HIST(MainDir2[mode]) + HIST("/h2dEtaVsPhi"), cluster.eta(), cluster.phi()); + histos.fill(HIST(MainDir2[mode]) + HIST("/h2dEtaVsPhi"), cluster.eta(), cluster.phi()); histos.fill(HIST(MainDir2[mode]) + HIST("/h2dTime"), gammapT, cluster.time()); histos.fill(HIST(MainDir2[mode]) + HIST("/hExotic"), cluster.isExotic()); - } // Function to fill QA histograms. mode = 0 (before selections, all v0s), 1 (photon candidates), 2 (lambda/alambda candidates) @@ -1876,12 +1873,12 @@ struct sigma0builder { histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 3.); if (cluster.energy() < EMCalPhotonSelections.MinEnergy || cluster.energy() > EMCalPhotonSelections.MaxEnergy) return false; - + // Eta histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 4.); if (TMath::Abs(cluster.eta()) > EMCalPhotonSelections.MaxEta) return false; - + // Timing histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 5.); if (cluster.time() < EMCalPhotonSelections.MinTime || cluster.time() > EMCalPhotonSelections.MaxTime) @@ -1889,20 +1886,20 @@ struct sigma0builder { // Exotic Clusters histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 6.); - if (cluster.isExotic() && EMCalPhotonSelections.RemoveExotic) { + if (cluster.isExotic() && EMCalPhotonSelections.RemoveExotic) { return false; } // Shower shape long axis histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 7.); - if (cluster.nCells() > 1) { // Only if we have more than one - if (cluster.m02() < EMCalPhotonSelections.MinM02 || cluster.m02() > EMCalPhotonSelections.MaxM02) { + if (cluster.nCells() > 1) { // Only if we have more than one + if (cluster.m02() < EMCalPhotonSelections.MinM02 || cluster.m02() > EMCalPhotonSelections.MaxM02) { return false; } } histos.fill(HIST("EMCalPhotonSel/hSelectionStatistics"), 8.); - + return true; } @@ -2266,12 +2263,12 @@ struct sigma0builder { // Same index rejection if (gamma.globalIndex() == lambda.globalIndex()) return false; - + // Checking if both V0s are made of the very same tracks if (gamma.posTrackExtraId() == lambda.posTrackExtraId() || gamma.negTrackExtraId() == lambda.negTrackExtraId() || gamma.posTrackExtraId() == lambda.negTrackExtraId() || - gamma.negTrackExtraId() == lambda.posTrackExtraId()) { + gamma.negTrackExtraId() == lambda.posTrackExtraId()) { return false; } @@ -2311,18 +2308,17 @@ struct sigma0builder { // MC properties if constexpr (requires { gamma.motherMCPartId(); lambda.motherMCPartId(); }) { auto sigma0MCInfo = getV0PairMCInfo(gamma, lambda, collision, mcparticles); - + sigma0mccores(sigma0MCInfo.V0PairMCParticleID, sigma0MCInfo.V0PairMCRadius, sigma0MCInfo.V0PairPDGCode, sigma0MCInfo.V0PairPDGCodeMother, sigma0MCInfo.V0PairMCProcess, sigma0MCInfo.fV0PairProducedByGenerator, sigma0MCInfo.V01MCpx, sigma0MCInfo.V01MCpy, sigma0MCInfo.V01MCpz, sigma0MCInfo.EMCalClusterAmplitude, sigma0MCInfo.V01PDGCodePos, sigma0MCInfo.V01PDGCodeNeg, sigma0MCInfo.fIsV01Primary, sigma0MCInfo.V01PDGCode, sigma0MCInfo.V01PDGCodeMother, sigma0MCInfo.fIsV01CorrectlyAssign, sigma0MCInfo.V02MCpx, sigma0MCInfo.V02MCpy, sigma0MCInfo.V02MCpz, sigma0MCInfo.V02PDGCodePos, sigma0MCInfo.V02PDGCodeNeg, sigma0MCInfo.fIsV02Primary, sigma0MCInfo.V02PDGCode, sigma0MCInfo.V02PDGCodeMother, sigma0MCInfo.fIsV02CorrectlyAssign); - } // Sigma0s -> stracollisions link histos.fill(HIST("SigmaSel/hSigma0DauDeltaIndex"), gamma.globalIndex() - lambda.globalIndex()); - sigma0CollRefs(collision.globalIndex()); + sigma0CollRefs(collision.globalIndex()); //_______________________________________________ // Photon extra properties @@ -2390,7 +2386,7 @@ struct sigma0builder { // Build sigma0 candidate template bool buildEMCalSigma0(TV0Object const& lambda, TEMCalClsObject const& gamma, TCollision const& collision, TMCParticles const& mcparticles, std::vector const& emcaltracksmatched) - { + { // calculate pT for cluster assuming they are photons (so no mass) float gammapT = sqrt(gamma.energy() * gamma.energy()) / std::cosh(gamma.eta()); @@ -2405,9 +2401,9 @@ struct sigma0builder { std::array pVecLambda{lambda.px(), lambda.py(), lambda.pz()}; auto arrMom = std::array{pVecPhotons, pVecLambda}; - float sigmaMass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassLambda0}); + float sigmaMass = RecoDecay::m(arrMom, std::array{o2::constants::physics::MassPhoton, o2::constants::physics::MassLambda0}); - // N.B. At this stage, we are only using the reconstructed rapidity (ideally with a very loose cut) + // N.B. At this stage, we are only using the reconstructed rapidity (ideally with a very loose cut) // A proper selection should be done in the sigmaanalysis float sigmaY = RecoDecay::y(std::array{gammapx + lambda.px(), gammapy + lambda.py(), gammapz + lambda.pz()}, o2::constants::physics::MassSigma0); @@ -2422,36 +2418,35 @@ struct sigma0builder { histos.fill(HIST("SigmaSel/hSigmaMassSelected"), sigmaMass); histos.fill(HIST("SigmaSel/hSelectionStatistics"), 3.); - + //_______________________________________________ // Calculate properties & Fill tables - sigma0cores(gamma.globalIndex(), lambda.globalIndex(), + sigma0cores(gamma.globalIndex(), lambda.globalIndex(), 0.0f, 0.0f, 0.0f, 0.0f, // N.B: Filling with dummy values for now - gammapx, gammapy, gammapz, 0.0f, + gammapx, gammapy, gammapz, 0.0f, lambda.px(), lambda.py(), lambda.pz(), lambda.mLambda(), lambda.mAntiLambda()); // MC properties - if constexpr (requires { gamma.mcParticleIds(); lambda.motherMCPartId(); }) { - + if constexpr (requires { gamma.mcParticleIds(); lambda.motherMCPartId(); }) { + auto sigma0MCInfo = getClusterV0PairMCInfo(gamma, lambda, collision, mcparticles); - + sigma0mccores(sigma0MCInfo.V0PairMCParticleID, sigma0MCInfo.V0PairMCRadius, sigma0MCInfo.V0PairPDGCode, sigma0MCInfo.V0PairPDGCodeMother, sigma0MCInfo.V0PairMCProcess, sigma0MCInfo.fV0PairProducedByGenerator, sigma0MCInfo.V01MCpx, sigma0MCInfo.V01MCpy, sigma0MCInfo.V01MCpz, sigma0MCInfo.EMCalClusterAmplitude, sigma0MCInfo.V01PDGCodePos, sigma0MCInfo.V01PDGCodeNeg, sigma0MCInfo.fIsV01Primary, sigma0MCInfo.V01PDGCode, sigma0MCInfo.V01PDGCodeMother, sigma0MCInfo.fIsV01CorrectlyAssign, sigma0MCInfo.V02MCpx, sigma0MCInfo.V02MCpy, sigma0MCInfo.V02MCpz, sigma0MCInfo.V02PDGCodePos, sigma0MCInfo.V02PDGCodeNeg, sigma0MCInfo.fIsV02Primary, sigma0MCInfo.V02PDGCode, sigma0MCInfo.V02PDGCodeMother, sigma0MCInfo.fIsV02CorrectlyAssign); - } - sigma0CollRefs(collision.globalIndex()); + sigma0CollRefs(collision.globalIndex()); //_______________________________________________ - // Photon extra properties + // Photon extra properties bool hasAssociatedTrack = emcaltracksmatched[gamma.globalIndex()]; sigmaEmCalPhotonExtras(gamma.id(), gamma.energy(), gamma.eta(), gamma.phi(), - gamma.m02(), gamma.m20(), gamma.nCells(), gamma.time(), - gamma.isExotic(), gamma.distanceToBadChannel(), gamma.nlm(), gamma.definition(), hasAssociatedTrack); - + gamma.m02(), gamma.m20(), gamma.nCells(), gamma.time(), + gamma.isExotic(), gamma.distanceToBadChannel(), gamma.nlm(), gamma.definition(), hasAssociatedTrack); + //_______________________________________________ // Lambda extra properties auto posTrackLambda = lambda.template posTrackExtra_as(); @@ -2624,31 +2619,31 @@ struct sigma0builder { // Auxiliary vectors to store best candidates std::vector bestGammasArray; - std::vector bestLambdasArray; + std::vector bestLambdasArray; std::vector bestKShortsArray; // Custom grouping - std::vector> v0grouped(collisions.size()); - std::vector> emclustersgrouped(collisions.size()); + std::vector> v0grouped(collisions.size()); + std::vector> emclustersgrouped(collisions.size()); std::vector emcaltracksgrouped; // Grouping step: for (const auto& v0 : fullV0s) { v0grouped[v0.straCollisionId()].push_back(v0.globalIndex()); } - - if constexpr (soa::is_table) { + + if constexpr (soa::is_table) { emcaltracksgrouped.resize(fullEMCalClusters.size(), false); for (const auto& cluster : fullEMCalClusters) { emclustersgrouped[cluster.collisionId()].push_back(cluster.globalIndex()); } - // Mapping emccluster to matched tracks + // Mapping emccluster to matched tracks for (const auto& calotrack : emcaltracks) { emcaltracksgrouped[calotrack.emcalclusterId()] = true; } } - + //_______________________________________________ // Collisions loop for (const auto& coll : collisions) { @@ -2657,15 +2652,15 @@ struct sigma0builder { if (!IsEventAccepted(coll, true)) continue; } - + if constexpr (soa::is_table) { - if (emclustersgrouped[coll.globalIndex()].size()==0 && eventSelections.fSkipEmptyEMCal) + if (emclustersgrouped[coll.globalIndex()].size() == 0 && eventSelections.fSkipEmptyEMCal) continue; - } + } // Clear vectors bestGammasArray.clear(); - bestLambdasArray.clear(); + bestLambdasArray.clear(); bestKShortsArray.clear(); float centrality = doPPAnalysis ? coll.centFT0M() : coll.centFT0C(); @@ -2678,55 +2673,56 @@ struct sigma0builder { if (fFillNoSelV0Histos) fillV0Histos<0>(v0, coll); // Filling "all V0s" histograms - - if (fUsePCMPhoton){ + + if (fUsePCMPhoton) { if (processPhotonCandidate(v0)) { // selecting photons if (fFillSelPhotonHistos) - fillV0Histos<1>(v0, coll); // QA histos + fillV0Histos<1>(v0, coll); // QA histos bestGammasArray.push_back(v0.globalIndex()); // Save indices of best gamma candidates } } - + if (processLambdaCandidate(v0, coll)) { // selecting lambdas if (fFillSelLambdaHistos) - fillV0Histos<2>(v0, coll); // QA histos + fillV0Histos<2>(v0, coll); // QA histos bestLambdasArray.push_back(v0.globalIndex()); // Save indices of best lambda candidates } if (processKShortCandidate(v0, coll)) { // selecting kshorts if (fFillSelKShortHistos) - fillV0Histos<3>(v0, coll); // QA histos + fillV0Histos<3>(v0, coll); // QA histos bestKShortsArray.push_back(v0.globalIndex()); // Save indices of best kshort candidates } } - + // If EMCalClusters is available, we don't use PCM - if constexpr (soa::is_table) { + if constexpr (soa::is_table) { for (size_t i = 0; i < emclustersgrouped[coll.globalIndex()].size(); i++) { auto cluster = fullEMCalClusters.rawIteratorAt(emclustersgrouped[coll.globalIndex()][i]); - fillEMCalHistos<0>(cluster); + fillEMCalHistos<0>(cluster); - if (processEMCalPhotonCandidate(cluster)) { // selecting photons - fillEMCalHistos<1>(cluster); // QA histos + if (processEMCalPhotonCandidate(cluster)) { // selecting photons + fillEMCalHistos<1>(cluster); // QA histos bestGammasArray.push_back(cluster.globalIndex()); // Save indices of best gamma candidates } } } - + //_______________________________________________ // Wrongly collision association study (MC-specific) if constexpr (requires { coll.straMCCollisionId(); }) { if (doAssocStudy) { - if constexpr (!soa::is_table) analyzeV0CollAssoc(coll, fullV0s, bestGammasArray, true); // Photon-analysis + if constexpr (!soa::is_table) + analyzeV0CollAssoc(coll, fullV0s, bestGammasArray, true); // Photon-analysis analyzeV0CollAssoc(coll, fullV0s, bestLambdasArray, false); // Lambda-analysis } } - + //_______________________________________________ // Photon-V0 nested loop int nSigma0Candidates = 0; for (size_t i = 0; i < bestGammasArray.size(); ++i) { - + //_______________________________________________ // Sigma0 loop if (fillSigma0Tables) { @@ -2734,52 +2730,52 @@ struct sigma0builder { auto lambda = fullV0s.rawIteratorAt(bestLambdasArray[j]); // Building sigma0 candidate & filling tables - if constexpr (soa::is_table) { // using EMCal photons + if constexpr (soa::is_table) { // using EMCal photons auto gamma1 = fullEMCalClusters.rawIteratorAt(bestGammasArray[i]); - if (!buildEMCalSigma0(lambda, gamma1, coll, mcparticles, emcaltracksgrouped)) continue; - } - else { // using PCM photons + if (!buildEMCalSigma0(lambda, gamma1, coll, mcparticles, emcaltracksgrouped)) + continue; + } else { // using PCM photons auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); - if (!buildPCMSigma0(lambda, gamma1, coll, mcparticles)) continue; - } - + if (!buildPCMSigma0(lambda, gamma1, coll, mcparticles)) + continue; + } + nSigma0Candidates++; } } //_______________________________________________ - // KStar loop - if constexpr (!soa::is_table){ // Don't use EMCal clusters here + // KStar loop + if constexpr (!soa::is_table) { // Don't use EMCal clusters here if (fillKStarTables) { auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); for (size_t j = 0; j < bestKShortsArray.size(); ++j) { auto kshort = fullV0s.rawIteratorAt(bestKShortsArray[j]); - // Building kstar candidate & filling tables + // Building kstar candidate & filling tables if (!buildKStar(kshort, gamma1, coll, mcparticles)) continue; - } } } - + //_______________________________________________ // pi0 loop - if constexpr (!soa::is_table){ // Don't use EMCal clusters here + if constexpr (!soa::is_table) { // Don't use EMCal clusters here if (fillPi0Tables) { auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); for (size_t j = i + 1; j < bestGammasArray.size(); ++j) { auto gamma2 = fullV0s.rawIteratorAt(bestGammasArray[j]); - // Building pi0 candidate & filling tables + // Building pi0 candidate & filling tables if (!buildPi0(gamma1, gamma2, coll, mcparticles)) continue; } } - } + } } - LOGF(info, "N. photons: %i, N. lambdas: %i, expected pairs: %i, got: %i", bestGammasArray.size(), bestLambdasArray.size(), bestGammasArray.size()*bestLambdasArray.size(), nSigma0Candidates); + LOGF(info, "N. photons: %i, N. lambdas: %i, expected pairs: %i, got: %i", bestGammasArray.size(), bestLambdasArray.size(), bestGammasArray.size() * bestLambdasArray.size(), nSigma0Candidates); } } @@ -2907,10 +2903,10 @@ struct sigma0builder { template void runPCMVsEMCalQA(TCollision const& collisions, TV0s const& fullV0s, TEMCal const& fullEMCalClusters, TEMCalTracks const& emcaltracks, TMCParticles const& mcparticles) { - + // Custom grouping - std::vector> v0grouped(collisions.size()); - std::vector> emclustersgrouped(collisions.size()); + std::vector> v0grouped(collisions.size()); + std::vector> emclustersgrouped(collisions.size()); std::vector emcaltracksgrouped(collisions.size(), false); // Grouping step: @@ -2922,7 +2918,7 @@ struct sigma0builder { emclustersgrouped[cluster.collisionId()].push_back(cluster.globalIndex()); } - // Mapping emccluster to matched tracks + // Mapping emccluster to matched tracks for (const auto& calotrack : emcaltracks) { emcaltracksgrouped[calotrack.emcalclusterId()] = true; } @@ -2954,27 +2950,26 @@ struct sigma0builder { auto v0MC = v0.template v0MCCore_as>(); auto mcv0Photon = mcparticles.rawIteratorAt(v0MC.particleIdMC()); - if (mcv0Photon.pdgCode() != PDG_t::kGamma || !mcv0Photon.isPhysicalPrimary() || TMath::Abs(mcv0Photon.y()) > 0.5) + if (mcv0Photon.pdgCode() != PDG_t::kGamma || !mcv0Photon.isPhysicalPrimary() || TMath::Abs(mcv0Photon.y()) > 0.5) continue; - + // MC pT histo histos.fill(HIST("PhotonMCQA/hPCMPhotonMCpT"), mcv0Photon.pt()); - + // pT resolution - histos.fill(HIST("PhotonMCQA/h2dPCMPhotonMCpTResolution"), mcv0Photon.pt(), 1 - v0.pt()/mcv0Photon.pt()); + histos.fill(HIST("PhotonMCQA/h2dPCMPhotonMCpTResolution"), mcv0Photon.pt(), 1 - v0.pt() / mcv0Photon.pt()); // photon from sigma0/asigma0 - auto const& v0photonMothers = mcv0Photon.template mothers_as(); + auto const& v0photonMothers = mcv0Photon.template mothers_as(); if (!v0photonMothers.empty()) { // Assumption: first mother is the physical one auto const& v0photonMother = v0photonMothers.front(); if (TMath::Abs(v0photonMother.pdgCode()) == PDG_t::kSigma0) { // Sigma0 or ASigma0 // For efficiency histos.fill(HIST("PhotonMCQA/hPCMSigma0PhotonMCpT"), mcv0Photon.pt()); - + // pT resolution - histos.fill(HIST("PhotonMCQA/h2dPCMSigma0PhotonMCpTResolution"), mcv0Photon.pt(), 1 - v0.pt()/mcv0Photon.pt()); - + histos.fill(HIST("PhotonMCQA/h2dPCMSigma0PhotonMCpTResolution"), mcv0Photon.pt(), 1 - v0.pt() / mcv0Photon.pt()); } } } @@ -2991,10 +2986,10 @@ struct sigma0builder { // --- Reco kinematics (assume photon mass = 0) // ============================================================ - float E = cluster.energy(); - float eta = cluster.eta(); + float E = cluster.energy(); + float eta = cluster.eta(); - float pt = E / std::cosh(eta); + float pt = E / std::cosh(eta); // ============================================================ // --- MC matching @@ -3017,8 +3012,7 @@ struct sigma0builder { // Case 1: particle itself is photon if (mcPart.pdgCode() == PDG_t::kGamma) { photonId = mcPart.globalIndex(); - } - else { + } else { // Case 2: climb ancestry to find photon int candidateId = aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain( @@ -3046,25 +3040,25 @@ struct sigma0builder { // Select TRUE + PRIMARY photons // ============================================================ - if (mcPhoton.pdgCode() != PDG_t::kGamma || !mcPhoton.isPhysicalPrimary() || TMath::Abs(mcPhoton.y()) > 0.5) + if (mcPhoton.pdgCode() != PDG_t::kGamma || !mcPhoton.isPhysicalPrimary() || TMath::Abs(mcPhoton.y()) > 0.5) continue; - + // ============================================================ // Fill QA histos // ============================================================ // MC pT histo histos.fill(HIST("PhotonMCQA/hEMCalPhotonMCpT"), mcPhoton.pt()); - + // pT resolution - histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCpTResolution"), mcPhoton.pt(), 1 - pt/mcPhoton.pt()); + histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCpTResolution"), mcPhoton.pt(), 1 - pt / mcPhoton.pt()); // Energy resolution - histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCEnergyResolution"), mcPhoton.e(), 1 - cluster.energy()/mcPhoton.e()); + histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCEnergyResolution"), mcPhoton.e(), 1 - cluster.energy() / mcPhoton.e()); // Eta resolution histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCEtaResolution"), mcPhoton.eta(), cluster.eta() - mcPhoton.eta()); - + // Phi resolution histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCPhiResolution"), mcPhoton.phi(), cluster.phi() - mcPhoton.phi()); @@ -3072,50 +3066,50 @@ struct sigma0builder { histos.fill(HIST("PhotonMCQA/h2dEMCalPhotonMCFractionEnergy"), mcPhoton.pt(), cluster.amplitudeA()[i]); // photon from sigma0/asigma0 - auto const& photonMothers = mcPhoton.template mothers_as(); + auto const& photonMothers = mcPhoton.template mothers_as(); if (!photonMothers.empty()) { // Assumption: first mother is the physical one auto const& photonMother = photonMothers.front(); if (TMath::Abs(photonMother.pdgCode()) == PDG_t::kSigma0) { // Sigma0 or ASigma0 // For efficiency histos.fill(HIST("PhotonMCQA/hEMCalSigma0PhotonMCpT"), mcPhoton.pt()); - + // pT resolution - histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCpTResolution"), mcPhoton.pt(), 1 - pt/mcPhoton.pt()); + histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCpTResolution"), mcPhoton.pt(), 1 - pt / mcPhoton.pt()); // Energy resolution - histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCEnergyResolution"), mcPhoton.e(), 1 - cluster.energy()/mcPhoton.e()); + histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCEnergyResolution"), mcPhoton.e(), 1 - cluster.energy() / mcPhoton.e()); // Eta resolution histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCEtaResolution"), mcPhoton.eta(), cluster.eta() - mcPhoton.eta()); - + // Phi resolution histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCPhiResolution"), mcPhoton.phi(), cluster.phi() - mcPhoton.phi()); // Fraction of energy Vs MC pT histos.fill(HIST("PhotonMCQA/h2dEMCalSigma0PhotonMCFractionEnergy"), mcPhoton.pt(), cluster.amplitudeA()[i]); } - } + } } } - + } // end of collisions loop - // Process MC generated photons + // Process MC generated photons for (const auto& mcpart : mcparticles) { - if (mcpart.pdgCode() != PDG_t::kGamma || !mcpart.isPhysicalPrimary() || TMath::Abs(mcpart.y()) > 0.5) + if (mcpart.pdgCode() != PDG_t::kGamma || !mcpart.isPhysicalPrimary() || TMath::Abs(mcpart.y()) > 0.5) continue; - + histos.fill(HIST("PhotonMCQA/hGenPhoton"), mcpart.pt()); // photon from sigma0/asigma0 - auto const& photonMothers = mcpart.template mothers_as(); + auto const& photonMothers = mcpart.template mothers_as(); if (!photonMothers.empty()) { // Assumption: first mother is the physical one auto const& photonMother = photonMothers.front(); if (TMath::Abs(photonMother.pdgCode()) == PDG_t::kSigma0) { // Sigma0 or ASigma0 histos.fill(HIST("PhotonMCQA/hGenSigma0Photon"), mcpart.pt()); - } + } } } } diff --git a/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx b/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx index 89cb5982da1..b41eae69fec 100644 --- a/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx @@ -65,7 +65,6 @@ using MCSigma0s = soa::Join; using MCSigma0sWithEMCal = soa::Join; - static const std::vector PhotonSels = {"NoSel", "V0Type", "DCADauToPV", "DCADau", "DauTPCCR", "TPCNSigmaEl", "V0pT", "Y", "V0Radius", "RZCut", "Armenteros", "CosPA", @@ -223,17 +222,17 @@ struct sigmaanalysis { struct : ConfigurableGroup { std::string prefix = "EMCalPhotonSelections"; // JSON group name Configurable definition{"definition", 13, "Cluster definitions to be accepted (e.g. 13 for kV3MostSplitLowSeed)"}; - Configurable MinCells{"MinCells", 1, "Min number of cells in cluster"}; - Configurable MinEnergy{"MinEnergy", 0.0, "Minimum energy of selected clusters (GeV)"}; - Configurable MaxEnergy{"MaxEnergy", 5, "Max energy of selected clusters (GeV)"}; - Configurable MaxEta{"MaxEta", 1.0, "Max absolute cluster Eta"}; - Configurable MinTime{"MinTime", -50, "Minimum time of selected clusters (ns)"}; - Configurable MaxTime{"MaxTime", 50, "Max time of selected clusters (ns)"}; - Configurable RemoveExotic{"RemoveExotic", false, "Flag to enable the removal of exotic clusters"}; - Configurable MinM02{"MinM02", -1., "Minimum shower shape long axis"}; - Configurable MaxM02{"MaxM02", 5., "Max shower shape long axis"}; - Configurable RemoveMatchedTrack{"RemoveMatchedTrack", true, "Flag to enable the removal of clusters matched to tracks"}; - + Configurable MinCells{"MinCells", 1, "Min number of cells in cluster"}; + Configurable MinEnergy{"MinEnergy", 0.0, "Minimum energy of selected clusters (GeV)"}; + Configurable MaxEnergy{"MaxEnergy", 5, "Max energy of selected clusters (GeV)"}; + Configurable MaxEta{"MaxEta", 1.0, "Max absolute cluster Eta"}; + Configurable MinTime{"MinTime", -50, "Minimum time of selected clusters (ns)"}; + Configurable MaxTime{"MaxTime", 50, "Max time of selected clusters (ns)"}; + Configurable RemoveExotic{"RemoveExotic", false, "Flag to enable the removal of exotic clusters"}; + Configurable MinM02{"MinM02", -1., "Minimum shower shape long axis"}; + Configurable MaxM02{"MaxM02", 5., "Max shower shape long axis"}; + Configurable RemoveMatchedTrack{"RemoveMatchedTrack", true, "Flag to enable the removal of clusters matched to tracks"}; + } EMCalPhotonSelections; struct : ConfigurableGroup { @@ -309,7 +308,7 @@ struct sigmaanalysis { void init(InitContext const&) { LOGF(info, "Initializing now: cross-checking correctness..."); - if ((doprocessRealData + doprocessRealDataWithEMCal + doprocessMonteCarlo + doprocessMonteCarloWithEMCal+ doprocessPi0RealData + doprocessPi0MonteCarlo > 1) || + if ((doprocessRealData + doprocessRealDataWithEMCal + doprocessMonteCarlo + doprocessMonteCarloWithEMCal + doprocessPi0RealData + doprocessPi0MonteCarlo > 1) || (doprocessGeneratedRun3 + doprocessPi0GeneratedRun3 > 1)) { LOGF(fatal, "You have enabled more than one process function. Please check your configuration! Aborting now."); } @@ -389,7 +388,6 @@ struct sigmaanalysis { histos.add(histodir + "/Photon/hMass", "hMass", kTH1D, {axisPhotonMass}); histos.add(histodir + "/Photon/h3dPhotonYSigma0Mass", "h3dPhotonYSigma0Mass", kTH3D, {axisRapidity, axisPt, axisSigmaMass}); histos.add(histodir + "/Photon/h3dPhotonRadiusSigma0Mass", "h3dPhotonRadiusSigma0Mass", kTH3D, {axisV0Radius, axisPt, axisSigmaMass}); - } if (doprocessRealDataWithEMCal || doprocessMonteCarloWithEMCal) { @@ -401,7 +399,7 @@ struct sigmaanalysis { histos.add(histodir + "/EMCalPhotonSel/h2dEtaVsPhi", "h2dEtaVsPhi", kTH2D, {axisRapidity, axisPhi}); histos.add(histodir + "/EMCalPhotonSel/hTime", "hTime", kTH1D, {axisClrTime}); histos.add(histodir + "/EMCalPhotonSel/hExotic", "hExotic", kTH1D, {{2, -0.5f, 1.5f}}); - histos.add(histodir + "/EMCalPhotonSel/hShape", "hShape", kTH1D, {axisClrShape}); + histos.add(histodir + "/EMCalPhotonSel/hShape", "hShape", kTH1D, {axisClrShape}); } histos.add(histodir + "/Lambda/hTrackCode", "hTrackCode", kTH1D, {{11, 0.5f, 11.5f}}); @@ -442,7 +440,7 @@ struct sigmaanalysis { histos.add(histodir + "/Sigma0/hRadius", "hRadius", kTH1D, {axisV0PairRadius}); histos.add(histodir + "/Sigma0/h2dRadiusVspT", "h2dRadiusVspT", kTH2D, {axisV0PairRadius, axisPt}); histos.add(histodir + "/Sigma0/hDCAPairDau", "hDCAPairDau", kTH1D, {axisDCAdau}); - histos.add(histodir + "/Sigma0/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisSigmaMass}); + histos.add(histodir + "/Sigma0/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisSigmaMass}); histos.add(histodir + "/Sigma0/h3dOPAngleVsMass", "h3dOPAngleVsMass", kTH3D, {{140, 0.0f, +7.0f}, axisPt, axisSigmaMass}); histos.add(histodir + "/ASigma0/hMass", "hMass", kTH1D, {axisSigmaMass}); @@ -451,7 +449,7 @@ struct sigmaanalysis { histos.add(histodir + "/ASigma0/hRadius", "hRadius", kTH1D, {axisV0PairRadius}); histos.add(histodir + "/ASigma0/h2dRadiusVspT", "h2dRadiusVspT", kTH2D, {axisV0PairRadius, axisPt}); histos.add(histodir + "/ASigma0/hDCAPairDau", "hDCAPairDau", kTH1D, {axisDCAdau}); - histos.add(histodir + "/ASigma0/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisSigmaMass}); + histos.add(histodir + "/ASigma0/h3dMass", "h3dMass", kTH3D, {axisCentrality, axisPt, axisSigmaMass}); histos.add(histodir + "/ASigma0/h3dOPAngleVsMass", "h3dOPAngleVsMass", kTH3D, {{140, 0.0f, +7.0f}, axisPt, axisSigmaMass}); // Process MC @@ -460,14 +458,14 @@ struct sigmaanalysis { if (fillSelhistos) { histos.add(histodir + "/MC/Photon/hV0ToCollAssoc", "hV0ToCollAssoc", kTH1D, {{2, 0.0f, 2.0f}}); histos.add(histodir + "/MC/Photon/hPt", "hPt", kTH1D, {axisPt}); - histos.add(histodir + "/MC/Photon/hMCPt", "hMCPt", kTH1D, {axisPt}); + histos.add(histodir + "/MC/Photon/hMCPt", "hMCPt", kTH1D, {axisPt}); - if (doprocessMonteCarlo){ + if (doprocessMonteCarlo) { histos.add(histodir + "/MC/Photon/h2dPAVsPt", "h2dPAVsPt", kTH2D, {axisPA, axisPt}); histos.add(histodir + "/MC/Photon/hPt_BadCollAssig", "hPt_BadCollAssig", kTH1D, {axisPt}); histos.add(histodir + "/MC/Photon/h2dPAVsPt_BadCollAssig", "h2dPAVsPt_BadCollAssig", kTH2D, {axisPA, axisPt}); } - + histos.add(histodir + "/MC/Lambda/hV0ToCollAssoc", "hV0ToCollAssoc", kTH1D, {{2, 0.0f, 2.0f}}); histos.add(histodir + "/MC/Lambda/hPt", "hPt", kTH1D, {axisPt}); histos.add(histodir + "/MC/Lambda/hMCPt", "hMCPt", kTH1D, {axisPt}); @@ -480,7 +478,7 @@ struct sigmaanalysis { histos.add(histodir + "/MC/ALambda/h3dTPCvsTOFNSigma_Pr", "h3dTPCvsTOFNSigma_Pr", kTH3D, {axisTPCNSigma, axisTOFNSigma, axisPt}); histos.add(histodir + "/MC/ALambda/h3dTPCvsTOFNSigma_Pi", "h3dTPCvsTOFNSigma_Pi", kTH3D, {axisTPCNSigma, axisTOFNSigma, axisPt}); } - + histos.add(histodir + "/MC/Sigma0/hPt", "hPt", kTH1D, {axisPt}); histos.add(histodir + "/MC/Sigma0/hMCPt", "hMCPt", kTH1D, {axisPt}); histos.add(histodir + "/MC/Sigma0/hMass", "hMass", kTH1D, {axisSigmaMass}); @@ -797,7 +795,7 @@ struct sigmaanalysis { groupedCollisions[coll.straMCCollisionId()].push_back(coll.globalIndex()); } - for (auto const& mcCollision : mcCollisions) { + for (auto const& mcCollision : mcCollisions) { int biggestNContribs = -1; int bestCollisionIndex = -1; for (size_t i = 0; i < groupedCollisions[mcCollision.globalIndex()].size(); i++) { @@ -852,13 +850,13 @@ struct sigmaanalysis { } histos.fill(HIST("Gen/hGenEvents"), mcCollision.multMCNParticlesEta05(), 0 /* all gen. events*/); - + // Check if there is at least one of the reconstructed collisions associated to this MC collision // If so, we consider it bool atLeastOne = false; int biggestNContribs = -1; float centrality = 100.5f; - int nCollisions = 0; + int nCollisions = 0; for (size_t i = 0; i < groupedCollisions[mcCollision.globalIndex()].size(); i++) { auto collision = collisions.rawIteratorAt(groupedCollisions[mcCollision.globalIndex()][i]); @@ -981,7 +979,7 @@ struct sigmaanalysis { int TrkCode = 10; // 1: TPC-only, 2: TPC+Something, 3: ITS-Only, 4: ITS+TPC + Something, 10: anything else if (isGamma) { - if constexpr (requires { sigma.photonV0Type();}) { + if constexpr (requires { sigma.photonV0Type(); }) { if (sigma.photonPosTrackCode() == 1 && sigma.photonNegTrackCode() == 1) TrkCode = 1; if ((sigma.photonPosTrackCode() != 1 && sigma.photonNegTrackCode() == 1) || (sigma.photonPosTrackCode() == 1 && sigma.photonNegTrackCode() != 1)) @@ -1012,9 +1010,9 @@ struct sigmaanalysis { //_______________________________________ // Gamma MC association if (sigma.photonPDGCode() == 22) { - if (sigma.photonmcpt() > 0) { - histos.fill(HIST("BeforeSel/MC/Reso/h2dGammaPtResolution"), sigma.photonmcpt(), sigma.photonPt() - sigma.photonmcpt()); // pT resolution - histos.fill(HIST("BeforeSel/MC/Reso/h2dGammaInvPtResolution"), 1.f / sigma.photonmcpt(), 1.f / sigma.photonPt() - 1.f / sigma.photonmcpt()); // pT resolution + if (sigma.photonmcpt() > 0) { + histos.fill(HIST("BeforeSel/MC/Reso/h2dGammaPtResolution"), sigma.photonmcpt(), sigma.photonPt() - sigma.photonmcpt()); // pT resolution + histos.fill(HIST("BeforeSel/MC/Reso/h2dGammaInvPtResolution"), 1.f / sigma.photonmcpt(), 1.f / sigma.photonPt() - 1.f / sigma.photonmcpt()); // pT resolution } } @@ -1042,10 +1040,10 @@ struct sigmaanalysis { // Sigma and AntiSigma MC association if (sigma.isSigma0()) { histos.fill(HIST("BeforeSel/MC/Reso/h2dSigma0RadiusResolution"), sigma.mcpt(), sigma.radius() - sigma.mcradius()); // pT resolution - if (sigma.mcpt() > 0){ - histos.fill(HIST("BeforeSel/MC/Reso/h2dSigma0PtResolution"), sigma.mcpt(), sigma.pt() - sigma.mcpt()); // pT resolution + if (sigma.mcpt() > 0) { + histos.fill(HIST("BeforeSel/MC/Reso/h2dSigma0PtResolution"), sigma.mcpt(), sigma.pt() - sigma.mcpt()); // pT resolution histos.fill(HIST("BeforeSel/MC/Reso/h2dSigma0InvPtResolution"), 1.f / sigma.mcpt(), 1.f / sigma.pt() - 1.f / sigma.mcpt()); // pT resolution - } + } } if (sigma.isAntiSigma0()) { histos.fill(HIST("BeforeSel/MC/Reso/h2dASigma0RadiusResolution"), sigma.mcpt(), sigma.radius() - sigma.mcradius()); // pT resolution @@ -1101,15 +1099,15 @@ struct sigmaanalysis { // Check whether it is before or after selections static constexpr std::string_view MainDir[] = {"BeforeSel", "AfterSel"}; - + float centrality = getCentralityRun3(collision); if (fillSelhistos) { - + // Get V0trackCode int LambdaTrkCode = retrieveV0TrackCode(sigma); - - if constexpr (requires { sigma.photonV0Type();}) { // Processing PCM photon + + if constexpr (requires { sigma.photonV0Type(); }) { // Processing PCM photon // Base properties int GammaTrkCode = retrieveV0TrackCode(sigma); @@ -1144,16 +1142,15 @@ struct sigmaanalysis { histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenteros"), sigma.photonAlpha(), sigma.photonQt()); histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h3dPhotonYSigma0Mass"), sigma.photonY(), sigma.pt(), sigma.sigma0Mass()); histos.fill(HIST(MainDir[mode]) + HIST("/Photon/h3dPhotonRadiusSigma0Mass"), sigma.photonRadius(), sigma.pt(), sigma.sigma0Mass()); - } - if constexpr (requires { sigma.photonDefinition();}) { // Processing EMCal photon - + if constexpr (requires { sigma.photonDefinition(); }) { // Processing EMCal photon + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hpT"), sigma.photonPt()); histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hDefinition"), sigma.photonDefinition()); histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hNCells"), sigma.photonNCells()); histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hEnergy"), sigma.photonEnergy()); - histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/h2dEtaVsPhi"), sigma.photonEMCEta(), sigma.photonEMCPhi()); + histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/h2dEtaVsPhi"), sigma.photonEMCEta(), sigma.photonEMCPhi()); histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hTime"), sigma.photonTime()); histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hExotic"), sigma.photonIsExotic()); histos.fill(HIST(MainDir[mode]) + HIST("/EMCalPhotonSel/hShape"), sigma.photonM02()); @@ -1176,11 +1173,11 @@ struct sigmaanalysis { histos.fill(HIST(MainDir[mode]) + HIST("/Lambda/hNegChi2PerNc"), sigma.lambdaNegChi2PerNcl()); histos.fill(HIST(MainDir[mode]) + HIST("/Lambda/hLifeTime"), sigma.lambdaLifeTime()); - histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenteros"), sigma.lambdaAlpha(), sigma.lambdaQt()); + histos.fill(HIST(MainDir[mode]) + HIST("/h2dArmenteros"), sigma.lambdaAlpha(), sigma.lambdaQt()); } //_______________________________________ - // Sigmas and Lambdas + // Sigmas and Lambdas if (sigma.lambdaAlpha() > 0) { if (fillSelhistos) { histos.fill(HIST(MainDir[mode]) + HIST("/Lambda/h2dTPCvsTOFNSigma_LambdaPr"), sigma.lambdaPosPrTPCNSigma(), sigma.lambdaPrTOFNSigma()); @@ -1198,7 +1195,7 @@ struct sigmaanalysis { histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/hRadius"), sigma.radius()); histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h2dRadiusVspT"), sigma.radius(), sigma.pt()); histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/hDCAPairDau"), sigma.dcadaughters()); - histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h3dMass"), centrality, sigma.pt(), sigma.sigma0Mass()); + histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h3dMass"), centrality, sigma.pt(), sigma.sigma0Mass()); histos.fill(HIST(MainDir[mode]) + HIST("/Sigma0/h3dOPAngleVsMass"), sigma.opAngle(), sigma.pt(), sigma.sigma0Mass()); } else { if (fillSelhistos) { @@ -1217,7 +1214,7 @@ struct sigmaanalysis { histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/hRadius"), sigma.radius()); histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h2dRadiusVspT"), sigma.radius(), sigma.pt()); histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/hDCAPairDau"), sigma.dcadaughters()); - histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h3dMass"), centrality, sigma.pt(), sigma.sigma0Mass()); + histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h3dMass"), centrality, sigma.pt(), sigma.sigma0Mass()); histos.fill(HIST(MainDir[mode]) + HIST("/ASigma0/h3dOPAngleVsMass"), sigma.opAngle(), sigma.pt(), sigma.sigma0Mass()); } @@ -1233,9 +1230,9 @@ struct sigmaanalysis { histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hV0ToCollAssoc"), sigma.photonIsCorrectlyAssoc()); histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hPt"), sigma.photonPt()); - histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hMCPt"), sigma.photonmcpt()); - - if constexpr (requires { sigma.photonV0Type();}) { // Processing PCM photon + histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/hMCPt"), sigma.photonmcpt()); + + if constexpr (requires { sigma.photonV0Type(); }) { // Processing PCM photon histos.fill(HIST(MainDir[mode]) + HIST("/MC/Photon/h2dPAVsPt"), TMath::ACos(sigma.photonCosPA()), sigma.photonmcpt()); if (!sigma.photonIsCorrectlyAssoc()) { @@ -1267,7 +1264,7 @@ struct sigmaanalysis { } //_______________________________________ // Sigma0 MC association - if (sigma.isSigma0()) { + if (sigma.isSigma0()) { histos.fill(HIST(MainDir[mode]) + HIST("/MC/Sigma0/hPt"), sigma.pt()); histos.fill(HIST(MainDir[mode]) + HIST("/MC/Sigma0/hMCPt"), sigma.mcpt()); @@ -1455,43 +1452,43 @@ struct sigmaanalysis { // Apply specific selections for photons template - bool selectEMCalPhoton(TClusterObject const& cand) + bool selectEMCalPhoton(TClusterObject const& cand) { - // Clusterizer + // Clusterizer if (cand.photonDefinition() != EMCalPhotonSelections.definition && EMCalPhotonSelections.definition > -1) return false; - // Number of Cells + // Number of Cells if (cand.photonNCells() < EMCalPhotonSelections.MinCells) return false; - // Energy + // Energy if (cand.photonEnergy() < EMCalPhotonSelections.MinEnergy || cand.photonEnergy() > EMCalPhotonSelections.MaxEnergy) return false; - - // Eta + + // Eta if (TMath::Abs(cand.photonEMCEta()) > EMCalPhotonSelections.MaxEta) return false; - - // Timing + + // Timing if (cand.photonTime() < EMCalPhotonSelections.MinTime || cand.photonTime() > EMCalPhotonSelections.MaxTime) return false; - // Exotic Clusters - if (cand.photonIsExotic() && EMCalPhotonSelections.RemoveExotic) { + // Exotic Clusters + if (cand.photonIsExotic() && EMCalPhotonSelections.RemoveExotic) { return false; } - // Shower shape long axis - if (cand.photonNCells() > 1) { // Only if we have more than one - if (cand.photonM02() < EMCalPhotonSelections.MinM02 || cand.photonM02() > EMCalPhotonSelections.MaxM02) { + // Shower shape long axis + if (cand.photonNCells() > 1) { // Only if we have more than one + if (cand.photonM02() < EMCalPhotonSelections.MinM02 || cand.photonM02() > EMCalPhotonSelections.MaxM02) { return false; } } // Has matched track? if (cand.photonHasAssocTrk() && EMCalPhotonSelections.RemoveMatchedTrack) - return false; + return false; } // Apply specific selections for lambdas @@ -1603,13 +1600,13 @@ struct sigmaanalysis { template bool processSigma0Candidate(TSigma0Object const& cand) { - // Photon specific selections - if constexpr (requires { cand.photonV0Type();}) { // Processing PCM photon + // Photon specific selections + if constexpr (requires { cand.photonV0Type(); }) { // Processing PCM photon if (!selectPhoton(cand)) return false; } - if constexpr (requires { cand.photonDefinition();}) { // Processing EMCal photon + if constexpr (requires { cand.photonDefinition(); }) { // Processing EMCal photon if (!selectEMCalPhoton(cand)) return false; } @@ -1678,7 +1675,7 @@ struct sigmaanalysis { fillHistos<0>(sigma0, coll); // Select sigma0 candidates - if (!processSigma0Candidate(sigma0)) + if (!processSigma0Candidate(sigma0)) continue; // Fill histos after all selections From 12023c9a3b80aa56819dd0bccb9e1edf8b295a57 Mon Sep 17 00:00:00 2001 From: gianniliveraro Date: Mon, 30 Mar 2026 11:20:25 -0300 Subject: [PATCH 4/4] Silly fix in return of bool function --- PWGLF/Tasks/Strangeness/sigmaanalysis.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx b/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx index b41eae69fec..fd6b48609e1 100644 --- a/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx @@ -1485,10 +1485,11 @@ struct sigmaanalysis { return false; } } - // Has matched track? if (cand.photonHasAssocTrk() && EMCalPhotonSelections.RemoveMatchedTrack) return false; + + return true; } // Apply specific selections for lambdas