diff --git a/DPG/Tasks/AOTTrack/PID/TPC/CMakeLists.txt b/DPG/Tasks/AOTTrack/PID/TPC/CMakeLists.txt index 01d69f64386..f804cea4c51 100644 --- a/DPG/Tasks/AOTTrack/PID/TPC/CMakeLists.txt +++ b/DPG/Tasks/AOTTrack/PID/TPC/CMakeLists.txt @@ -24,3 +24,8 @@ o2physics_add_dpl_workflow(pid-tpc-qa-mc SOURCES qaPIDTPCMC.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(tree-creator-pid-tpc-qa + SOURCES treeCreatorPidTpcQa.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB + COMPONENT_NAME Analysis) diff --git a/DPG/Tasks/AOTTrack/PID/TPC/treeCreatorPidTpcQa.cxx b/DPG/Tasks/AOTTrack/PID/TPC/treeCreatorPidTpcQa.cxx new file mode 100644 index 00000000000..130e275ec17 --- /dev/null +++ b/DPG/Tasks/AOTTrack/PID/TPC/treeCreatorPidTpcQa.cxx @@ -0,0 +1,297 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file treeCreatorPidTpcQa.cxx +/// \brief Creates trees with PID QA variables along with variables used for NN training +/// +/// \author Ana Marin +/// \author Oleksii Lubynets + +#include "treeCreatorPidTpcQa.h" + +#include "Common/CCDB/RCTSelectionFlags.h" +#include "Common/CCDB/ctpRateFetcher.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "DPG/Tasks/TPC/tpcSkimsTableCreator.h" +#include "DPG/Tasks/TPC/utilsTpcSkimsTableCreator.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::track; +using namespace o2::dpg_tpcskimstablecreator; + +#define DO_FOR_ALL_PARTICLES(MACRO) \ + MACRO(El, Electron) \ + MACRO(Mu, Muon) \ + MACRO(Pi, Pion) \ + MACRO(Ka, Kaon) \ + MACRO(Pr, Proton) \ + MACRO(De, Deuteron) \ + MACRO(Tr, Triton) \ + MACRO(He, Helium3) \ + MACRO(Al, Alpha) + +struct TreeCreatorPidTpcQa { + Produces rowPidTpcQa; + + Configurable applyEvSel{"applyEvSel", 2, "Flag to apply event selection: 0 -> no event selection, 1 -> Run 2 event selection, 2 -> Run 3 event selection"}; + Configurable cutVtxZ{"cutVtxZ", 10.f, "Cut on vertex Z position [cm]"}; + Configurable trackSelection{"trackSelection", 1, "Track selection: 0 -> No Cut, 1 -> kGlobalTrack, 2 -> kGlobalTrackWoPtEta, 3 -> kGlobalTrackWoDCA, 4 -> kQualityTracks, 5 -> kInAcceptanceTracks"}; + Configurable requireGlobalTrack{"requireGlobalTrack", true, "Skip non-global tracks"}; + Configurable requireIts{"requireIts", true, "Skip tracks without ITS"}; + Configurable cutMinTPCNcls{"cutMinTPCNcls", 0, "Minimum number or TPC Clusters for tracks"}; + Configurable cutRapidity{"cutRapidity", 999.f, "Rapidity cut"}; + Configurable nClNorm{"nClNorm", 152.f, "Number of cluster normalization. Run 2: 159, Run 3 152"}; + // Configurable for the path of CCDB General Run Parameters LHC Interface information + Configurable ccdbPathGrpLhcIf{"ccdbPathGrpLhcIf", "GLO/Config/GRPLHCIF", "Path on the CCDB for the GRPLHCIF object"}; + // Configurables for output tables reservation size + Configurable reserveRatio{"reserveRatio", 1.f, "Ratio of how many rows expected in the output table to the input Tracks table size"}; + Configurable saveReserveQaHisto{"saveReserveQaHisto", true, "Flag to save the DF-wise ratio of output table size to that of input table"}; + // Configurables for run condtion table + Configurable rctLabel{"rctLabel", "CBT_hadronPID", "select 1 [CBT, CBT_hadronPID, CBT_muon_glo] see O2Physics/Common/CCDB/RCTSelectionFlags.h"}; + Configurable checkZdc{"checkZdc", false, "set ZDC flag for PbPb"}; + Configurable treatLimitedAcceptanceAsBad{"treatLimitedAcceptanceAsBad", false, "reject all events where the detectors relevant for the specified Runlist are flagged as LimitedAcceptance"}; + Configurable requireGoodRct{"requireGoodRct", false, "require good detector flag in run condtion table"}; + +#define DECLARE_PARTICLE_WISE_CONFIGURABLES(ParticleNameShort, ParticleNameLong) \ + Configurable cutTpcInnerParameterMin##ParticleNameLong{"cutTpcInnerParameterMin" #ParticleNameLong, 0.f, "Lower-value cut on tpcInnerParam for " #ParticleNameLong}; /* o2-linter: disable=name/configurable (Configurable defined in macro)*/ \ + Configurable cutTpcInnerParameterMax##ParticleNameLong{"cutTpcInnerParameterMax" #ParticleNameLong, 999.f, "Upper-value cut on tpcInnerParam for " #ParticleNameLong}; /* o2-linter: disable=name/configurable (Configurable defined in macro)*/ \ + Configurable cutNSigmaTpcAbs##ParticleNameLong{"cutNSigmaTpcAbs" #ParticleNameLong, 999.f, "Cut on absolute value of nSigmaTpc for " #ParticleNameLong}; // o2-linter: disable=name/configurable (Configurable defined in macro) + + DO_FOR_ALL_PARTICLES(DECLARE_PARTICLE_WISE_CONFIGURABLES) +#undef DECLARE_PARTICLE_WISE_CONFIGURABLES + +#define PACK_CONFIGURABLES_TO_ARRAY(ParticleNameShort, ParticleNameLong) &cutTpcInnerParameterMin##ParticleNameLong, + std::array*, PID::Alpha + 1> cutTpcInnerParameterMin{ + DO_FOR_ALL_PARTICLES(PACK_CONFIGURABLES_TO_ARRAY)}; +#undef PACK_CONFIGURABLES_TO_ARRAY + +#define PACK_CONFIGURABLES_TO_ARRAY(ParticleNameShort, ParticleNameLong) &cutTpcInnerParameterMax##ParticleNameLong, + std::array*, PID::Alpha + 1> cutTpcInnerParameterMax{ + DO_FOR_ALL_PARTICLES(PACK_CONFIGURABLES_TO_ARRAY)}; +#undef PACK_CONFIGURABLES_TO_ARRAY + +#define PACK_CONFIGURABLES_TO_ARRAY(ParticleNameShort, ParticleNameLong) &cutNSigmaTpcAbs##ParticleNameLong, + std::array*, PID::Alpha + 1> cutNSigmaTpcAbs{ + DO_FOR_ALL_PARTICLES(PACK_CONFIGURABLES_TO_ARRAY)}; +#undef PACK_CONFIGURABLES_TO_ARRAY + + HistogramRegistry registry{"registry", {}}; + + Service ccdb{}; + + ctpRateFetcher mRateFetcher{}; + + o2::aod::rctsel::RCTFlagsChecker rctChecker{}; + + using CollisionsExtra = soa::Join; + using TrackCandidates = soa::Join; + + Preslice perCollisionTracks = aod::track::collisionId; + + int mEnabledParticles{0}; + int mProcessedParticles{0}; + + template + bool initPerParticle() + { + static_assert(ParticleId >= 0 && ParticleId <= PID::Alpha && "Particle index outside limits"); + int enabledProcesses{0}; + + switch (ParticleId) { +#define INIT_PARTICLE(ParticleNameShort, ParticleNameLong) \ + case PID::ParticleNameLong: \ + if (!doprocess##ParticleNameLong && !doprocessFull##ParticleNameLong && !doprocessFullWithTOF##ParticleNameLong) { \ + return false; \ + } \ + if (doprocess##ParticleNameLong) { \ + ++enabledProcesses; \ + } \ + if (doprocessFull##ParticleNameLong) { \ + ++enabledProcesses; \ + } \ + if (doprocessFullWithTOF##ParticleNameLong) { \ + ++enabledProcesses; \ + } \ + LOG(info) << "Enabled TPC QA for " << #ParticleNameLong; \ + break; + + DO_FOR_ALL_PARTICLES(INIT_PARTICLE) +#undef INIT_PARTICLE + } + if (enabledProcesses != 1) { + LOG(fatal) << "Cannot enable more than one process function per particle, check and retry!"; + } + return true; + } + + void init(o2::framework::InitContext&) + { + static_for<0, PID::Alpha>([&](auto ParticleId) { + mEnabledParticles += static_cast(initPerParticle()); + }); + + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setFatalWhenNull(false); + + rctChecker.init(rctLabel, checkZdc, treatLimitedAcceptanceAsBad); + + if (saveReserveQaHisto) { + registry.add("hOutputRatio", "Table out/in ratio;Table out/in ratio;Entries", {HistType::kTH1F, {{100, 0, reserveRatio}}}); + } + } + + template + void processSingleParticle(CollisionsExtra const& collisions, + TrackType const& tracks) + { + if (mProcessedParticles == 0) { + rowPidTpcQa.reserve(tracks.size() * reserveRatio); + } + + std::string irSource{}; + float sqrtSNN{}; // placeholder to satisfy evaluateIrSourceAndSqrtSnn's signature + bool isFirstCollision{true}; + for (const auto& collision : collisions) { + if (!isEventSelected(collision, applyEvSel) || (std::abs(collision.posZ()) > cutVtxZ)) { + continue; + } + + const bool isGoodRctEvent = rctChecker.checkTable(collision); + if (requireGoodRct && !isGoodRctEvent) { + continue; + } + + const auto bc = collision.bc_as(); + if (isFirstCollision) { + evaluateIrSourceAndSqrtSnn(ccdb, ccdbPathGrpLhcIf, bc.timestamp(), irSource, sqrtSNN); + } + isFirstCollision = false; + const float ft0Occ = collision.ft0cOccupancyInTimeRange(); + const float multTPC = collision.multTPC() / MultiplicityNorm; + const auto hadronicRate = !irSource.empty() ? mRateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), irSource) * OneToKilo : 0.; + + const auto tracksFromCollision = tracks.sliceBy(perCollisionTracks, static_cast(collision.globalIndex())); + + for (const auto& track : tracksFromCollision) { + bool isGoodTrack = isTrackSelected(track, trackSelection); + isGoodTrack &= (!requireGlobalTrack || track.isGlobalTrack()); + isGoodTrack &= (!requireIts || track.hasITS()); + isGoodTrack &= track.hasTPC(); + isGoodTrack &= (track.tpcNClsFound() >= cutMinTPCNcls); + + const float rapidity = track.rapidity(PID::getMass(ParticleId)); + const float momentum = track.p(); + const float nClNormalized = std::sqrt(nClNorm / track.tpcNClsFound()); + const auto nclPID = static_cast(track.tpcNClsPID()); + const float phi = track.phi(); + const float tgl = track.tgl(); + const float tpcInnerParam = track.tpcInnerParam(); + const float signed1Pt = track.signed1Pt(); + const float nSigmaTpc = o2::aod::pidutils::tpcNSigma(track); + + isGoodTrack &= (std::fabs(rapidity) <= cutRapidity); + isGoodTrack &= (tpcInnerParam >= *cutTpcInnerParameterMin.at(ParticleId)); + isGoodTrack &= (tpcInnerParam <= *cutTpcInnerParameterMax.at(ParticleId)); + isGoodTrack &= (std::fabs(nSigmaTpc) <= *cutNSigmaTpcAbs.at(ParticleId)); + + if (!isGoodTrack) { + continue; + } + + float dedxDiff{UndefValueFloat}; + float dedxExpected{UndefValueFloat}; + float expSigma{UndefValueFloat}; + + if constexpr (IsFullTable) { + dedxDiff = o2::aod::pidutils::tpcExpSignalDiff(track); + dedxExpected = track.tpcSignal() - dedxDiff; + expSigma = o2::aod::pidutils::tpcExpSigma(track); + } + + float nSigmaTof{UndefValueFloat}; + + if constexpr (IsTofTable) { + nSigmaTof = o2::aod::pidutils::tofNSigma(track); + } + + rowPidTpcQa(isGoodRctEvent, ParticleId, ft0Occ, hadronicRate, multTPC, nClNormalized, nclPID, phi, tgl, tpcInnerParam, rapidity, momentum, signed1Pt, nSigmaTpc, dedxExpected, dedxDiff, expSigma, nSigmaTof); + } // tracksFromCollision + } // collisions + ++mProcessedParticles; + if (mProcessedParticles == mEnabledParticles) { + mProcessedParticles = 0; + if (saveReserveQaHisto) { + registry.fill(HIST("hOutputRatio"), static_cast((rowPidTpcQa.lastIndex() + 1)) / tracks.size()); + } + } + } + +#define MAKE_PROCESS_FUNCTIONS(ParticleNameShort, ParticleNameLong) \ + void process##ParticleNameLong(CollisionsExtra const& collisions, \ + soa::Join const& tracks, \ + aod::BCsWithTimestamps const&) \ + { \ + processSingleParticle(collisions, tracks); \ + } \ + PROCESS_SWITCH(TreeCreatorPidTpcQa, process##ParticleNameLong, Form("Process for the %s hypothesis for TPC NSigma QA", #ParticleNameLong), false); \ + \ + void processFull##ParticleNameLong(CollisionsExtra const& collisions, \ + soa::Join const& tracks, \ + aod::BCsWithTimestamps const&) \ + { \ + processSingleParticle(collisions, tracks); \ + } \ + PROCESS_SWITCH(TreeCreatorPidTpcQa, processFull##ParticleNameLong, Form("Process for the %s hypothesis for full TPC PID QA", #ParticleNameLong), false); \ + \ + void processFullWithTOF##ParticleNameLong(CollisionsExtra const& collisions, \ + soa::Join const& tracks, \ + aod::BCsWithTimestamps const&) \ + { \ + processSingleParticle(collisions, tracks); \ + } \ + PROCESS_SWITCH(TreeCreatorPidTpcQa, processFullWithTOF##ParticleNameLong, Form("Process for the %s hypothesis for full TPC PID QA with the TOF info added", #ParticleNameLong), false); + + DO_FOR_ALL_PARTICLES(MAKE_PROCESS_FUNCTIONS) +#undef MAKE_PROCESS_FUNCTIONS +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} +#undef DO_FOR_ALL_PARTICLES diff --git a/DPG/Tasks/AOTTrack/PID/TPC/treeCreatorPidTpcQa.h b/DPG/Tasks/AOTTrack/PID/TPC/treeCreatorPidTpcQa.h new file mode 100644 index 00000000000..f34a4c52a8f --- /dev/null +++ b/DPG/Tasks/AOTTrack/PID/TPC/treeCreatorPidTpcQa.h @@ -0,0 +1,57 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file treeCreatorPidTpcQa.h +/// \author Ana Marin +/// \author Oleksii Lubynets +/// \brief Creates trees with PID QA variables along with variables used for NN training + +#ifndef DPG_TASKS_AOTTRACK_PID_TPC_TREECREATORPIDTPCQA_H_ +#define DPG_TASKS_AOTTRACK_PID_TPC_TREECREATORPIDTPCQA_H_ + +#include "DPG/Tasks/TPC/tpcSkimsTableCreator.h" + +#include +#include + +namespace o2::aod +{ +namespace dpg_tpcpidqa +{ +DECLARE_SOA_COLUMN(NSigmaTpc, nSigmaTpc, float); +DECLARE_SOA_COLUMN(DedxExpected, dedxExpected, float); +DECLARE_SOA_COLUMN(DedxDiff, dedxDiff, float); +DECLARE_SOA_COLUMN(ExpSigma, expSigma, float); +DECLARE_SOA_COLUMN(NSigmaTof, nSigmaTof, float); +} // namespace dpg_tpcpidqa + +DECLARE_SOA_TABLE(QaPidTpc, "AOD", "QAPIDTPC", + tpcskims::IsGoodRct, + tpcskims::PidIndex, + tpcskims::Ft0Occ, + tpcskims::HadronicRate, + tpcskims::NormMultTPC, + tpcskims::NormNClustersTPC, + tpcskims::NormNClustersTPCPID, + o2::aod::track::Phi, + o2::aod::track::Tgl, + o2::aod::track::TPCInnerParam, + o2::aod::track::Y, + o2::aod::track::P, + o2::aod::track::Signed1Pt, + dpg_tpcpidqa::NSigmaTpc, + dpg_tpcpidqa::DedxExpected, + dpg_tpcpidqa::DedxDiff, + dpg_tpcpidqa::ExpSigma, + dpg_tpcpidqa::NSigmaTof) +} // namespace o2::aod + +#endif // DPG_TASKS_AOTTRACK_PID_TPC_TREECREATORPIDTPCQA_H_