diff --git a/PWGLF/DataModel/LFSigmaTables.h b/PWGLF/DataModel/LFSigmaTables.h index 9d728bab596..6c8c9c3a232 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, @@ -502,6 +500,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(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); +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); +DECLARE_SOA_COLUMN(PhotonHasAssocTrk, photonHasAssocTrk, bool); + +} // namespace sigma0EMPhoton + +DECLARE_SOA_TABLE(Sigma0EMPhotons, "AOD", "SIGMA0EMPHOTON", + sigma0EMPhoton::PhotonID, + sigma0EMPhoton::PhotonEnergy, + sigma0EMPhoton::PhotonEMCEta, + sigma0EMPhoton::PhotonEMCPhi, + sigma0EMPhoton::PhotonM02, + sigma0EMPhoton::PhotonM20, + sigma0EMPhoton::PhotonNCells, + sigma0EMPhoton::PhotonTime, + sigma0EMPhoton::PhotonIsExotic, + sigma0EMPhoton::PhotonDistToBad, + sigma0EMPhoton::PhotonNLM, + sigma0EMPhoton::PhotonDefinition, + sigma0EMPhoton::PhotonHasAssocTrk); + // For Lambda extra info namespace sigma0LambdaExtra { @@ -567,6 +601,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); @@ -576,6 +611,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); @@ -585,6 +623,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); @@ -691,15 +731,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, @@ -728,11 +770,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 { @@ -923,12 +960,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/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..bad2b681d13 100644 --- a/PWGLF/TableProducer/Strangeness/sigma0builder.cxx +++ b/PWGLF/TableProducer/Strangeness/sigma0builder.cxx @@ -20,6 +20,8 @@ // 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" @@ -66,8 +68,10 @@ using V0DerivedMCDatas = soa::Join; using V0TOFDerivedMCDatas = soa::Join; +using EMCalMCClusters = soa::Join; + static const std::vector DirList = {"V0BeforeSel", "PhotonSel", "LambdaSel", "KShortSel"}; -; +static const std::vector DirList2 = {"EMCalPhotonBeforeSel", "EMCalPhotonSel"}; struct sigma0builder { Service ccdb; @@ -75,7 +79,6 @@ struct sigma0builder { //___________________________________________________ // KStar Specific - Produces kstarcores; // kstar candidates info for analysis Produces kshortExtras; // lambdas from sigma0 candidates info Produces kstarPhotonExtras; // photons from kstar candidates info @@ -86,15 +89,14 @@ 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 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 sigma0Gens; // Generated sigma0s + Produces sigma0GenCollRefs; // references collisions from sigma0Gens //__________________________________________________ // Pi0 specific @@ -161,8 +163,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 +185,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 +214,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 +239,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)"}; @@ -315,16 +338,22 @@ struct sigma0builder { ConfigurableAxis axisIRBinning{"axisIRBinning", {151, -10, 1500}, "Binning for the interaction rate (kHz)"}; ConfigurableAxis axisLifetime{"axisLifetime", {200, 0, 50}, "Lifetime"}; - // For manual sliceBy (necessary to calculate the correction factors) - PresliceUnsorted> perMcCollision = aod::v0data::straMCCollisionId; + // 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"}; void init(InitContext const&) { LOGF(info, "Initializing now: cross-checking correctness..."); if (doprocessRealData + doprocessRealDataWithTOF + + doprocessRealDataWithEMCal + doprocessMonteCarlo + doprocessMonteCarloWithTOF + + doprocessMonteCarloWithEMCal + doprocessV0QA + doprocessV0MCQA > 1) { @@ -373,9 +402,11 @@ struct sigma0builder { } } + bool fUsePCMPhoton = !doprocessRealDataWithEMCal && !doprocessMonteCarloWithEMCal && !doprocessPCMVsEMCalQA; + for (const auto& histodir : DirList) { if ((histodir == "V0BeforeSel" && !fFillNoSelV0Histos) || - (histodir == "PhotonSel" && !fFillSelPhotonHistos) || + (histodir == "PhotonSel" && !fFillSelPhotonHistos && !fUsePCMPhoton) || (histodir == "LambdaSel" && !fFillSelLambdaHistos) || (histodir == "KShortSel" && !fFillSelKShortHistos)) { continue; @@ -424,20 +455,43 @@ 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 || 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"); + 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"); + + 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"); + + } 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}); histos.get(HIST("LambdaSel/hSelectionStatistics"))->GetXaxis()->SetBinLabel(1, "No Sel"); @@ -471,13 +525,14 @@ 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"); 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}); @@ -501,7 +556,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}}); @@ -529,6 +584,30 @@ 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 @@ -705,6 +784,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; @@ -720,6 +803,7 @@ struct sigma0builder { float V02MCpy = -999.f; float V02MCpz = -999.f; float V0PairMCRadius = -999.f; + float EMCalClusterAmplitude = -999.f; }; // ______________________________________________________ @@ -831,6 +915,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(); @@ -933,6 +1021,169 @@ struct sigma0builder { return MCinfo; } + template + V0PairMCInfo getClusterV0PairMCInfo(TEMCalCls const& cluster, + TV0 const& v0, + TCollision const& collision, + TMCParticles const& mcparticles) + { + // Output container + V0PairMCInfo MCinfo; + + // ============================================================ + // 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 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< + soa::Join>(); + MCinfo.fIsV02CorrectlyAssign = (v0MC.straMCCollisionId() == MCCollision.globalIndex()); + } + + // Retrieve Lambda mothers + auto const& lambdaMothers = mcLambda.template mothers_as(); + if (lambdaMothers.empty()) { + // No ancestry -> cannot match to any parent + return MCinfo; + } + + // 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; + + // Get Sigma0 index from photon mother + auto mothers = mcPhoton.mothersIds(); + if (mothers.empty()) // No mothers? Weird + continue; + + 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 + } + } + + // ============================================================ + // 3) --- SUCCESS: true Lambda–photon pair from same Sigma0 --- + // ============================================================ + + if (matchedPhotonId >= 0 && matchedMotherIndex >= 0) { + + 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(); + } + + // --- Photon MC info + MCinfo.V01MCpx = mcPhoton.px(); + MCinfo.V01MCpy = mcPhoton.py(); + MCinfo.V01MCpz = mcPhoton.pz(); + MCinfo.fIsV01Primary = mcPhoton.isPhysicalPrimary(); + MCinfo.V01PDGCode = mcPhoton.pdgCode(); + + if (!mcPhoton.mothersIds().empty()) { + auto mcMother = mcparticles.iteratorAt(mcPhoton.mothersIds()[0]); + MCinfo.V01PDGCodeMother = mcMother.pdgCode(); + } + + 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; + } + // ______________________________________________________ // Check whether the collision passes our collision selections // Should work with collisions, mccollisions, stracollisions and stramccollisions tables! @@ -1079,12 +1330,21 @@ struct sigma0builder { std::vector getListOfRecoCollIndices(TMCollisions 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) { - auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex()); 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; @@ -1108,6 +1368,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) { @@ -1125,14 +1393,15 @@ 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; @@ -1501,9 +1770,28 @@ 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("/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("/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) 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 +1855,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 +2255,20 @@ 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) { //_______________________________________________ + // 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.negTrackExtraId() == lambda.negTrackExtraId() || + gamma.posTrackExtraId() == lambda.negTrackExtraId() || + gamma.negTrackExtraId() == lambda.posTrackExtraId()) { return false; } @@ -1959,26 +2302,23 @@ 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); - 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); } // Sigma0s -> stracollisions link histos.fill(HIST("SigmaSel/hSigma0DauDeltaIndex"), gamma.globalIndex() - lambda.globalIndex()); sigma0CollRefs(collision.globalIndex()); - sigmaIndices(gamma.globalIndex(), lambda.globalIndex()); //_______________________________________________ // Photon extra properties @@ -2043,6 +2383,111 @@ struct sigma0builder { return true; } + // 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()); + + // 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}); + + // 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; + + 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(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(); }) { + + 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()); + + //_______________________________________________ + // 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); + + //_______________________________________________ + // 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 +2610,40 @@ 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, TEMCalTracks const& emcaltracks, 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 bestKShortsArray; // Custom grouping 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) { + 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; + } + } + //_______________________________________________ // Collisions loop for (const auto& coll : collisions) { @@ -2193,11 +2653,14 @@ struct sigma0builder { continue; } + if constexpr (soa::is_table) { + if (emclustersgrouped[coll.globalIndex()].size() == 0 && eventSelections.fSkipEmptyEMCal) + continue; + } + // Clear vectors bestGammasArray.clear(); bestLambdasArray.clear(); - - bestKStarGammasArray.clear(); bestKShortsArray.clear(); float centrality = doPPAnalysis ? coll.centFT0M() : coll.centFT0C(); @@ -2209,40 +2672,56 @@ 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 - if (processPhotonCandidate(v0)) { // selecting photons - if (fFillSelPhotonHistos) - fillHistos<1>(v0, coll); // QA histos - bestGammasArray.push_back(v0.globalIndex()); // Save indices of best gamma candidates + 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 + int nSigma0Candidates = 0; for (size_t i = 0; i < bestGammasArray.size(); ++i) { - auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); //_______________________________________________ // Sigma0 loop @@ -2251,35 +2730,52 @@ 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, emcaltracksgrouped)) + continue; + } else { // using PCM photons + auto gamma1 = fullV0s.rawIteratorAt(bestGammasArray[i]); + if (!buildPCMSigma0(lambda, gamma1, coll, mcparticles)) + continue; + } + + nSigma0Candidates++; } } //_______________________________________________ // 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; + 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; + } } } } + + LOGF(info, "N. photons: %i, N. lambdas: %i, expected pairs: %i, got: %i", bestGammasArray.size(), bestLambdasArray.size(), bestGammasArray.size() * bestLambdasArray.size(), nSigma0Candidates); } } @@ -2317,7 +2813,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 +2823,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 +2833,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 +2848,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()); @@ -2403,25 +2899,250 @@ 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); + dataProcess(collisions, fullV0s, nullptr, nullptr, nullptr); } void processRealDataWithTOF(soa::Join const& collisions, V0TOFStandardDerivedDatas const& fullV0s, dauTracks const&) { - dataProcess(collisions, fullV0s, nullptr); + dataProcess(collisions, fullV0s, nullptr, nullptr, nullptr); + } + + void processRealDataWithEMCal(soa::Join const& collisions, V0StandardDerivedDatas const& fullV0s, dauTracks const&, aod::EMCALClusters const& fullEMCalClusters, aod::EMCALMatchedTracks const& emcmatchedtracks) + { + 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, 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, 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, aod::EMCALMatchedTracks const& emcmatchedtracks) + { + dataProcess(collisions, fullV0s, fullEMCalMCClusters, emcmatchedtracks, mcParticles); } void processGeneratedRun3(aod::McParticles const& mcParticles) @@ -2445,14 +3166,22 @@ 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", 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); 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..fd6b48609e1 100644 --- a/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx @@ -60,8 +60,10 @@ 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 +90,6 @@ struct sigmaanalysis { ctpRateFetcher rateFetcher; //__________________________________________________ - // For manual sliceBy - // SliceCache cache; - PresliceUnsorted> perMcCollision = aod::v0data::straMCCollisionId; - HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; // Event level @@ -220,6 +218,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 +292,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 +308,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 +355,52 @@ 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}); @@ -400,8 +441,6 @@ struct sigmaanalysis { 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/h3dOPAngleVsMass", "h3dOPAngleVsMass", kTH3D, {{140, 0.0f, +7.0f}, axisPt, axisSigmaMass}); histos.add(histodir + "/ASigma0/hMass", "hMass", kTH1D, {axisSigmaMass}); @@ -411,22 +450,21 @@ struct sigmaanalysis { 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/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}); + + 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}); @@ -441,8 +479,6 @@ struct sigmaanalysis { 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 +501,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 +512,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 +787,21 @@ struct sigmaanalysis { std::vector getListOfRecoCollIndices(TMCollisions 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) { - auto groupedCollisions = collisions.sliceBy(perMcCollision, mcCollision.globalIndex()); 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 +826,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) { @@ -798,14 +851,14 @@ 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) { + 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 +979,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; @@ -956,9 +1011,8 @@ 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 + 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 +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"), 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 @@ -1044,39 +1100,61 @@ 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 +1172,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()); - if (sigma.lambdaAlpha() > 0) { if (fillSelhistos) { histos.fill(HIST(MainDir[mode]) + HIST("/Lambda/h2dTPCvsTOFNSigma_LambdaPr"), sigma.lambdaPosPrTPCNSigma(), sigma.lambdaPrTOFNSigma()); @@ -1119,8 +1196,6 @@ struct sigmaanalysis { 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/h3dOPAngleVsMass"), sigma.opAngle(), sigma.pt(), sigma.sigma0Mass()); } else { if (fillSelhistos) { @@ -1140,8 +1215,6 @@ struct sigmaanalysis { 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/h3dOPAngleVsMass"), sigma.opAngle(), sigma.pt(), sigma.sigma0Mass()); } @@ -1154,17 +1227,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 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()); + 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()); + } } } @@ -1191,9 +1265,6 @@ 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()); - histos.fill(HIST(MainDir[mode]) + HIST("/MC/Sigma0/hPt"), sigma.pt()); histos.fill(HIST(MainDir[mode]) + HIST("/MC/Sigma0/hMCPt"), sigma.mcpt()); @@ -1211,9 +1282,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 +1450,48 @@ 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; + + return true; + } + // Apply specific selections for lambdas template bool selectLambda(TV0Object const& cand) @@ -1492,8 +1602,15 @@ struct sigmaanalysis { bool processSigma0Candidate(TSigma0Object const& cand) { // Photon specific selections - if (!selectPhoton(cand)) - return false; + 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 +1656,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 @@ -1721,11 +1838,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 +1878,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);