diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index a2a3dcc005a..e4e79192278 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -1046,6 +1046,21 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h } } } + if (subGroupStr.Contains("emcal")) { + hm->AddHistogram(histClass, "EMCalEnergy", "EMCal cluster energy", false, 500, 0.0, 50.0, VarManager::kEMCalEnergy); + hm->AddHistogram(histClass, "EMCalCoreEnergy", "EMCal cluster core energy", false, 500, 0.0, 50.0, VarManager::kEMCalCoreEnergy); + hm->AddHistogram(histClass, "EMCalEta_Phi", "EMCal cluster #eta vs #varphi", false, 100, -1.0, 1.0, VarManager::kEMCalEta, 180, 0.0, TMath::TwoPi(), VarManager::kEMCalPhi); + hm->AddHistogram(histClass, "EMCalM02", "EMCal cluster M02", false, 200, 0.0, 2.0, VarManager::kEMCalM02); + hm->AddHistogram(histClass, "EMCalM20", "EMCal cluster M20", false, 200, 0.0, 2.0, VarManager::kEMCalM20); + hm->AddHistogram(histClass, "EMCalNCells", "EMCal cluster no. cells", false, 50, -0.5, 49.5, VarManager::kEMCalNCells); + hm->AddHistogram(histClass, "EMCalTime", "EMCal cluster time", false, 400, -100.0, 100.0, VarManager::kEMCalTime); + hm->AddHistogram(histClass, "EMCalEnergy_Time", "EMCal cluster energy vs time", false, 100, 0.0, 50.0, VarManager::kEMCalEnergy, 200, -100.0, 100.0, VarManager::kEMCalTime); + if (subGroupStr.Contains("emcalmatch")) { + hm->AddHistogram(histClass, "EMCalEoverP", "EMCal E/p", false, 300, 0.0, 3.0, VarManager::kEMCalEoverP); + hm->AddHistogram(histClass, "Pt_EMCalEoverP", "p_{T} vs EMCal E/p", false, 200, 0.0, 20.0, VarManager::kPt, 300, 0.0, 3.0, VarManager::kEMCalEoverP); + hm->AddHistogram(histClass, "EMCalMatchDeltaEta_DeltaPhi", "EMCal matching #Delta#eta vs #Delta#varphi", false, 100, -0.1, 0.1, VarManager::kEMCalMatchDeltaEta, 100, -0.1, 0.1, VarManager::kEMCalMatchDeltaPhi); + } + } if (subGroupStr.Contains("muon")) { if (!subGroupStr.Contains("ambiguity")) { hm->AddHistogram(histClass, "MuonNClusters", "", false, 100, 0.0, 10.0, VarManager::kMuonNClusters); diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 0ae36125960..28b3c909bc5 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -999,6 +999,38 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kTRDsignal] = ""; fgVariableNames[kTOFbeta] = "TOF #beta"; fgVariableUnits[kTOFbeta] = ""; + fgVariableNames[kEMCalEnergy] = "EMCal cluster energy"; + fgVariableUnits[kEMCalEnergy] = "GeV"; + fgVariableNames[kEMCalCoreEnergy] = "EMCal cluster core energy"; + fgVariableUnits[kEMCalCoreEnergy] = "GeV"; + fgVariableNames[kEMCalRawEnergy] = "EMCal cluster raw energy"; + fgVariableUnits[kEMCalRawEnergy] = "GeV"; + fgVariableNames[kEMCalEta] = "EMCal cluster #eta"; + fgVariableUnits[kEMCalEta] = ""; + fgVariableNames[kEMCalPhi] = "EMCal cluster #varphi"; + fgVariableUnits[kEMCalPhi] = "rad."; + fgVariableNames[kEMCalM02] = "EMCal cluster M02"; + fgVariableUnits[kEMCalM02] = ""; + fgVariableNames[kEMCalM20] = "EMCal cluster M20"; + fgVariableUnits[kEMCalM20] = ""; + fgVariableNames[kEMCalNCells] = "EMCal cluster no. cells"; + fgVariableUnits[kEMCalNCells] = ""; + fgVariableNames[kEMCalTime] = "EMCal cluster time"; + fgVariableUnits[kEMCalTime] = "ns"; + fgVariableNames[kEMCalIsExotic] = "EMCal cluster is exotic"; + fgVariableUnits[kEMCalIsExotic] = ""; + fgVariableNames[kEMCalDistanceToBadChannel] = "EMCal cluster distance to bad channel"; + fgVariableUnits[kEMCalDistanceToBadChannel] = ""; + fgVariableNames[kEMCalNLM] = "EMCal cluster no. local maxima"; + fgVariableUnits[kEMCalNLM] = ""; + fgVariableNames[kEMCalDefinition] = "EMCal cluster definition"; + fgVariableUnits[kEMCalDefinition] = ""; + fgVariableNames[kEMCalEoverP] = "E_{EMCal}/p_{track}"; + fgVariableUnits[kEMCalEoverP] = ""; + fgVariableNames[kEMCalMatchDeltaEta] = "#Delta#eta (track-cluster)"; + fgVariableUnits[kEMCalMatchDeltaEta] = ""; + fgVariableNames[kEMCalMatchDeltaPhi] = "#Delta#varphi (track-cluster)"; + fgVariableUnits[kEMCalMatchDeltaPhi] = ""; fgVariableNames[kTrackLength] = "track length"; fgVariableUnits[kTrackLength] = "cm"; fgVariableNames[kTrackDCAxy] = "DCA_{xy}"; @@ -2272,6 +2304,22 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kTRDsignal"] = kTRDsignal; fgVarNamesMap["kTRDPattern"] = kTRDPattern; fgVarNamesMap["kTOFbeta"] = kTOFbeta; + fgVarNamesMap["kEMCalEnergy"] = kEMCalEnergy; + fgVarNamesMap["kEMCalCoreEnergy"] = kEMCalCoreEnergy; + fgVarNamesMap["kEMCalRawEnergy"] = kEMCalRawEnergy; + fgVarNamesMap["kEMCalEta"] = kEMCalEta; + fgVarNamesMap["kEMCalPhi"] = kEMCalPhi; + fgVarNamesMap["kEMCalM02"] = kEMCalM02; + fgVarNamesMap["kEMCalM20"] = kEMCalM20; + fgVarNamesMap["kEMCalNCells"] = kEMCalNCells; + fgVarNamesMap["kEMCalTime"] = kEMCalTime; + fgVarNamesMap["kEMCalIsExotic"] = kEMCalIsExotic; + fgVarNamesMap["kEMCalDistanceToBadChannel"] = kEMCalDistanceToBadChannel; + fgVarNamesMap["kEMCalNLM"] = kEMCalNLM; + fgVarNamesMap["kEMCalDefinition"] = kEMCalDefinition; + fgVarNamesMap["kEMCalEoverP"] = kEMCalEoverP; + fgVarNamesMap["kEMCalMatchDeltaEta"] = kEMCalMatchDeltaEta; + fgVarNamesMap["kEMCalMatchDeltaPhi"] = kEMCalMatchDeltaPhi; fgVarNamesMap["kTrackLength"] = kTrackLength; fgVarNamesMap["kTrackDCAxy"] = kTrackDCAxy; fgVarNamesMap["kTrackDCAxyProng1"] = kTrackDCAxyProng1; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 7079269d946..dfa486bb504 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -152,7 +152,8 @@ class VarManager : public TObject MFTCov = BIT(27), TrackTOFService = BIT(28), ParticleMC = BIT(29), - MuonDca = BIT(30) + MuonDca = BIT(30), + TrackEMCal = BIT(31) // NOTE: last free bit of the uint32 track fill maps; used both for aod::EMCALClusters and ReducedEMCals }; enum PairCandidateType { @@ -193,6 +194,11 @@ class VarManager : public TObject kMuonIsPropagated = 7 // whether the muon was propagated already }; + enum EMCalClusterFilteringBits { + kEMCalClusterIsMatched = 0, // cluster has at least one matched track (from EMCALMatchedTracks) + kEMCalClusterUserCutsBits = 1 // first bit for the user EMCal cluster selections + }; + public: enum Variables { kNothing = -1, @@ -633,6 +639,22 @@ class VarManager : public TObject kIsDalitzLeg, // Up to 8 dalitz selections kBarrelNAssocsInBunch = kIsDalitzLeg + 8, // number of in bunch collision associations kBarrelNAssocsOutOfBunch, // number of out of bunch collision associations + kEMCalEnergy, // EMCal cluster energy + kEMCalCoreEnergy, // EMCal cluster core energy + kEMCalRawEnergy, // EMCal cluster raw energy + kEMCalEta, // EMCal cluster pseudorapidity + kEMCalPhi, // EMCal cluster azimuthal angle + kEMCalM02, // EMCal shower shape long axis + kEMCalM20, // EMCal shower shape short axis + kEMCalNCells, // number of cells in the EMCal cluster + kEMCalTime, // EMCal cluster time + kEMCalIsExotic, // exotic EMCal cluster flag + kEMCalDistanceToBadChannel, // distance of the EMCal cluster to the closest bad channel + kEMCalNLM, // number of local maxima of the EMCal cluster + kEMCalDefinition, // EMCal cluster definition + kEMCalEoverP, // EMCal cluster energy over matched track momentum + kEMCalMatchDeltaEta, // eta residual of the matched EMCal cluster w.r.t. the track + kEMCalMatchDeltaPhi, // phi residual of the matched EMCal cluster w.r.t. the track kNBarrelTrackVariables, // Muon track variables @@ -1407,6 +1429,8 @@ class VarManager : public TObject static void FillTrack(T const& track, float* values = nullptr); template static void FillPhoton(T const& photon, float* values = nullptr); + template + static void FillTrackEMCal(T const& cluster, float trackP = -1.0f, float deltaEta = -999.0f, float deltaPhi = -999.0f, float* values = nullptr); template static void FillTrackCollision(T const& track, C const& collision, float* values = nullptr); template @@ -3440,6 +3464,35 @@ void VarManager::FillPhoton(T const& track, float* values) } } +template +void VarManager::FillTrackEMCal(T const& cluster, float trackP, float deltaEta, float deltaPhi, float* values) +{ + // Fill the EMCal cluster quantities; the cluster can be either an aod::EMCALCluster (skimming) + // or an aod::ReducedEMCal (analysis on skimmed data), the column schemas are identical. + // trackP: momentum of the matched track, used to compute E/p when called in a track-cluster matching context; + // the default (negative) value leaves E/p at the -999 sentinel, as for tracks without a matched cluster. + if (!values) { + values = fgValues; + } + + values[kEMCalEnergy] = cluster.energy(); + values[kEMCalCoreEnergy] = cluster.coreEnergy(); + values[kEMCalRawEnergy] = cluster.rawEnergy(); + values[kEMCalEta] = cluster.eta(); + values[kEMCalPhi] = cluster.phi(); + values[kEMCalM02] = cluster.m02(); + values[kEMCalM20] = cluster.m20(); + values[kEMCalNCells] = cluster.nCells(); + values[kEMCalTime] = cluster.time(); + values[kEMCalIsExotic] = cluster.isExotic(); + values[kEMCalDistanceToBadChannel] = cluster.distanceToBadChannel(); + values[kEMCalNLM] = cluster.nlm(); + values[kEMCalDefinition] = cluster.definition(); + values[kEMCalEoverP] = (trackP > 0.0f ? cluster.energy() / trackP : -999.0f); + values[kEMCalMatchDeltaEta] = deltaEta; + values[kEMCalMatchDeltaPhi] = deltaPhi; +} + template void VarManager::FillTrackMC(const U& mcStack, T const& track, float* values) { diff --git a/PWGDQ/DataModel/ReducedInfoTables.h b/PWGDQ/DataModel/ReducedInfoTables.h index 781942d5fe7..3e32af90097 100644 --- a/PWGDQ/DataModel/ReducedInfoTables.h +++ b/PWGDQ/DataModel/ReducedInfoTables.h @@ -366,6 +366,39 @@ DECLARE_SOA_TABLE(ReducedFITs, "AOD", "REDUCEDFIT", //! FIT detector information using ReducedFIT = ReducedFITs::iterator; +namespace reducedemcal +{ +// EMCal cluster information (from aod::EMCALClusters, see PWGJE/DataModel/EMCALClusters.h) +DECLARE_SOA_INDEX_COLUMN(ReducedEvent, reducedevent); //! +DECLARE_SOA_BITMAP_COLUMN(FilteringFlags, filteringFlags, 64); //! Bitmap of the cluster selections +DECLARE_SOA_COLUMN(Energy, energy, float); //! cluster energy (GeV) +DECLARE_SOA_COLUMN(CoreEnergy, coreEnergy, float); //! cluster core energy (GeV) +DECLARE_SOA_COLUMN(RawEnergy, rawEnergy, float); //! raw cluster energy (GeV) +DECLARE_SOA_COLUMN(Eta, eta, float); //! cluster pseudorapidity (calculated using vertex) +DECLARE_SOA_COLUMN(Phi, phi, float); //! cluster azimuthal angle (calculated using vertex) +DECLARE_SOA_COLUMN(M02, m02, float); //! shower shape long axis +DECLARE_SOA_COLUMN(M20, m20, float); //! shower shape short axis +DECLARE_SOA_COLUMN(NCells, nCells, int); //! number of cells in cluster +DECLARE_SOA_COLUMN(Time, time, float); //! cluster time (ns) +DECLARE_SOA_COLUMN(IsExotic, isExotic, bool); //! flag to mark cluster as exotic +DECLARE_SOA_COLUMN(DistanceToBadChannel, distanceToBadChannel, float); //! distance to bad channel +DECLARE_SOA_COLUMN(NLM, nlm, int); //! number of local maxima +DECLARE_SOA_COLUMN(Definition, definition, int); //! cluster definition, see EMCALClusterDefinition.h +DECLARE_SOA_COLUMN(IsAmbiguous, isAmbiguous, bool); //! true for clusters from BCs with no or multiple collisions +} // namespace reducedemcal + +// Skimmed EMCal cluster information; standalone table (not joinable to the barrel track tables), +// the matching to barrel tracks is stored as an index column in the ReducedTracksBarrelEMCal table +DECLARE_SOA_TABLE(ReducedEMCals, "AOD", "REDUCEDEMCAL", //! + o2::soa::Index<>, reducedemcal::ReducedEventId, reducedemcal::FilteringFlags, + reducedemcal::Energy, reducedemcal::CoreEnergy, reducedemcal::RawEnergy, + reducedemcal::Eta, reducedemcal::Phi, reducedemcal::M02, reducedemcal::M20, + reducedemcal::NCells, reducedemcal::Time, reducedemcal::IsExotic, + reducedemcal::DistanceToBadChannel, reducedemcal::NLM, reducedemcal::Definition, + reducedemcal::IsAmbiguous); + +using ReducedEMCal = ReducedEMCals::iterator; + namespace reducedtrack { // basic track information @@ -399,6 +432,9 @@ DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, //! [](float pt, float eta) -> float { return pt * std::sinh(eta); }); DECLARE_SOA_DYNAMIC_COLUMN(P, p, //! [](float pt, float eta) -> float { return pt * std::cosh(eta); }); +DECLARE_SOA_INDEX_COLUMN(ReducedEMCal, matchedEMCalCluster); //! matching index pointing to the ReducedEMCals table, -1 if no matched cluster +DECLARE_SOA_COLUMN(EMCalMatchDeltaEta, emcalMatchDeltaEta, float); //! eta residual of the matched EMCal cluster w.r.t. the track +DECLARE_SOA_COLUMN(EMCalMatchDeltaPhi, emcalMatchDeltaPhi, float); //! phi residual of the matched EMCal cluster w.r.t. the track } // namespace reducedtrack // basic track information @@ -447,11 +483,18 @@ DECLARE_SOA_TABLE(ReducedTracksBarrelPID, "AOD", "RTBARRELPID", //! DECLARE_SOA_TABLE(ReducedTracksBarrelInfo, "AOD", "RTBARRELINFO", reducedtrack::CollisionId, collision::PosX, collision::PosY, collision::PosZ, reducedtrack::TrackId); +// barrel track EMCal matching information (joined with ReducedTracks); +// one entry per skimmed barrel track, the index is -1 for tracks without a matched cluster +DECLARE_SOA_TABLE(ReducedTracksBarrelEMCal, "AOD", "RTBARRELEMCAL", //! + reducedtrack::ReducedEMCalId, + reducedtrack::EMCalMatchDeltaEta, reducedtrack::EMCalMatchDeltaPhi); + using ReducedTrack = ReducedTracks::iterator; using ReducedTrackBarrel = ReducedTracksBarrel::iterator; using ReducedTrackBarrelCov = ReducedTracksBarrelCov::iterator; using ReducedTrackBarrelPID = ReducedTracksBarrelPID::iterator; using ReducedTrackBarrelInfo = ReducedTracksBarrelInfo::iterator; +using ReducedTrackBarrelEMCal = ReducedTracksBarrelEMCal::iterator; namespace reducedtrackMC { diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index 32749d294bf..72a31272e1b 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -26,6 +26,7 @@ #include "PWGDQ/Core/MuonMatchingMlResponse.h" #include "PWGDQ/Core/VarManager.h" #include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "PWGJE/DataModel/EMCALClusters.h" #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/RCTSelectionFlags.h" @@ -155,6 +156,7 @@ constexpr static uint32_t gkEventFillMapWithMultsExtra = VarManager::ObjTypes::B // constexpr static uint32_t gkEventFillMapWithCentRun2 = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionCentRun2; // Unused variable // constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackPID | VarManager::ObjTypes::TrackPIDExtra; constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackCov | VarManager::ObjTypes::TrackPID | VarManager::ObjTypes::TrackPIDExtra; +constexpr static uint32_t gkTrackFillMapWithCovWithEMCal = gkTrackFillMapWithCov | VarManager::ObjTypes::TrackEMCal; constexpr static uint32_t gkTrackFillMapWithV0Bits = gkTrackFillMapWithCov | VarManager::ObjTypes::TrackV0Bits; constexpr static uint32_t gkTrackFillMapWithV0BitsNoTOF = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackV0Bits | VarManager::ObjTypes::TrackTPCPID; constexpr static uint32_t gkTrackFillMapNoTOF = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackTPCPID; @@ -176,7 +178,8 @@ enum SkimStatsHists { kStatsMuons, kStatsOrphanTracks, kStatsZorroInfo, - kStatsZorroSel + kStatsZorroSel, + kStatsEMCal }; namespace dqtablemaker_helpers @@ -204,6 +207,8 @@ struct TableMaker { Produces trackBarrelCov; Produces trackBarrelPID; Produces trackBarrelAssoc; + Produces emcal; + Produces trackBarrelEMCal; Produces muonBasic; Produces muonExtra; Produces muonCov; @@ -230,6 +235,7 @@ struct TableMaker { Configurable fConfigEventCutsJSON{"cfgEventCutsJSON", "", "Additional event selection in JSON format"}; Configurable fConfigTrackCutsJSON{"cfgBarrelTrackCutsJSON", "", "Additional list of barrel track cuts in JSON format"}; Configurable fConfigMuonCutsJSON{"cfgMuonCutsJSON", "", "Additional list of muon cuts in JSON format"}; + Configurable fConfigEMCalCutsJSON{"cfgEMCalClusterCutsJSON", "", "List of EMCal cluster cuts in JSON format; if empty, all clusters are written"}; } fConfigCuts; // RCT selection @@ -255,6 +261,7 @@ struct TableMaker { Configurable fConfigAddEventHistogram{"cfgAddEventHistogram", "", "Comma separated list of histograms"}; Configurable fConfigAddTrackHistogram{"cfgAddTrackHistogram", "", "Comma separated list of histograms"}; Configurable fConfigAddMuonHistogram{"cfgAddMuonHistogram", "", "Comma separated list of histograms"}; + Configurable fConfigAddEMCalHistogram{"cfgAddEMCalHistogram", "emcal", "Comma separated list of histogram subgroups for the EMCal cluster classes"}; Configurable fConfigAddJSONHistograms{"cfgAddJSONHistograms", "", "Histograms in JSON format"}; Configurable fConfigIrEstimator{"cfgIrEstimator", "", "Estimator of the interaction rate (pp,OO --> T0VTX, Pb-Pb --> ZNC hadronic), to be used with cfgFillBcStat"}; } fConfigHistOutput; @@ -334,6 +341,7 @@ struct TableMaker { AnalysisCompositeCut* fEventCut = nullptr; //! Event selection cut std::vector fTrackCuts; //! Barrel track cuts std::vector fMuonCuts; //! Muon track cuts + std::vector fEMCalCuts; //! EMCal cluster cuts bool fDoDetailedQA = false; // Bool to set detailed QA true, if QA is set true int fCurrentRun = -1; // needed to detect if the run changed and trigger update of calibrations etc. @@ -345,6 +353,16 @@ struct TableMaker { std::map fFwdTrackIndexMapReversed; // key: new fwd-track global index, value: fwd-track global index std::map fFwdTrackFilterMap; // key: fwd-track global index, value: fwd-track filter map std::map fMftIndexMap; // key: MFT tracklet global index, value: new MFT tracklet global index + std::map fEmcalIndexMap; // key: EMCal cluster global index, value: skimmed cluster index + + // best (smallest delta-R) matched skimmed EMCal cluster for a barrel track + struct EMCalMatch { + int32_t clusterIdx = -1; + float deltaEta = -999.0f; + float deltaPhi = -999.0f; + float deltaR2 = 999.0f; + }; + std::map fTrackEMCalMatchMap; // key: track global index, value: best matched skimmed EMCal cluster std::map fBestMatch; std::unordered_map map_mfttrackcovs; @@ -371,6 +389,8 @@ struct TableMaker { Preslice trackIndicesPerCollision = aod::track_association::collisionId; Preslice fwdtrackIndicesPerCollision = aod::track_association::collisionId; Preslice mfttrackIndicesPerCollision = aod::track_association::collisionId; + Preslice emcalClustersPerCollision = aod::emcalcluster::collisionId; + Preslice emcalMatchedTracksPerCluster = aod::emcalclustercell::emcalclusterId; Preslice preslice = aod::track::collisionId; Partition tracksPos = (((aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)) && (aod::track::tgl > static_cast(0.05))); @@ -484,9 +504,11 @@ struct TableMaker { // Check whether we have to define barrel or muon histograms bool enableBarrelHistos = (context.mOptions.get("processPP") || context.mOptions.get("processPPWithFilter") || context.mOptions.get("processPPWithFilterBarrelOnly") || context.mOptions.get("processPPBarrelOnly") || context.mOptions.get("processPbPb") || context.mOptions.get("processPbPbBarrelOnly") || context.mOptions.get("processPbPbBarrelOnlyWithV0Bits") || context.mOptions.get("processPbPbBarrelOnlyWithV0BitsNoTOF")) || - context.mOptions.get("processPbPbWithFilterBarrelOnly") || context.mOptions.get("processPPBarrelOnlyWithV0s") || context.mOptions.get("processPbPbBarrelOnlyNoTOF"); + context.mOptions.get("processPbPbWithFilterBarrelOnly") || context.mOptions.get("processPPBarrelOnlyWithV0s") || context.mOptions.get("processPbPbBarrelOnlyNoTOF") || context.mOptions.get("processPPBarrelOnlyWithEMCal") || context.mOptions.get("processPPWithEMCal"); + + bool enableEMCalHistos = (context.mOptions.get("processPPBarrelOnlyWithEMCal") || context.mOptions.get("processPPWithEMCal")); - bool enableMuonHistos = (context.mOptions.get("processPP") || context.mOptions.get("processPPWithFilter") || context.mOptions.get("processPPWithFilterMuonOnly") || context.mOptions.get("processPPWithFilterMuonMFT") || context.mOptions.get("processPPMuonOnly") || context.mOptions.get("processPPRealignedMuonOnly") || context.mOptions.get("processPPMuonMFT") || context.mOptions.get("processPPMuonMFTWithMultsExtra") || + bool enableMuonHistos = (context.mOptions.get("processPP") || context.mOptions.get("processPPWithEMCal") || context.mOptions.get("processPPWithFilter") || context.mOptions.get("processPPWithFilterMuonOnly") || context.mOptions.get("processPPWithFilterMuonMFT") || context.mOptions.get("processPPMuonOnly") || context.mOptions.get("processPPRealignedMuonOnly") || context.mOptions.get("processPPMuonMFT") || context.mOptions.get("processPPMuonMFTWithMultsExtra") || context.mOptions.get("processPbPb") || context.mOptions.get("processPbPbMuonOnly") || context.mOptions.get("processPbPbWithFilterMuonOnly") || context.mOptions.get("processPbPbStreamMuonOnly") || context.mOptions.get("processPbPbRealignedMuonOnly") || context.mOptions.get("processPbPbMuonMFT")); if (enableBarrelHistos) { @@ -507,6 +529,18 @@ struct TableMaker { histClasses += "TrackBarrel_PostCalibProton;"; } } + if (enableEMCalHistos) { + // EMCal cluster histograms, before selections + if (fDoDetailedQA) { + histClasses += "EMCalClusters_BeforeCuts;"; + } + if (fConfigHistOutput.fConfigQA) { + // EMCal cluster histograms after selections; one histogram directory for each user specified selection + for (const auto& cut : fEMCalCuts) { + histClasses += Form("EMCalClusters_%s;", cut->GetName()); + } + } + } if (enableMuonHistos) { // Muon tracks before cuts and MFT tracks if (fDoDetailedQA) { @@ -583,6 +617,15 @@ struct TableMaker { } } + // EMCal cluster cuts, via JSON only (there are no named EMCal cluster cuts in the CutsLibrary yet) + TString addEMCalCutsStr = fConfigCuts.fConfigEMCalCutsJSON.value; + if (addEMCalCutsStr != "") { + std::vector addEMCalCuts = dqcuts::GetCutsFromJSON(addEMCalCutsStr.Data()); + for (const auto& t : addEMCalCuts) { + fEMCalCuts.push_back(dynamic_cast(t)); + } + } + // Muon cuts cutNamesStr = fConfigCuts.fConfigMuonCuts.value; if (!cutNamesStr.IsNull()) { @@ -647,6 +690,13 @@ struct TableMaker { } } + TString histEMCalName = fConfigHistOutput.fConfigAddEMCalHistogram.value; + if (classStr.Contains("EMCalClusters")) { + if (fConfigHistOutput.fConfigQA) { + dqhistograms::DefineHistograms(fHistMan, objArray->At(iclass)->GetName(), "track", histEMCalName); + } + } + TString histMftName = fConfigHistOutput.fConfigAddMuonHistogram.value; if (classStr.Contains("Mft")) { if (fConfigHistOutput.fConfigDetailedQA) { @@ -714,6 +764,15 @@ struct TableMaker { TH2D* histZorroSel = new TH2D("ZorroSel", "trigger of interested", 1, -0.5, 0.5, 1, -0.5, 0.5); fStatsList->AddAt(histZorroSel, kStatsZorroSel); + + // EMCal cluster statistics: one bin per cluster selection, plus one bin counting all skimmed clusters + TH1D* histEMCal = new TH1D("EMCalClusterStats", "EMCal cluster statistics", fEMCalCuts.size() + 1, -0.5, fEMCalCuts.size() + 0.5); + ib = 1; + for (auto cut = fEMCalCuts.begin(); cut != fEMCalCuts.end(); cut++, ib++) { + histEMCal->GetXaxis()->SetBinLabel(ib, (*cut)->GetName()); + } + histEMCal->GetXaxis()->SetBinLabel(fEMCalCuts.size() + 1, "Skimmed clusters"); + fStatsList->AddAt(histEMCal, kStatsEMCal); } template @@ -1423,6 +1482,17 @@ struct TableMaker { -999.0); } + // write the matched EMCal cluster index (or -1 if the track has no matched cluster); + // this table has exactly one entry per skimmed track, so it is joinable to the ReducedTracks tables + if constexpr (static_cast(TTrackFillMap & VarManager::ObjTypes::TrackEMCal)) { + auto match = fTrackEMCalMatchMap.find(track.globalIndex()); + if (match != fTrackEMCalMatchMap.end()) { + outTables.trackBarrelEMCal(match->second.clusterIdx, match->second.deltaEta, match->second.deltaPhi); + } else { + outTables.trackBarrelEMCal(-1, -999.0f, -999.0f); + } + } + fTrackIndexMap[track.globalIndex()] = outTables.trackBasic.lastIndex(); // write the skimmed collision - track association @@ -1430,6 +1500,71 @@ struct TableMaker { } // end loop over associations } // end skimTracks + template + void skimEMCal(TEMCals const& clusters, TEMCalMatches const& matchedTracks) + { + // Skim the EMCal clusters of all the selected collisions and record, for each track matched to a cluster, + // the closest (smallest delta-R) skimmed cluster. This has to run before skimTracks, such that the + // matching indices are available at the time the tracks are written. + // NOTE: only the unique clusters (aod::EMCALClusters) are skimmed here; the ambiguous clusters + // (aod::EMCALAmbiguousClusters, from BCs with none or multiple reconstructed collisions) are not + // associated to a collision and have no track matching, so they are left for a dedicated treatment. + for (auto const& [origIdx, skimIdx] : fCollIndexMap) { + auto groupedClusters = clusters.sliceBy(emcalClustersPerCollision, origIdx); + for (const auto& cluster : groupedClusters) { + VarManager::FillTrackEMCal(cluster); + if (fDoDetailedQA) { + fHistMan->FillHistClass("EMCalClusters_BeforeCuts", dqtablemaker_helpers::varValues()); + } + + // apply the cluster selections and fill the filtering tag and the stats histogram + auto clusterTempFilterMap = static_cast(0); + int i = 0; + for (auto cut = fEMCalCuts.begin(); cut != fEMCalCuts.end(); cut++, i++) { + if ((*cut)->IsSelected(dqtablemaker_helpers::varValues())) { + clusterTempFilterMap |= (static_cast(1) << i); + if (fConfigHistOutput.fConfigQA) { + fHistMan->FillHistClass(Form("EMCalClusters_%s", (*cut)->GetName()), dqtablemaker_helpers::varValues()); + } + (dynamic_cast(fStatsList->At(kStatsEMCal)))->Fill(static_cast(i)); + } + } + // if cluster selections are specified, write only the clusters which fulfill at least one of them + if (!fEMCalCuts.empty() && clusterTempFilterMap == 0) { + continue; + } + (dynamic_cast(fStatsList->At(kStatsEMCal)))->Fill(static_cast(fEMCalCuts.size())); + auto clusterFilteringTag = (static_cast(clusterTempFilterMap) << VarManager::kEMCalClusterUserCutsBits); + + // flag the clusters which have at least one matched track, such that a clean + // "EMCal only" cluster sample (e.g. for photons) can be selected on the skimmed data + auto clusterMatches = matchedTracks.sliceBy(emcalMatchedTracksPerCluster, cluster.globalIndex()); + if (clusterMatches.size() > 0) { + clusterFilteringTag |= (static_cast(1) << VarManager::kEMCalClusterIsMatched); + } + + outTables.emcal(skimIdx, clusterFilteringTag, + cluster.energy(), cluster.coreEnergy(), cluster.rawEnergy(), + cluster.eta(), cluster.phi(), cluster.m02(), cluster.m20(), + cluster.nCells(), cluster.time(), cluster.isExotic(), + cluster.distanceToBadChannel(), cluster.nlm(), cluster.definition(), + false); + fEmcalIndexMap[cluster.globalIndex()] = outTables.emcal.lastIndex(); + + // record the track matches of this cluster; for each track keep only the closest cluster + for (const auto& match : clusterMatches) { + float deltaEta = match.deltaEta(); + float deltaPhi = match.deltaPhi(); + float deltaR2 = deltaEta * deltaEta + deltaPhi * deltaPhi; + auto existing = fTrackEMCalMatchMap.find(match.trackId()); + if (existing == fTrackEMCalMatchMap.end() || deltaR2 < existing->second.deltaR2) { + fTrackEMCalMatchMap[match.trackId()] = EMCalMatch{static_cast(outTables.emcal.lastIndex()), deltaEta, deltaPhi, deltaR2}; + } + } + } // end loop over clusters + } // end loop over skimmed collisions + } // end skimEMCal + template void skimMFT(TEvent const& collision, TBCs const& /*bcs*/, MFTTracks const& /*mfts*/, MFTTrackAssoc const& mftAssocs) { @@ -1689,7 +1824,8 @@ struct TableMaker { typename TTrackAssoc, typename TFwdTrackAssoc, typename TMFTTrackAssoc, typename TMFTCov, typename TFt0s, typename TFv0as, typename TFdds> void fullSkimming(TEvents const& collisions, TBCs const& bcs, TZdcs const& zdcs, TTracks const& tracksBarrel, TMuons const& muons, TMFTTracks const& mftTracks, - TTrackAssoc const& trackAssocs, TFwdTrackAssoc const& fwdTrackAssocs, TMFTTrackAssoc const& mftAssocs, TMFTCov const& mftCovs, TFt0s const& ft0s, TFv0as const& fv0as, TFdds const& fdds) + TTrackAssoc const& trackAssocs, TFwdTrackAssoc const& fwdTrackAssocs, TMFTTrackAssoc const& mftAssocs, TMFTCov const& mftCovs, TFt0s const& ft0s, TFv0as const& fv0as, TFdds const& fdds, + aod::EMCALClusters const* emcalClusters = nullptr, aod::EMCALMatchedTracks const* emcalMatchedTracks = nullptr) { if (bcs.size() > 0 && fCurrentRun != bcs.begin().runNumber()) { @@ -1772,6 +1908,16 @@ struct TableMaker { outTables.trackBarrelAssoc.reserve(trackAssocs.size()); } + // skim the EMCal clusters of the selected collisions; this runs before the track skimming + // such that the track to cluster matching indices are available when the tracks are written + if constexpr (static_cast(TTrackFillMap & VarManager::ObjTypes::TrackEMCal)) { + fEmcalIndexMap.clear(); + fTrackEMCalMatchMap.clear(); + outTables.emcal.reserve(emcalClusters->size()); + outTables.trackBarrelEMCal.reserve(tracksBarrel.size()); + skimEMCal(*emcalClusters, *emcalMatchedTracks); + } + if constexpr (static_cast(TMFTFillMap)) { fMftIndexMap.clear(); map_mfttrackcovs.clear(); @@ -1860,6 +2006,17 @@ struct TableMaker { fullSkimming(collisions, bcs, nullptr, tracksBarrel, muons, nullptr, trackAssocs, fwdTrackAssocs, nullptr, nullptr, nullptr, nullptr, nullptr); } + // produce the full (barrel + muon) DQ skimmed data model with EMCal clusters, typically for pp (e.g. for e-mu correlations with EMCal electron ID); + // requires the emcal-correction-task upstream in the workflow (for the EMCALClusters and EMCALMatchedTracks tables) + void processPPWithEMCal(MyEventsWithMultsExtra const& collisions, MyBCs const& bcs, + MyBarrelTracksWithCov const& tracksBarrel, + MyMuonsWithCov const& muons, + aod::EMCALClusters const& emcalClusters, aod::EMCALMatchedTracks const& emcalMatchedTracks, + TrackAssoc const& trackAssocs, FwdTrackAssoc const& fwdTrackAssocs) + { + fullSkimming(collisions, bcs, nullptr, tracksBarrel, muons, nullptr, trackAssocs, fwdTrackAssocs, nullptr, nullptr, nullptr, nullptr, nullptr, &emcalClusters, &emcalMatchedTracks); + } + // produce the barrel-only DQ skimmed data model typically for pp/p-Pb or UPC Pb-Pb (no centrality), subscribe to the DQ event filter (filter-pp or filter-PbPb) void processPPWithFilterBarrelOnly(MyEventsWithMultsAndFilter const& collisions, MyBCs const& bcs, aod::Zdcs const& zdcs, MyBarrelTracksWithCov const& tracksBarrel, @@ -1891,6 +2048,16 @@ struct TableMaker { fullSkimming(collisions, bcs, zdcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr); } + // produce the barrel-only DQ skimmed data model with EMCal clusters, typically for pp (no centrality); + // requires the emcal-correction-task upstream in the workflow (for the EMCALClusters and EMCALMatchedTracks tables) + void processPPBarrelOnlyWithEMCal(MyEventsWithMults const& collisions, MyBCs const& bcs, aod::Zdcs const& zdcs, + MyBarrelTracksWithCov const& tracksBarrel, + aod::EMCALClusters const& emcalClusters, aod::EMCALMatchedTracks const& emcalMatchedTracks, + TrackAssoc const& trackAssocs) + { + fullSkimming(collisions, bcs, zdcs, tracksBarrel, nullptr, nullptr, trackAssocs, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, &emcalClusters, &emcalMatchedTracks); + } + // produce the barrel-only DQ skimmed barrel data model, with V0 tagged tracks void processPPBarrelOnlyWithV0s(MyEventsWithMults const& collisions, MyBCs const& bcs, MyBarrelTracksWithV0BitsNoTOF const& tracksBarrel, @@ -2047,11 +2214,13 @@ struct TableMaker { } PROCESS_SWITCH(TableMaker, processPP, "Build full DQ skimmed data model for pp/p-Pb w/o event filtering (use Zorro)", false); + PROCESS_SWITCH(TableMaker, processPPWithEMCal, "Build full DQ skimmed data model for pp/p-Pb w/o event filtering, with EMCal clusters", false); PROCESS_SWITCH(TableMaker, processPPWithFilter, "Build full DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb, w/ event filtering", false); PROCESS_SWITCH(TableMaker, processPPWithFilterBarrelOnly, "Build barrel only DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb, w/ event filtering", false); PROCESS_SWITCH(TableMaker, processPPWithFilterMuonOnly, "Build muon only DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb, w/ event filtering", false); PROCESS_SWITCH(TableMaker, processPPWithFilterMuonMFT, "Build muon + mft DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb, w/ event filtering", false); PROCESS_SWITCH(TableMaker, processPPBarrelOnly, "Build barrel only DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb", false); + PROCESS_SWITCH(TableMaker, processPPBarrelOnlyWithEMCal, "Build barrel only DQ skimmed data model with EMCal clusters, typically for pp", false); PROCESS_SWITCH(TableMaker, processPPBarrelOnlyWithV0s, "Build barrel only DQ skimmed data model, pp like, with V0 tagged tracks", false); PROCESS_SWITCH(TableMaker, processPPMuonOnly, "Build muon only DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb", false); PROCESS_SWITCH(TableMaker, processPPRealignedMuonOnly, "Build realigned muon only DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb", false); diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index eaaea676093..3b339ce0d33 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -247,6 +247,8 @@ using MyBarrelTracksWithAmbiguities = soa::Join; using MyBarrelTracksWithCovWithAmbiguities = soa::Join; using MyBarrelTracksWithCovWithAmbiguitiesWithColl = soa::Join; +using MyBarrelTracksWithEMCal = soa::Join; +using MyBarrelTracksWithCovWithEMCal = soa::Join; using MyDielectronCandidates = soa::Join; using MyDitrackCandidates = soa::Join; using MyDimuonCandidates = soa::Join; @@ -274,6 +276,8 @@ constexpr static uint32_t gkEventFillMapWithQvectorCentr = VarManager::ObjTypes: constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID; constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID; constexpr static uint32_t gkTrackFillMapWithCovWithColl = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID | VarManager::ObjTypes::ReducedTrackCollInfo; +constexpr static uint32_t gkTrackFillMapWithEMCal = gkTrackFillMap | VarManager::ObjTypes::TrackEMCal; +constexpr static uint32_t gkTrackFillMapWithCovWithEMCal = gkTrackFillMapWithCov | VarManager::ObjTypes::TrackEMCal; // constexpr static uint32_t gkMuonFillMap = VarManager::ObjTypes::ReducedMuon | VarManager::ObjTypes::ReducedMuonExtra; constexpr static uint32_t gkMuonFillMapWithCov = VarManager::ObjTypes::ReducedMuon | VarManager::ObjTypes::ReducedMuonExtra | VarManager::ObjTypes::ReducedMuonCov; @@ -797,6 +801,14 @@ struct AnalysisTrackSelection { if (fPropTrack) { VarManager::FillTrackCollision(track, event); } + // fill the quantities of the matched EMCal cluster (e.g. E/p); + // tracks without a matched cluster keep the sentinel values set by ResetValues + if constexpr (static_cast(TTrackFillMap & VarManager::ObjTypes::TrackEMCal)) { + if (track.has_matchedEMCalCluster()) { + auto cluster = track.template matchedEMCalCluster_as(); + VarManager::FillTrackEMCal(cluster, track.p(), track.emcalMatchDeltaEta(), track.emcalMatchDeltaPhi()); + } + } if (fConfigQA) { fHistMan->FillHistClass("TrackBarrel_BeforeCuts", dqtablereader_helpers::varValues()); } @@ -892,6 +904,14 @@ struct AnalysisTrackSelection { { runTrackSelection(assocs, events, tracks); } + void processSkimmedWithEMCal(ReducedTracksAssoc const& assocs, MyEventsSelected const& events, MyBarrelTracksWithEMCal const& tracks, aod::ReducedEMCals const& /*emcals*/) + { + runTrackSelection(assocs, events, tracks); + } + void processSkimmedWithCovWithEMCal(ReducedTracksAssoc const& assocs, MyEventsVtxCovSelected const& events, MyBarrelTracksWithCovWithEMCal const& tracks, aod::ReducedEMCals const& /*emcals*/) + { + runTrackSelection(assocs, events, tracks); + } void processDummy(MyEventsBasic const&) { // do nothing @@ -900,6 +920,8 @@ struct AnalysisTrackSelection { PROCESS_SWITCH(AnalysisTrackSelection, processSkimmed, "Run barrel track selection on DQ skimmed track associations", false); PROCESS_SWITCH(AnalysisTrackSelection, processSkimmedWithMultExtra, "Run barrel track selection on DQ skimmed track associations, with extra multiplicity tables", false); PROCESS_SWITCH(AnalysisTrackSelection, processSkimmedWithCov, "Run barrel track selection on DQ skimmed tracks w/ cov matrix associations", false); + PROCESS_SWITCH(AnalysisTrackSelection, processSkimmedWithEMCal, "Run barrel track selection on DQ skimmed track associations, with EMCal matching info", false); + PROCESS_SWITCH(AnalysisTrackSelection, processSkimmedWithCovWithEMCal, "Run barrel track selection on DQ skimmed tracks w/ cov matrix associations, with EMCal matching info", false); PROCESS_SWITCH(AnalysisTrackSelection, processDummy, "Dummy function", true); };