From b71a8f8d163af5cb58f7b5651807038f81ac2a18 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Fri, 13 Dec 2024 11:25:29 +0100 Subject: [PATCH 1/2] [DPG] [TOF] add plots for delta (#8764) --- DPG/Tasks/AOTTrack/PID/TOF/qaPIDTOFMC.cxx | 447 +++++++++++----------- 1 file changed, 226 insertions(+), 221 deletions(-) diff --git a/DPG/Tasks/AOTTrack/PID/TOF/qaPIDTOFMC.cxx b/DPG/Tasks/AOTTrack/PID/TOF/qaPIDTOFMC.cxx index e841e057497..73d1c72ad1e 100644 --- a/DPG/Tasks/AOTTrack/PID/TOF/qaPIDTOFMC.cxx +++ b/DPG/Tasks/AOTTrack/PID/TOF/qaPIDTOFMC.cxx @@ -12,13 +12,14 @@ /// /// \file qaPIDTOFMC.cxx /// \author Nicolò Jacazio -/// \brief Task to produce QA output of the PID with TOF running on the MC. +/// \brief Task to produce QA output of the PID with TOF running on the MC e.g. to compute purity. /// // O2 includes #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "Framework/StaticFor.h" +#include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponse.h" #include "Framework/runDataProcessing.h" @@ -27,169 +28,71 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::track; +static constexpr int Np = 9; +static constexpr int NpNp = Np * Np; + +std::array, Np> hParticlePt; +std::array, Np> hParticleP; +std::array, Np> hParticleEta; + +std::array, Np> hTrackPt; +std::array, Np> hTrackP; +std::array, Np> hTrackEta; +std::array, Np> hTrackLength; + +std::array, Np> hSignalMC; +std::array, Np> hSignalMCprm; +std::array, Np> hSignalMCstr; +std::array, Np> hSignalMCmat; + +std::array, Np> hNSigma; +std::array, Np> hNSigmaprm; +std::array, Np> hNSigmastr; +std::array, Np> hNSigmamat; + +std::array, NpNp> hNSigmaMC; +std::array, NpNp> hNSigmaMCprm; +std::array, NpNp> hNSigmaMCstr; +std::array, NpNp> hNSigmaMCmat; + +std::array, NpNp> hDeltaMCEvTime; +std::array, NpNp> hDeltaMCEvTimeTrueGoodEv; +std::array, NpNp> hDeltaMCEvTimeTrueBadEv; +std::array, NpNp> hDeltaMCEvTimeprm; +std::array, NpNp> hDeltaMCEvTimestr; +std::array, NpNp> hDeltaMCEvTimemat; + +std::array, NpNp> hDeltaMCEvTimeMC; +std::array, NpNp> hDeltaMCEvTimeMCprm; +std::array, NpNp> hDeltaMCEvTimeMCstr; +std::array, NpNp> hDeltaMCEvTimeMCmat; +std::array, NpNp> hDeltaMCEvTimeMCGoodMatch; +std::array, NpNp> hDeltaMCEvTimeMCBadMatch; + +std::array, NpNp> hDeltaExpTimeMC; + +template +double calculateTimeOfFlight(double momentum, double length) +{ + static constexpr float mass = pid_constants::sMasses[id]; + static constexpr float c = o2::constants::physics::LightSpeedCm2PS; + // Calculate the Lorentz factor gamma + const float gamma = std::sqrt(1 + std::pow(momentum / (mass * c), 2)); + + // Calculate velocity + double velocity = momentum / (gamma * mass); + + // Calculate time of flight + double timeOfFlight = length / velocity; + + return timeOfFlight; +} + /// Task to produce the TOF QA plots struct pidTofQaMc { SliceCache cache; - - static constexpr int Np = 9; - static constexpr int NpNp = Np * Np; - static constexpr std::string_view hparticlept[Np] = {"particlept/El", "particlept/Mu", "particlept/Pi", - "particlept/Ka", "particlept/Pr", "particlept/De", - "particlept/Tr", "particlept/He", "particlept/Al"}; - static constexpr std::string_view hparticlep[Np] = {"particlep/El", "particlep/Mu", "particlep/Pi", - "particlep/Ka", "particlep/Pr", "particlep/De", - "particlep/Tr", "particlep/He", "particlep/Al"}; - static constexpr std::string_view hparticleeta[Np] = {"particleeta/El", "particleeta/Mu", "particleeta/Pi", - "particleeta/Ka", "particleeta/Pr", "particleeta/De", - "particleeta/Tr", "particleeta/He", "particleeta/Al"}; - static constexpr std::string_view htrackpt[Np] = {"trackpt/El", "trackpt/Mu", "trackpt/Pi", - "trackpt/Ka", "trackpt/Pr", "trackpt/De", - "trackpt/Tr", "trackpt/He", "trackpt/Al"}; - static constexpr std::string_view htrackp[Np] = {"trackp/El", "trackp/Mu", "trackp/Pi", - "trackp/Ka", "trackp/Pr", "trackp/De", - "trackp/Tr", "trackp/He", "trackp/Al"}; - static constexpr std::string_view htracketa[Np] = {"tracketa/El", "tracketa/Mu", "tracketa/Pi", - "tracketa/Ka", "tracketa/Pr", "tracketa/De", - "tracketa/Tr", "tracketa/He", "tracketa/Al"}; - static constexpr std::string_view htracklength[Np] = {"tracklength/El", "tracklength/Mu", "tracklength/Pi", - "tracklength/Ka", "tracklength/Pr", "tracklength/De", - "tracklength/Tr", "tracklength/He", "tracklength/Al"}; - // Signal - static constexpr std::string_view hsignalMC[Np] = {"signalMC/El", "signalMC/Mu", "signalMC/Pi", - "signalMC/Ka", "signalMC/Pr", "signalMC/De", - "signalMC/Tr", "signalMC/He", "signalMC/Al"}; - static constexpr std::string_view hsignalMCprm[Np] = {"signalMCprm/El", "signalMCprm/Mu", "signalMCprm/Pi", - "signalMCprm/Ka", "signalMCprm/Pr", "signalMCprm/De", - "signalMCprm/Tr", "signalMCprm/He", "signalMCprm/Al"}; - static constexpr std::string_view hsignalMCstr[Np] = {"signalMCstr/El", "signalMCstr/Mu", "signalMCstr/Pi", - "signalMCstr/Ka", "signalMCstr/Pr", "signalMCstr/De", - "signalMCstr/Tr", "signalMCstr/He", "signalMCstr/Al"}; - static constexpr std::string_view hsignalMCmat[Np] = {"signalMCmat/El", "signalMCmat/Mu", "signalMCmat/Pi", - "signalMCmat/Ka", "signalMCmat/Pr", "signalMCmat/De", - "signalMCmat/Tr", "signalMCmat/He", "signalMCmat/Al"}; - // Nsigma - static constexpr std::string_view hnsigma[Np] = {"nsigma/El", "nsigma/Mu", "nsigma/Pi", - "nsigma/Ka", "nsigma/Pr", "nsigma/De", - "nsigma/Tr", "nsigma/He", "nsigma/Al"}; - static constexpr std::string_view hnsigmaprm[Np] = {"nsigmaprm/El", "nsigmaprm/Mu", "nsigmaprm/Pi", - "nsigmaprm/Ka", "nsigmaprm/Pr", "nsigmaprm/De", - "nsigmaprm/Tr", "nsigmaprm/He", "nsigmaprm/Al"}; - static constexpr std::string_view hnsigmastr[Np] = {"nsigmastr/El", "nsigmastr/Mu", "nsigmastr/Pi", - "nsigmastr/Ka", "nsigmastr/Pr", "nsigmastr/De", - "nsigmastr/Tr", "nsigmastr/He", "nsigmastr/Al"}; - static constexpr std::string_view hnsigmamat[Np] = {"nsigmamat/El", "nsigmamat/Mu", "nsigmamat/Pi", - "nsigmamat/Ka", "nsigmamat/Pr", "nsigmamat/De", - "nsigmamat/Tr", "nsigmamat/He", "nsigmamat/Al"}; - static constexpr std::string_view hnsigmaMC[NpNp] = {"nsigmaMC/El/El", "nsigmaMC/El/Mu", "nsigmaMC/El/Pi", - "nsigmaMC/El/Ka", "nsigmaMC/El/Pr", "nsigmaMC/El/De", - "nsigmaMC/El/Tr", "nsigmaMC/El/He", "nsigmaMC/El/Al", - "nsigmaMC/Mu/El", "nsigmaMC/Mu/Mu", "nsigmaMC/Mu/Pi", - "nsigmaMC/Mu/Ka", "nsigmaMC/Mu/Pr", "nsigmaMC/Mu/De", - "nsigmaMC/Mu/Tr", "nsigmaMC/Mu/He", "nsigmaMC/Mu/Al", - "nsigmaMC/Pi/El", "nsigmaMC/Pi/Mu", "nsigmaMC/Pi/Pi", - "nsigmaMC/Pi/Ka", "nsigmaMC/Pi/Pr", "nsigmaMC/Pi/De", - "nsigmaMC/Pi/Tr", "nsigmaMC/Pi/He", "nsigmaMC/Pi/Al", - "nsigmaMC/Ka/El", "nsigmaMC/Ka/Mu", "nsigmaMC/Ka/Pi", - "nsigmaMC/Ka/Ka", "nsigmaMC/Ka/Pr", "nsigmaMC/Ka/De", - "nsigmaMC/Ka/Tr", "nsigmaMC/Ka/He", "nsigmaMC/Ka/Al", - "nsigmaMC/Pr/El", "nsigmaMC/Pr/Mu", "nsigmaMC/Pr/Pi", - "nsigmaMC/Pr/Ka", "nsigmaMC/Pr/Pr", "nsigmaMC/Pr/De", - "nsigmaMC/Pr/Tr", "nsigmaMC/Pr/He", "nsigmaMC/Pr/Al", - "nsigmaMC/De/El", "nsigmaMC/De/Mu", "nsigmaMC/De/Pi", - "nsigmaMC/De/Ka", "nsigmaMC/De/Pr", "nsigmaMC/De/De", - "nsigmaMC/De/Tr", "nsigmaMC/De/He", "nsigmaMC/De/Al", - "nsigmaMC/Tr/El", "nsigmaMC/Tr/Mu", "nsigmaMC/Tr/Pi", - "nsigmaMC/Tr/Ka", "nsigmaMC/Tr/Pr", "nsigmaMC/Tr/De", - "nsigmaMC/Tr/Tr", "nsigmaMC/Tr/He", "nsigmaMC/Tr/Al", - "nsigmaMC/He/El", "nsigmaMC/He/Mu", "nsigmaMC/He/Pi", - "nsigmaMC/He/Ka", "nsigmaMC/He/Pr", "nsigmaMC/He/De", - "nsigmaMC/He/Tr", "nsigmaMC/He/He", "nsigmaMC/He/Al", - "nsigmaMC/Al/El", "nsigmaMC/Al/Mu", "nsigmaMC/Al/Pi", - "nsigmaMC/Al/Ka", "nsigmaMC/Al/Pr", "nsigmaMC/Al/De", - "nsigmaMC/Al/Tr", "nsigmaMC/Al/He", "nsigmaMC/Al/Al"}; - static constexpr std::string_view hnsigmaMCstr[NpNp] = {"nsigmaMCstr/El/El", "nsigmaMCstr/El/Mu", "nsigmaMCstr/El/Pi", - "nsigmaMCstr/El/Ka", "nsigmaMCstr/El/Pr", "nsigmaMCstr/El/De", - "nsigmaMCstr/El/Tr", "nsigmaMCstr/El/He", "nsigmaMCstr/El/Al", - "nsigmaMCstr/Mu/El", "nsigmaMCstr/Mu/Mu", "nsigmaMCstr/Mu/Pi", - "nsigmaMCstr/Mu/Ka", "nsigmaMCstr/Mu/Pr", "nsigmaMCstr/Mu/De", - "nsigmaMCstr/Mu/Tr", "nsigmaMCstr/Mu/He", "nsigmaMCstr/Mu/Al", - "nsigmaMCstr/Pi/El", "nsigmaMCstr/Pi/Mu", "nsigmaMCstr/Pi/Pi", - "nsigmaMCstr/Pi/Ka", "nsigmaMCstr/Pi/Pr", "nsigmaMCstr/Pi/De", - "nsigmaMCstr/Pi/Tr", "nsigmaMCstr/Pi/He", "nsigmaMCstr/Pi/Al", - "nsigmaMCstr/Ka/El", "nsigmaMCstr/Ka/Mu", "nsigmaMCstr/Ka/Pi", - "nsigmaMCstr/Ka/Ka", "nsigmaMCstr/Ka/Pr", "nsigmaMCstr/Ka/De", - "nsigmaMCstr/Ka/Tr", "nsigmaMCstr/Ka/He", "nsigmaMCstr/Ka/Al", - "nsigmaMCstr/Pr/El", "nsigmaMCstr/Pr/Mu", "nsigmaMCstr/Pr/Pi", - "nsigmaMCstr/Pr/Ka", "nsigmaMCstr/Pr/Pr", "nsigmaMCstr/Pr/De", - "nsigmaMCstr/Pr/Tr", "nsigmaMCstr/Pr/He", "nsigmaMCstr/Pr/Al", - "nsigmaMCstr/De/El", "nsigmaMCstr/De/Mu", "nsigmaMCstr/De/Pi", - "nsigmaMCstr/De/Ka", "nsigmaMCstr/De/Pr", "nsigmaMCstr/De/De", - "nsigmaMCstr/De/Tr", "nsigmaMCstr/De/He", "nsigmaMCstr/De/Al", - "nsigmaMCstr/Tr/El", "nsigmaMCstr/Tr/Mu", "nsigmaMCstr/Tr/Pi", - "nsigmaMCstr/Tr/Ka", "nsigmaMCstr/Tr/Pr", "nsigmaMCstr/Tr/De", - "nsigmaMCstr/Tr/Tr", "nsigmaMCstr/Tr/He", "nsigmaMCstr/Tr/Al", - "nsigmaMCstr/He/El", "nsigmaMCstr/He/Mu", "nsigmaMCstr/He/Pi", - "nsigmaMCstr/He/Ka", "nsigmaMCstr/He/Pr", "nsigmaMCstr/He/De", - "nsigmaMCstr/He/Tr", "nsigmaMCstr/He/He", "nsigmaMCstr/He/Al", - "nsigmaMCstr/Al/El", "nsigmaMCstr/Al/Mu", "nsigmaMCstr/Al/Pi", - "nsigmaMCstr/Al/Ka", "nsigmaMCstr/Al/Pr", "nsigmaMCstr/Al/De", - "nsigmaMCstr/Al/Tr", "nsigmaMCstr/Al/He", "nsigmaMCstr/Al/Al"}; - static constexpr std::string_view hnsigmaMCmat[NpNp] = {"nsigmaMCmat/El/El", "nsigmaMCmat/El/Mu", "nsigmaMCmat/El/Pi", - "nsigmaMCmat/El/Ka", "nsigmaMCmat/El/Pr", "nsigmaMCmat/El/De", - "nsigmaMCmat/El/Tr", "nsigmaMCmat/El/He", "nsigmaMCmat/El/Al", - "nsigmaMCmat/Mu/El", "nsigmaMCmat/Mu/Mu", "nsigmaMCmat/Mu/Pi", - "nsigmaMCmat/Mu/Ka", "nsigmaMCmat/Mu/Pr", "nsigmaMCmat/Mu/De", - "nsigmaMCmat/Mu/Tr", "nsigmaMCmat/Mu/He", "nsigmaMCmat/Mu/Al", - "nsigmaMCmat/Pi/El", "nsigmaMCmat/Pi/Mu", "nsigmaMCmat/Pi/Pi", - "nsigmaMCmat/Pi/Ka", "nsigmaMCmat/Pi/Pr", "nsigmaMCmat/Pi/De", - "nsigmaMCmat/Pi/Tr", "nsigmaMCmat/Pi/He", "nsigmaMCmat/Pi/Al", - "nsigmaMCmat/Ka/El", "nsigmaMCmat/Ka/Mu", "nsigmaMCmat/Ka/Pi", - "nsigmaMCmat/Ka/Ka", "nsigmaMCmat/Ka/Pr", "nsigmaMCmat/Ka/De", - "nsigmaMCmat/Ka/Tr", "nsigmaMCmat/Ka/He", "nsigmaMCmat/Ka/Al", - "nsigmaMCmat/Pr/El", "nsigmaMCmat/Pr/Mu", "nsigmaMCmat/Pr/Pi", - "nsigmaMCmat/Pr/Ka", "nsigmaMCmat/Pr/Pr", "nsigmaMCmat/Pr/De", - "nsigmaMCmat/Pr/Tr", "nsigmaMCmat/Pr/He", "nsigmaMCmat/Pr/Al", - "nsigmaMCmat/De/El", "nsigmaMCmat/De/Mu", "nsigmaMCmat/De/Pi", - "nsigmaMCmat/De/Ka", "nsigmaMCmat/De/Pr", "nsigmaMCmat/De/De", - "nsigmaMCmat/De/Tr", "nsigmaMCmat/De/He", "nsigmaMCmat/De/Al", - "nsigmaMCmat/Tr/El", "nsigmaMCmat/Tr/Mu", "nsigmaMCmat/Tr/Pi", - "nsigmaMCmat/Tr/Ka", "nsigmaMCmat/Tr/Pr", "nsigmaMCmat/Tr/De", - "nsigmaMCmat/Tr/Tr", "nsigmaMCmat/Tr/He", "nsigmaMCmat/Tr/Al", - "nsigmaMCmat/He/El", "nsigmaMCmat/He/Mu", "nsigmaMCmat/He/Pi", - "nsigmaMCmat/He/Ka", "nsigmaMCmat/He/Pr", "nsigmaMCmat/He/De", - "nsigmaMCmat/He/Tr", "nsigmaMCmat/He/He", "nsigmaMCmat/He/Al", - "nsigmaMCmat/Al/El", "nsigmaMCmat/Al/Mu", "nsigmaMCmat/Al/Pi", - "nsigmaMCmat/Al/Ka", "nsigmaMCmat/Al/Pr", "nsigmaMCmat/Al/De", - "nsigmaMCmat/Al/Tr", "nsigmaMCmat/Al/He", "nsigmaMCmat/Al/Al"}; - static constexpr std::string_view hnsigmaMCprm[NpNp] = {"nsigmaMCprm/El/El", "nsigmaMCprm/El/Mu", "nsigmaMCprm/El/Pi", - "nsigmaMCprm/El/Ka", "nsigmaMCprm/El/Pr", "nsigmaMCprm/El/De", - "nsigmaMCprm/El/Tr", "nsigmaMCprm/El/He", "nsigmaMCprm/El/Al", - "nsigmaMCprm/Mu/El", "nsigmaMCprm/Mu/Mu", "nsigmaMCprm/Mu/Pi", - "nsigmaMCprm/Mu/Ka", "nsigmaMCprm/Mu/Pr", "nsigmaMCprm/Mu/De", - "nsigmaMCprm/Mu/Tr", "nsigmaMCprm/Mu/He", "nsigmaMCprm/Mu/Al", - "nsigmaMCprm/Pi/El", "nsigmaMCprm/Pi/Mu", "nsigmaMCprm/Pi/Pi", - "nsigmaMCprm/Pi/Ka", "nsigmaMCprm/Pi/Pr", "nsigmaMCprm/Pi/De", - "nsigmaMCprm/Pi/Tr", "nsigmaMCprm/Pi/He", "nsigmaMCprm/Pi/Al", - "nsigmaMCprm/Ka/El", "nsigmaMCprm/Ka/Mu", "nsigmaMCprm/Ka/Pi", - "nsigmaMCprm/Ka/Ka", "nsigmaMCprm/Ka/Pr", "nsigmaMCprm/Ka/De", - "nsigmaMCprm/Ka/Tr", "nsigmaMCprm/Ka/He", "nsigmaMCprm/Ka/Al", - "nsigmaMCprm/Pr/El", "nsigmaMCprm/Pr/Mu", "nsigmaMCprm/Pr/Pi", - "nsigmaMCprm/Pr/Ka", "nsigmaMCprm/Pr/Pr", "nsigmaMCprm/Pr/De", - "nsigmaMCprm/Pr/Tr", "nsigmaMCprm/Pr/He", "nsigmaMCprm/Pr/Al", - "nsigmaMCprm/De/El", "nsigmaMCprm/De/Mu", "nsigmaMCprm/De/Pi", - "nsigmaMCprm/De/Ka", "nsigmaMCprm/De/Pr", "nsigmaMCprm/De/De", - "nsigmaMCprm/De/Tr", "nsigmaMCprm/De/He", "nsigmaMCprm/De/Al", - "nsigmaMCprm/Tr/El", "nsigmaMCprm/Tr/Mu", "nsigmaMCprm/Tr/Pi", - "nsigmaMCprm/Tr/Ka", "nsigmaMCprm/Tr/Pr", "nsigmaMCprm/Tr/De", - "nsigmaMCprm/Tr/Tr", "nsigmaMCprm/Tr/He", "nsigmaMCprm/Tr/Al", - "nsigmaMCprm/He/El", "nsigmaMCprm/He/Mu", "nsigmaMCprm/He/Pi", - "nsigmaMCprm/He/Ka", "nsigmaMCprm/He/Pr", "nsigmaMCprm/He/De", - "nsigmaMCprm/He/Tr", "nsigmaMCprm/He/He", "nsigmaMCprm/He/Al", - "nsigmaMCprm/Al/El", "nsigmaMCprm/Al/Mu", "nsigmaMCprm/Al/Pi", - "nsigmaMCprm/Al/Ka", "nsigmaMCprm/Al/Pr", "nsigmaMCprm/Al/De", - "nsigmaMCprm/Al/Tr", "nsigmaMCprm/Al/He", "nsigmaMCprm/Al/Al"}; - static constexpr const char* pT[Np] = {"e", "#mu", "#pi", "K", "p", "d", "t", "^{3}He", "#alpha"}; + static constexpr const char* pName[Np] = {"El", "Mu", "Pi", "Ka", "Pr", "De", "Tr", "He", "Al"}; static constexpr int PDGs[Np] = {11, 13, 211, 321, 2212, 1000010020, 1000010030, 1000020030}; HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -207,6 +110,7 @@ struct pidTofQaMc { Configurable doAl{"doAl", 0, "Process alpha"}; ConfigurableAxis binsPt{"binsPt", {2000, 0.f, 20.f}, "Binning of the pT axis"}; ConfigurableAxis binsNsigma{"binsNsigma", {2000, -50.f, 50.f}, "Binning of the NSigma axis"}; + ConfigurableAxis binsDelta{"binsDelta", {2000, -500.f, 500.f}, "Binning of the Delta axis"}; ConfigurableAxis binsSignal{"binsSignal", {6000, 0, 2000}, "Binning of the TPC signal axis"}; ConfigurableAxis binsLength{"binsLength", {1000, 0, 3000}, "Binning of the Length axis"}; ConfigurableAxis binsEta{"binsEta", {100, -4, 4}, "Binning of the Eta axis"}; @@ -271,33 +175,50 @@ struct pidTofQaMc { const AxisSpec lengthAxis{binsLength, "Track length (cm)"}; const AxisSpec etaAxis{binsEta, "#it{#eta}"}; const AxisSpec nSigmaAxis{binsNsigma, Form("N_{#sigma}^{TOF}(%s)", pT[mcID])}; + const AxisSpec deltaAxis{binsDelta, Form("#Delta^{TOF}(%s)", pT[mcID])}; + const AxisSpec deltaLengthAxis{binsDelta, Form("t_{exp}(%s)-t_{exp}^{*}(%s)", pT[mcID], pT[mcID])}; // Particle info - histos.add(hparticlept[mcID].data(), "", kTH1F, {pAxis}); - histos.add(hparticlep[mcID].data(), "", kTH1F, {pAxis}); - histos.add(hparticleeta[mcID].data(), "", kTH1F, {pAxis}); + hParticlePt[mcID] = histos.add(Form("particlept/%s", pName[mcID]), pT[mcID], kTH1D, {ptAxis}); + hParticleP[mcID] = histos.add(Form("particlep/%s", pName[mcID]), pT[mcID], kTH1D, {pAxis}); + hParticleEta[mcID] = histos.add(Form("particleeta/%s", pName[mcID]), pT[mcID], kTH1D, {etaAxis}); // Track info - histos.add(htrackpt[mcID].data(), "", kTH1F, {pAxis}); - histos.add(htrackp[mcID].data(), "", kTH1F, {pAxis}); - histos.add(htracketa[mcID].data(), "", kTH1F, {pAxis}); - histos.add(htracklength[mcID].data(), "", kTH1F, {pAxis}); + hTrackPt[mcID] = histos.add(Form("trackpt/%s", pName[mcID]), pT[mcID], kTH1D, {ptAxis}); + hTrackP[mcID] = histos.add(Form("trackp/%s", pName[mcID]), pT[mcID], kTH1D, {pAxis}); + hTrackEta[mcID] = histos.add(Form("tracketa/%s", pName[mcID]), pT[mcID], kTH1D, {etaAxis}); + hTrackLength[mcID] = histos.add(Form("tracklength/%s", pName[mcID]), pT[mcID], kTH1D, {lengthAxis}); // NSigma - histos.add(hnsigma[mcID].data(), pT[mcID], HistType::kTH2F, {ptAxis, nSigmaAxis}); - histos.add(hsignalMC[mcID].data(), pT[mcID], HistType::kTH2F, {pAxis, signalAxis}); + hSignalMC[mcID] = histos.add(Form("signalMC/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {pAxis, signalAxis}); + hNSigma[mcID] = histos.add(Form("nsigma/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, nSigmaAxis}); + hDeltaMCEvTime[mcID] = histos.add(Form("deltamcevtime/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + hDeltaMCEvTimeTrueGoodEv[mcID] = histos.add(Form("deltamcevtimegoodev/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + hDeltaMCEvTimeTrueBadEv[mcID] = histos.add(Form("deltamcevtimebadev/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + hDeltaMCEvTimeMC[mcID] = histos.add(Form("deltamcevtimeMC/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + hDeltaMCEvTimeMCGoodMatch[mcID] = histos.add(Form("deltamcevtimeMCgm/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + hDeltaMCEvTimeMCBadMatch[mcID] = histos.add(Form("deltamcevtimeMCbm/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + + hDeltaExpTimeMC[mcID] = histos.add(Form("deltaexptimeMC/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaLengthAxis}); if (!checkPrimaries) { return; } + hSignalMCprm[mcID] = histos.add(Form("signalMCprm/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {pAxis, signalAxis}); + hSignalMCstr[mcID] = histos.add(Form("signalMCstr/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {pAxis, signalAxis}); + hSignalMCmat[mcID] = histos.add(Form("signalMCmat/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {pAxis, signalAxis}); + + hNSigmaprm[mcID] = histos.add(Form("nsigmaprm/%s", pName[mcID]), Form("Primary %s", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); + hNSigmastr[mcID] = histos.add(Form("nsigmastr/%s", pName[mcID]), Form("Secondary %s from decay", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); + hNSigmamat[mcID] = histos.add(Form("nsigmamat/%s", pName[mcID]), Form("Secondary %s from material", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); - histos.add(hnsigmaprm[mcID].data(), Form("Primary %s", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); - histos.add(hnsigmastr[mcID].data(), Form("Secondary %s from decay", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); - histos.add(hnsigmamat[mcID].data(), Form("Secondary %s from material", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); + hDeltaMCEvTimeprm[mcID] = histos.add(Form("deltamcevtimeprm/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + hDeltaMCEvTimestr[mcID] = histos.add(Form("deltamcevtimestr/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + hDeltaMCEvTimemat[mcID] = histos.add(Form("deltamcevtimemat/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); - histos.add(hsignalMCprm[mcID].data(), Form("Primary %s", pT[mcID]), HistType::kTH2F, {pAxis, signalAxis}); - histos.add(hsignalMCstr[mcID].data(), Form("Secondary %s from decay", pT[mcID]), HistType::kTH2F, {pAxis, signalAxis}); - histos.add(hsignalMCmat[mcID].data(), Form("Secondary %s from material", pT[mcID]), HistType::kTH2F, {pAxis, signalAxis}); + hDeltaMCEvTimeMCprm[mcID] = histos.add(Form("deltamcevtimeMCprm/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + hDeltaMCEvTimeMCstr[mcID] = histos.add(Form("deltamcevtimeMCstr/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); + hDeltaMCEvTimeMCmat[mcID] = histos.add(Form("deltamcevtimeMCmat/%s", pName[mcID]), pT[mcID], HistType::kTH2F, {ptAxis, deltaAxis}); } template @@ -355,11 +276,11 @@ struct pidTofQaMc { const AxisSpec nSigmaAxis{binsNsigma, Form("N_{#sigma}^{TOF}(%s)", pT[massID])}; - histos.add(hnsigmaMC[mcID * Np + massID].data(), Form("True %s", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); + hNSigmaMC[mcID * Np + massID] = histos.add(Form("nsigmaMC/%s/%s", pName[mcID], pName[massID]), pT[mcID], HistType::kTH2F, {ptAxis, nSigmaAxis}); if (checkPrimaries) { - histos.add(hnsigmaMCprm[mcID * Np + massID].data(), Form("True Primary %s", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); - histos.add(hnsigmaMCstr[mcID * Np + massID].data(), Form("True Secondary %s from decay", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); - histos.add(hnsigmaMCmat[mcID * Np + massID].data(), Form("True Secondary %s from material", pT[mcID]), HistType::kTH2F, {ptAxis, nSigmaAxis}); + hNSigmaMCprm[mcID * Np + massID] = histos.add(Form("nsigmaMCprm/%s/%s", pName[mcID], pName[massID]), pT[mcID], HistType::kTH2F, {ptAxis, nSigmaAxis}); + hNSigmaMCstr[mcID * Np + massID] = histos.add(Form("nsigmaMCstr/%s/%s", pName[mcID], pName[massID]), pT[mcID], HistType::kTH2F, {ptAxis, nSigmaAxis}); + hNSigmaMCmat[mcID * Np + massID] = histos.add(Form("nsigmaMCmat/%s/%s", pName[mcID], pName[massID]), pT[mcID], HistType::kTH2F, {ptAxis, nSigmaAxis}); } } @@ -376,6 +297,8 @@ struct pidTofQaMc { histos.add("event/T0", ";Tracks with TOF;T0 (ps);Counts", HistType::kTH2F, {{1000, 0, 1000}, {1000, -1000, 1000}}); histos.add("event/vertexz", ";Vtx_{z} (cm);Entries", kTH1F, {{100, -20, 20}}); + histos.add("event/tofbadmatch", "", kTH1F, {ptAxis}); + histos.add("event/tofgoodmatch", "", kTH1F, {ptAxis}); static_for<0, 8>([&](auto i) { static_for<0, 8>([&](auto j) { @@ -401,7 +324,7 @@ struct pidTofQaMc { { switch (pdgSign.value) { case 0: - if (abs(particle.pdgCode()) != PDGs[mcID]) { + if (std::abs(particle.pdgCode()) != PDGs[mcID]) { return; } break; @@ -468,9 +391,9 @@ struct pidTofQaMc { LOG(fatal) << "Can't interpret index"; } - histos.fill(HIST(hparticlep[mcID]), particle.p()); - histos.fill(HIST(hparticlept[mcID]), particle.pt()); - histos.fill(HIST(hparticleeta[mcID]), particle.eta()); + hParticlePt[mcID]->Fill(particle.pt()); + hParticleP[mcID]->Fill(particle.p()); + hParticleEta[mcID]->Fill(particle.p()); } template @@ -529,23 +452,61 @@ struct pidTofQaMc { const float nsigma = o2::aod::pidutils::tofNSigma(track); // Fill for all - histos.fill(HIST(hnsigma[mcID]), track.pt(), nsigma); + hNSigma[mcID]->Fill(track.pt(), nsigma); + float expTime = 0.f; + switch (mcID) { + case 0: + expTime = track.tofExpTimeEl(); + break; + case 1: + expTime = track.tofExpTimeMu(); + break; + case 2: + expTime = track.tofExpTimePi(); + break; + case 3: + expTime = track.tofExpTimeKa(); + break; + case 4: + expTime = track.tofExpTimePr(); + break; + case 5: + expTime = track.tofExpTimeDe(); + break; + case 6: + expTime = track.tofExpTimeTr(); + break; + case 7: + expTime = track.tofExpTimeHe(); + break; + case 8: + expTime = track.tofExpTimeAl(); + break; + default: + break; + } + const float delta = track.tofSignal() - expTime - particle.mcCollision().t() * 1000.f; + + hDeltaMCEvTime[mcID]->Fill(track.pt(), delta); if (checkPrimaries) { if (!particle.isPhysicalPrimary()) { if (particle.getProcess() == 4) { - histos.fill(HIST(hnsigmastr[mcID]), track.pt(), nsigma); + hNSigmastr[mcID]->Fill(track.pt(), nsigma); + hDeltaMCEvTimestr[mcID]->Fill(track.pt(), delta); } else { - histos.fill(HIST(hnsigmamat[mcID]), track.pt(), nsigma); + hNSigmamat[mcID]->Fill(track.pt(), nsigma); + hDeltaMCEvTimemat[mcID]->Fill(track.pt(), delta); } } else { - histos.fill(HIST(hnsigmaprm[mcID]), track.pt(), nsigma); + hNSigmaprm[mcID]->Fill(track.pt(), nsigma); + hDeltaMCEvTimeprm[mcID]->Fill(track.pt(), delta); } } switch (pdgSign.value) { case 0: - if (abs(particle.pdgCode()) != PDGs[mcID]) { + if (std::abs(particle.pdgCode()) != PDGs[mcID]) { return; } break; @@ -564,22 +525,48 @@ struct pidTofQaMc { } // Track info - histos.fill(HIST(htrackp[mcID]), track.p()); - histos.fill(HIST(htrackpt[mcID]), track.pt()); - histos.fill(HIST(htracketa[mcID]), track.eta()); - histos.fill(HIST(htracklength[mcID]), track.length()); + hTrackPt[mcID]->Fill(track.pt()); + hTrackP[mcID]->Fill(track.p()); + hTrackEta[mcID]->Fill(track.eta()); + hTrackLength[mcID]->Fill(track.length()); // PID info - histos.fill(HIST(hsignalMC[mcID]), track.p(), track.beta()); + const float beta = track.beta(); + // const float beta = track.tofBeta(); + hSignalMC[mcID]->Fill(track.p(), beta); + hDeltaMCEvTimeMC[mcID]->Fill(track.pt(), delta); + + if (track.mcMask() & (0x1 << 15) && track.mcMask() & (0x1 << 13)) { + hDeltaMCEvTimeMCBadMatch[mcID]->Fill(track.pt(), delta); + } else { + hDeltaMCEvTimeMCGoodMatch[mcID]->Fill(track.pt(), delta); + } + + // Check that the track collision and the particle collisions match + if (particle.isPhysicalPrimary()) { + if (track.template collision_as().mcCollision().globalIndex() != particle.mcCollision().globalIndex()) { + hDeltaMCEvTimeTrueBadEv[mcID]->Fill(track.pt(), delta); + } else { + hDeltaMCEvTimeTrueGoodEv[mcID]->Fill(track.pt(), delta); + } + } + + const float mcExpTime = calculateTimeOfFlight(track.tpcInnerParam(), track.length()); + // const float mcExpTime = calculateTimeOfFlight(particle.p(), track.length()); + hDeltaExpTimeMC[mcID]->Fill(track.pt(), expTime - mcExpTime); + if (checkPrimaries) { if (!particle.isPhysicalPrimary()) { if (particle.getProcess() == 4) { - histos.fill(HIST(hsignalMCstr[mcID]), track.p(), track.beta()); + hSignalMCstr[mcID]->Fill(track.p(), beta); + hDeltaMCEvTimeMCstr[mcID]->Fill(track.pt(), delta); } else { - histos.fill(HIST(hsignalMCmat[mcID]), track.p(), track.beta()); + hSignalMCmat[mcID]->Fill(track.p(), beta); + hDeltaMCEvTimeMCmat[mcID]->Fill(track.pt(), delta); } } else { - histos.fill(HIST(hsignalMCprm[mcID]), track.p(), track.beta()); + hSignalMCprm[mcID]->Fill(track.p(), beta); + hDeltaMCEvTimeMCprm[mcID]->Fill(track.pt(), delta); } } } @@ -639,7 +626,7 @@ struct pidTofQaMc { switch (pdgSign.value) { case 0: - if (abs(particle.pdgCode()) != PDGs[mcID]) { + if (std::abs(particle.pdgCode()) != PDGs[mcID]) { return; } break; @@ -659,29 +646,31 @@ struct pidTofQaMc { const float nsigmaMassID = o2::aod::pidutils::tofNSigma(track); - histos.fill(HIST(hnsigmaMC[mcID * Np + massID]), track.pt(), nsigmaMassID); + hNSigmaMC[mcID * Np + massID]->Fill(track.pt(), nsigmaMassID); if (checkPrimaries) { if (!particle.isPhysicalPrimary()) { if (particle.getProcess() == 4) { - histos.fill(HIST(hnsigmaMCstr[mcID * Np + massID]), track.pt(), nsigmaMassID); + hNSigmaMCstr[mcID * Np + massID]->Fill(track.pt(), nsigmaMassID); } else { - histos.fill(HIST(hnsigmaMCmat[mcID * Np + massID]), track.pt(), nsigmaMassID); + hNSigmaMCmat[mcID * Np + massID]->Fill(track.pt(), nsigmaMassID); } } else { - histos.fill(HIST(hnsigmaMCprm[mcID * Np + massID]), track.pt(), nsigmaMassID); + hNSigmaMCprm[mcID * Np + massID]->Fill(track.pt(), nsigmaMassID); } } } - Preslice perCol = aod::track::collisionId; - Preslice perMCCol = aod::mcparticle::mcCollisionId; - - void process(soa::Join const& collisions, - soa::Join& tracks, + aod::TOFSignal, aod::McTrackLabels, aod::pidTOFbeta>; + using Colls = soa::Join; + Preslice perCol = aod::track::collisionId; + Preslice perMCCol = aod::mcparticle::mcCollisionId; + + void process(Colls const& collisions, + Trks& tracks, aod::McParticles& mcParticles, aod::McCollisions&) { @@ -689,10 +678,12 @@ struct pidTofQaMc { if (collision.numContrib() < nMinNumberOfContributors) { return; } + if (!collision.sel8()) { + continue; + } if (!collision.has_mcCollision()) { continue; } - const auto tracksInCollision = tracks.sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); const auto particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, collision.mcCollision().globalIndex(), cache); for (const auto& p : particlesInCollision) { @@ -701,52 +692,66 @@ struct pidTofQaMc { }); } + const auto& tracksInCollision = tracks.sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + // tracksInCollision.bindExternalIndices(&mcParticles); + // const auto& tracksWithPid = soa::Attach(tracksInCollision); + // tracksInCollision.copyIndexBindings(tracksWithPid); const float collisionTime_ps = collision.collisionTime() * 1000.f; unsigned int nTracksWithTOF = 0; - for (const auto& t : tracksInCollision) { + for (const auto& track : tracksInCollision) { // - if (!t.hasTOF()) { // Skipping tracks without TOF + if (!track.hasTOF()) { // Skipping tracks without TOF continue; } - if (t.eta() < minEta || t.eta() > maxEta) { + if (track.eta() < minEta || track.eta() > maxEta) { continue; } nTracksWithTOF++; // Fill for all - histos.fill(HIST("event/tofbeta"), t.p(), t.beta()); - if (!t.has_mcParticle()) { + // const float beta = track.tofBeta(); + const float beta = track.beta(); + histos.fill(HIST("event/tofbeta"), track.p(), beta); + if (!track.has_mcParticle()) { continue; } - const auto& particle = t.mcParticle(); + std::bitset<16> bits(track.mcMask()); + // LOG(info) << "Using bitset: " << bits; + if (bits[15]) { + histos.fill(HIST("event/tofbadmatch"), track.pt()); + } else { + histos.fill(HIST("event/tofgoodmatch"), track.pt()); + } + + const auto& particle = track.mcParticle(); if (checkPrimaries) { if (!particle.isPhysicalPrimary()) { - // histos.fill(HIST("event/tofbetaSec"), t.p(), t.beta()); + // histos.fill(HIST("event/tofbetaSec"), track.p(), beta); if (particle.getProcess() == 4) { - histos.fill(HIST("event/tofbetaStr"), t.tpcInnerParam(), t.tpcSignal()); + histos.fill(HIST("event/tofbetaStr"), track.tpcInnerParam(), track.tpcSignal()); } else { - histos.fill(HIST("event/tofbetaMat"), t.tpcInnerParam(), t.tpcSignal()); + histos.fill(HIST("event/tofbetaMat"), track.tpcInnerParam(), track.tpcSignal()); } } else { - histos.fill(HIST("event/tofbetaPrm"), t.p(), t.beta()); + histos.fill(HIST("event/tofbetaPrm"), track.p(), beta); } } // Fill with PDG codes static_for<0, 8>([&](auto i) { static_for<0, 8>([&](auto j) { - fillPIDInfoForPdg(t, particle); + fillPIDInfoForPdg(track, particle); }); - fillTrackInfoForPdg(t, particle); + fillTrackInfoForPdg(track, particle); }); } // track loop histos.fill(HIST("event/T0"), nTracksWithTOF, collisionTime_ps); histos.fill(HIST("event/vertexz"), collision.posZ()); } // collision loop - } // process() + } // process() }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; } From bdbef20d954e4407c5fcd744f55510cf17e88d66 Mon Sep 17 00:00:00 2001 From: Nicolas Strangmann <77485327+nstrangm@users.noreply.github.com> Date: Fri, 13 Dec 2024 11:32:42 +0100 Subject: [PATCH 2/2] [PWGEM/PhotonMeson] Fix abs in emcalQC task (#8957) Co-authored-by: Nicolas Strangmann --- PWGEM/PhotonMeson/Tasks/emcalQC.cxx | 97 +++++++++++++++-------------- 1 file changed, 50 insertions(+), 47 deletions(-) diff --git a/PWGEM/PhotonMeson/Tasks/emcalQC.cxx b/PWGEM/PhotonMeson/Tasks/emcalQC.cxx index 0b5be17d212..8ed00c1f50e 100644 --- a/PWGEM/PhotonMeson/Tasks/emcalQC.cxx +++ b/PWGEM/PhotonMeson/Tasks/emcalQC.cxx @@ -9,10 +9,14 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -// ======================== -// -// This code runs loop over EMCal clusters for EMCal QC. -// Please write to: nicolas.strangmann@cern.ch +/// EMCAL QC Task +/// +/// \file emcalQC.cxx +/// +/// \brief Task that runs basic EMCal cluster QA for derived data in the EM format +/// +/// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) Goethe University Frankfurt +/// #include #include @@ -50,10 +54,10 @@ using MyCollision = MyCollisions::iterator; using MyEMCClusters = soa::Join; using MyEMCCluster = MyEMCClusters::iterator; -struct emcalQC { +struct EmcalQC { Configurable cfgDo2DQA{"cfgDo2DQA", true, "perform 2 dimensional cluster QA"}; - ConfigurableAxis ConfVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; + ConfigurableAxis confVtxBins{"confVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; EMPhotonEventCut fEMEventCut; struct : ConfigurableGroup { @@ -78,19 +82,19 @@ struct emcalQC { std::string prefix = "emccut_group"; Configurable clusterDefinition{"clusterDefinition", "kV3Default", "Clusterizer to be selected, e.g. V3Default"}; Configurable minOpenAngle{"minOpenAngle", 0.0202, "apply min opening angle"}; - Configurable EMC_minTime{"EMC_minTime", -20., "Minimum cluster time for EMCal time cut"}; - Configurable EMC_maxTime{"EMC_maxTime", +25., "Maximum cluster time for EMCal time cut"}; - Configurable EMC_minM02{"EMC_minM02", 0.1, "Minimum M02 for EMCal M02 cut"}; - Configurable EMC_maxM02{"EMC_maxM02", 0.7, "Maximum M02 for EMCal M02 cut"}; - Configurable EMC_minE{"EMC_minE", 0.7, "Minimum cluster energy for EMCal energy cut"}; - Configurable EMC_minNCell{"EMC_minNCell", 1, "Minimum number of cells per cluster for EMCal NCell cut"}; - Configurable> EMC_TM_Eta{"EMC_TM_Eta", {0.01f, 4.07f, -2.5f}, "|eta| <= [0]+(pT+[1])^[2] for EMCal track matching"}; - Configurable> EMC_TM_Phi{"EMC_TM_Phi", {0.015f, 3.65f, -2.f}, "|phi| <= [0]+(pT+[1])^[2] for EMCal track matching"}; - Configurable EMC_Eoverp{"EMC_Eoverp", 1.75, "Minimum cluster energy over track momentum for EMCal track matching"}; - Configurable EMC_UseExoticCut{"EMC_UseExoticCut", true, "FLag to use the EMCal exotic cluster cut"}; + Configurable minClusterTime{"minClusterTime", -20., "Minimum cluster time for EMCal time cut"}; + Configurable maxClusterTime{"maxClusterTime", +25., "Maximum cluster time for EMCal time cut"}; + Configurable minM02{"minM02", 0.1, "Minimum M02 for EMCal M02 cut"}; + Configurable maxM02{"maxM02", 0.7, "Maximum M02 for EMCal M02 cut"}; + Configurable minClusterE{"minClusterE", 0.7, "Minimum cluster energy for EMCal energy cut"}; + Configurable minNCell{"minNCell", 1, "Minimum number of cells per cluster for EMCal NCell cut"}; + Configurable> tmEta{"tmEta", {0.01f, 4.07f, -2.5f}, "|eta| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable> tmPhi{"tmPhi", {0.015f, 3.65f, -2.f}, "|phi| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable tmEoverP{"tmEoverP", 1.75, "Minimum cluster energy over track momentum for EMCal track matching"}; + Configurable useExoticCut{"useExoticCut", true, "FLag to use the EMCal exotic cluster cut"}; } emccuts; - void DefineEMEventCut() + void defineEMEventCut() { fEMEventCut = EMPhotonEventCut("fEMEventCut", "fEMEventCut"); fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); @@ -105,40 +109,40 @@ struct emcalQC { fEMEventCut.SetRequireEMCHardwareTriggered(eventcuts.cfgRequireEMCHardwareTriggered); } - void DefineEMCCut() + void defineEMCCut() { - const float a = emccuts.EMC_TM_Eta->at(0); - const float b = emccuts.EMC_TM_Eta->at(1); - const float c = emccuts.EMC_TM_Eta->at(2); + const float a = emccuts.tmEta->at(0); + const float b = emccuts.tmEta->at(1); + const float c = emccuts.tmEta->at(2); - const float d = emccuts.EMC_TM_Phi->at(0); - const float e = emccuts.EMC_TM_Phi->at(1); - const float f = emccuts.EMC_TM_Phi->at(2); + const float d = emccuts.tmPhi->at(0); + const float e = emccuts.tmPhi->at(1); + const float f = emccuts.tmPhi->at(2); LOGF(info, "EMCal track matching parameters : a = %f, b = %f, c = %f, d = %f, e = %f, f = %f", a, b, c, d, e, f); fEMCCut.SetClusterizer(emccuts.clusterDefinition); - fEMCCut.SetMinE(emccuts.EMC_minE); - fEMCCut.SetMinNCell(emccuts.EMC_minNCell); - fEMCCut.SetM02Range(emccuts.EMC_minM02, emccuts.EMC_maxM02); - fEMCCut.SetTimeRange(emccuts.EMC_minTime, emccuts.EMC_maxTime); + fEMCCut.SetMinE(emccuts.minClusterE); + fEMCCut.SetMinNCell(emccuts.minNCell); + fEMCCut.SetM02Range(emccuts.minM02, emccuts.maxM02); + fEMCCut.SetTimeRange(emccuts.minClusterTime, emccuts.maxClusterTime); - fEMCCut.SetTrackMatchingEta([a, b, c](float pT) { return a + pow(pT + b, c); }); - fEMCCut.SetTrackMatchingPhi([d, e, f](float pT) { return d + pow(pT + e, f); }); + fEMCCut.SetTrackMatchingEta([a, b, c](float pT) { return a + std::pow(pT + b, c); }); + fEMCCut.SetTrackMatchingPhi([d, e, f](float pT) { return d + std::pow(pT + e, f); }); - fEMCCut.SetMinEoverP(emccuts.EMC_Eoverp); - fEMCCut.SetUseExoticCut(emccuts.EMC_UseExoticCut); + fEMCCut.SetMinEoverP(emccuts.tmEoverP); + fEMCCut.SetUseExoticCut(emccuts.useExoticCut); } HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; - std::vector zvtx_bin_edges; + std::vector zVtxBinEdges; void init(InitContext&) { - zvtx_bin_edges = std::vector(ConfVtxBins.value.begin(), ConfVtxBins.value.end()); - zvtx_bin_edges.erase(zvtx_bin_edges.begin()); + zVtxBinEdges = std::vector(confVtxBins.value.begin(), confVtxBins.value.end()); + zVtxBinEdges.erase(zVtxBinEdges.begin()); - DefineEMCCut(); - DefineEMEventCut(); + defineEMCCut(); + defineEMEventCut(); o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(&fRegistry); auto hEMCCollisionCounter = fRegistry.add("Event/hEMCCollisionCounter", "Number of collisions after event cuts", HistType::kTH1F, {{7, 0.5, 7.5}}, false); @@ -156,16 +160,16 @@ struct emcalQC { void processQC(MyCollisions const& collisions, MyEMCClusters const& clusters) { - for (auto& collision : collisions) { + for (const auto& collision : collisions) { - if (eventcuts.onlyKeepWeightedEvents && fabs(collision.weight() - 1.) < 1E-10) { + if (eventcuts.onlyKeepWeightedEvents && std::fabs(collision.weight() - 1.) < 1E-10) { continue; } fRegistry.fill(HIST("Event/hEMCCollisionCounter"), 1); if (collision.selection_bit(o2::aod::evsel::kIsTriggerTVX)) { fRegistry.fill(HIST("Event/hEMCCollisionCounter"), 2); - if (abs(collision.posZ()) < eventcuts.cfgZvtxMax) { + if (std::abs(collision.posZ()) < eventcuts.cfgZvtxMax) { fRegistry.fill(HIST("Event/hEMCCollisionCounter"), 3); if (collision.sel8()) { fRegistry.fill(HIST("Event/hEMCCollisionCounter"), 4); @@ -193,11 +197,11 @@ struct emcalQC { fRegistry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted fRegistry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted - auto clusters_per_coll = clusters.sliceBy(perCollision, collision.collisionId()); - fRegistry.fill(HIST("Cluster/before/hNgamma"), clusters_per_coll.size(), collision.weight()); + auto clustersPerColl = clusters.sliceBy(perCollision, collision.collisionId()); + fRegistry.fill(HIST("Cluster/before/hNgamma"), clustersPerColl.size(), collision.weight()); int ngBefore = 0; int ngAfter = 0; - for (auto& cluster : clusters_per_coll) { + for (const auto& cluster : clustersPerColl) { // Fill the cluster properties before applying any cuts if (!fEMCCut.IsSelectedEMCal(EMCPhotonCut::EMCPhotonCuts::kDefinition, cluster)) continue; @@ -236,12 +240,11 @@ struct emcalQC { void processDummy(MyCollisions const&) {} - PROCESS_SWITCH(emcalQC, processQC, "run EMCal QC", false); - PROCESS_SWITCH(emcalQC, processDummy, "Dummy function", true); + PROCESS_SWITCH(EmcalQC, processQC, "run EMCal QC", false); + PROCESS_SWITCH(EmcalQC, processDummy, "Dummy function", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{ - adaptAnalysisTask(cfgc, TaskName{"emcal-qc"})}; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; }