From 674113f4987c61d165ce994bae95758f09896036 Mon Sep 17 00:00:00 2001 From: Ana Marin Date: Tue, 24 Jun 2025 14:30:09 +0200 Subject: [PATCH 1/4] [DPG][Common] Test of dE/dx correction for occupancy --- Common/TableProducer/PID/CMakeLists.txt | 2 +- Common/TableProducer/PID/pidTPC.cxx | 124 +++++- Common/TableProducer/PID/pidTPCBase.cxx | 114 ++++- Common/TableProducer/PID/pidTPCBase.h | 46 ++ DPG/Tasks/TPC/tpcSkimsTableCreator.cxx | 555 +++++++++++++++++++++++- 5 files changed, 796 insertions(+), 45 deletions(-) diff --git a/Common/TableProducer/PID/CMakeLists.txt b/Common/TableProducer/PID/CMakeLists.txt index d28a3268954..afcda3e1be6 100644 --- a/Common/TableProducer/PID/CMakeLists.txt +++ b/Common/TableProducer/PID/CMakeLists.txt @@ -45,7 +45,7 @@ o2physics_add_dpl_workflow(pid-tof-full o2physics_add_dpl_workflow(pid-tpc-base SOURCES pidTPCBase.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(pid-tpc diff --git a/Common/TableProducer/PID/pidTPC.cxx b/Common/TableProducer/PID/pidTPC.cxx index 64500a47cb7..6cdca389a65 100644 --- a/Common/TableProducer/PID/pidTPC.cxx +++ b/Common/TableProducer/PID/pidTPC.cxx @@ -18,10 +18,10 @@ /// \brief Task to produce PID tables for TPC split for each particle. /// Only the tables for the mass hypotheses requested are filled, and only for the requested table size ("Full" or "Tiny"). The others are sent empty. /// -#include #include #include #include +#include #include // ROOT includes #include "TFile.h" @@ -29,21 +29,23 @@ #include "TSystem.h" // O2 includes +#include "MetadataHelper.h" +#include "TableHelper.h" +#include "pidTPCBase.h" + +#include "Common/Core/PID/TPCPIDResponse.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Tools/ML/model.h" + #include "CCDB/BasicCCDBManager.h" +#include "CCDB/CcdbApi.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" -#include "Framework/ASoAHelpers.h" #include "ReconstructionDataFormats/Track.h" -#include "CCDB/CcdbApi.h" -#include "Common/DataModel/PIDResponseTPC.h" -#include "Common/Core/PID/TPCPIDResponse.h" -#include "Framework/AnalysisDataModel.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/EventSelection.h" -#include "TableHelper.h" -#include "Tools/ML/model.h" -#include "pidTPCBase.h" -#include "MetadataHelper.h" using namespace o2; using namespace o2::framework; @@ -155,8 +157,10 @@ struct tpcPid { void init(o2::framework::InitContext& initContext) { // Protection for process flags - if ((doprocessStandard && doprocessMcTuneOnData) || (!doprocessStandard && !doprocessMcTuneOnData)) { - LOG(fatal) << "pid-tpc must have only one of the options 'processStandard' OR 'processMcTuneOnData' enabled. Please check your configuration."; + if (!((doprocessStandard && !doprocessStandard2 && !doprocessMcTuneOnData) || + (!doprocessStandard && doprocessStandard2 && !doprocessMcTuneOnData) || + (!doprocessStandard && !doprocessStandard2 && doprocessMcTuneOnData))) { + LOG(fatal) << "pid-tpc must have only one of the options 'processStandard', 'processStandard2', 'processMcTuneOnData' enabled. Please check your configuration."; } response = new o2::pid::tpc::Response(); // Checking the tables are requested in the workflow and enabling them @@ -552,6 +556,98 @@ struct tpcPid { Partition mcnotTPCStandaloneTracks = (aod::track::tpcNClsFindable > static_cast(0)) && ((aod::track::itsClusterSizes > static_cast(0)) || (aod::track::trdPattern > static_cast(0)) || (aod::track::tofExpMom > 0.f && aod::track::tofChi2 > 0.f)); // To count number of tracks for use in NN array Partition mctracksWithTPC = (aod::track::tpcNClsFindable > (uint8_t)0); + void processStandard2(Coll const& collisions, Trks const& tracks, aod::DEdxsCorrected const& dedxscorrected, aod::BCsWithTimestamps const& bcs) + { + const uint64_t outTable_size = tracks.size(); + const uint64_t dedxscorrected_size = dedxscorrected.size(); + + if (dedxscorrected_size != outTable_size) { + LOG(fatal) << "Size of dEdx corrected table does not match size of tracks! dEdx size: " << dedxscorrected_size << ", tracks size: " << outTable_size; + } + + auto reserveTable = [&outTable_size](const Configurable& flag, auto& table) { + if (flag.value != 1) { + return; + } + table.reserve(outTable_size); + }; + + // Prepare memory for enabled tables + reserveTable(pidFullEl, tablePIDFullEl); + reserveTable(pidFullMu, tablePIDFullMu); + reserveTable(pidFullPi, tablePIDFullPi); + reserveTable(pidFullKa, tablePIDFullKa); + reserveTable(pidFullPr, tablePIDFullPr); + reserveTable(pidFullDe, tablePIDFullDe); + reserveTable(pidFullTr, tablePIDFullTr); + reserveTable(pidFullHe, tablePIDFullHe); + reserveTable(pidFullAl, tablePIDFullAl); + + reserveTable(pidTinyEl, tablePIDTinyEl); + reserveTable(pidTinyMu, tablePIDTinyMu); + reserveTable(pidTinyPi, tablePIDTinyPi); + reserveTable(pidTinyKa, tablePIDTinyKa); + reserveTable(pidTinyPr, tablePIDTinyPr); + reserveTable(pidTinyDe, tablePIDTinyDe); + reserveTable(pidTinyTr, tablePIDTinyTr); + reserveTable(pidTinyHe, tablePIDTinyHe); + reserveTable(pidTinyAl, tablePIDTinyAl); + + const uint64_t tracksForNet_size = (skipTPCOnly) ? notTPCStandaloneTracks.size() : tracksWithTPC.size(); + std::vector network_prediction; + + if (useNetworkCorrection) { + network_prediction = createNetworkPrediction(collisions, tracks, bcs, tracksForNet_size); + } + + uint64_t count_tracks = 0; + uint64_t count_tracks2 = 0; + + for (auto const& trk : tracks) { + // Loop on Tracks + + const auto& bc = trk.has_collision() ? collisions.iteratorAt(trk.collisionId()).bc_as() : bcs.begin(); + auto dedx_corr = dedxscorrected.iteratorAt(count_tracks2); + count_tracks2++; + if (useCCDBParam && ccdbTimestamp.value == 0 && !ccdb->isCachedObjectValid(ccdbPath.value, bc.timestamp())) { // Updating parametrisation only if the initial timestamp is 0 + if (recoPass.value == "") { + LOGP(info, "Retrieving latest TPC response object for timestamp {}:", bc.timestamp()); + } else { + LOGP(info, "Retrieving TPC Response for timestamp {} and recoPass {}:", bc.timestamp(), recoPass.value); + } + response = ccdb->getSpecific(ccdbPath.value, bc.timestamp(), metadata); + headers = ccdbApi.retrieveHeaders(ccdbPath.value, metadata, bc.timestamp()); + if (!response) { + LOGP(warning, "!! Could not find a valid TPC response object for specific pass name {}! Falling back to latest uploaded object.", metadata["RecoPassName"]); + response = ccdb->getForTimeStamp(ccdbPath.value, bc.timestamp()); + headers = ccdbApi.retrieveHeaders(ccdbPath.value, nullmetadata, bc.timestamp()); + if (!response) { + LOGP(fatal, "Could not find ANY TPC response object for the timestamp {}!", bc.timestamp()); + } + } + LOG(info) << "Successfully retrieved TPC PID object from CCDB for timestamp " << bc.timestamp() << ", period " << headers["LPMProductionTag"] << ", recoPass " << headers["RecoPassName"]; + response->PrintAll(); + } + auto makePidTablesDefault = [&trk, &dedx_corr, &collisions, &network_prediction, &count_tracks, &tracksForNet_size, this](const int flagFull, auto& tableFull, const int flagTiny, auto& tableTiny, const o2::track::PID::ID pid) { + makePidTables(flagFull, tableFull, flagTiny, tableTiny, pid, dedx_corr.tpcSignalCorrected(), trk, collisions, network_prediction, count_tracks, tracksForNet_size); + }; + + makePidTablesDefault(pidFullEl, tablePIDFullEl, pidTinyEl, tablePIDTinyEl, o2::track::PID::Electron); + makePidTablesDefault(pidFullMu, tablePIDFullMu, pidTinyMu, tablePIDTinyMu, o2::track::PID::Muon); + makePidTablesDefault(pidFullPi, tablePIDFullPi, pidTinyPi, tablePIDTinyPi, o2::track::PID::Pion); + makePidTablesDefault(pidFullKa, tablePIDFullKa, pidTinyKa, tablePIDTinyKa, o2::track::PID::Kaon); + makePidTablesDefault(pidFullPr, tablePIDFullPr, pidTinyPr, tablePIDTinyPr, o2::track::PID::Proton); + makePidTablesDefault(pidFullDe, tablePIDFullDe, pidTinyDe, tablePIDTinyDe, o2::track::PID::Deuteron); + makePidTablesDefault(pidFullTr, tablePIDFullTr, pidTinyTr, tablePIDTinyTr, o2::track::PID::Triton); + makePidTablesDefault(pidFullHe, tablePIDFullHe, pidTinyHe, tablePIDTinyHe, o2::track::PID::Helium3); + makePidTablesDefault(pidFullAl, tablePIDFullAl, pidTinyAl, tablePIDTinyAl, o2::track::PID::Alpha); + + if (trk.hasTPC() && (!skipTPCOnly || trk.hasITS() || trk.hasTRD() || trk.hasTOF())) { + count_tracks++; // Increment network track counter only if track has TPC, and (not skipping TPConly) or (is not TPConly) + } + } + } + PROCESS_SWITCH(tpcPid, processStandard2, "Creating PID tables with Corrected dEdx", false); void processMcTuneOnData(CollMC const& collisionsMc, TrksMC const& tracksMc, aod::BCsWithTimestamps const& bcs, aod::McParticles const&) { gRandom->SetSeed(0); // Ensure unique seed from UUID for each process call diff --git a/Common/TableProducer/PID/pidTPCBase.cxx b/Common/TableProducer/PID/pidTPCBase.cxx index cf82ea0c9d5..3ba7e881eff 100644 --- a/Common/TableProducer/PID/pidTPCBase.cxx +++ b/Common/TableProducer/PID/pidTPCBase.cxx @@ -15,18 +15,21 @@ /// \brief Base to build tasks for TPC PID tasks. /// +#include #include #include -#include // O2 includes -#include "CCDB/BasicCCDBManager.h" -#include "Framework/AnalysisTask.h" -#include "ReconstructionDataFormats/Track.h" -#include "Common/DataModel/FT0Corrected.h" #include "TableHelper.h" #include "pidTPCBase.h" + +#include "Common/DataModel/FT0Corrected.h" +#include "Common/CCDB/ctpRateFetcher.h" + +#include "CCDB/BasicCCDBManager.h" +#include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" using namespace o2; using namespace o2::framework; @@ -76,7 +79,106 @@ struct PidMultiplicity { PROCESS_SWITCH(PidMultiplicity, processStandard, "Process with tracks, needs propagated tracks", true); }; +struct DeDxCorrection { + Produces dEdxCorrected; + using BCsRun3 = soa::Join; + using ColEvSels = soa::Join; + using FullTracksIU = soa::Join; + + uint64_t minGlobalBC = 0; + Service ccdb; + ctpRateFetcher mRateFetcher; + + Str_dEdx_correction str_dedx_correction; + + // void init(InitContext& initContext) + void init(o2::framework::InitContext&) + { + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + str_dedx_correction.init(); + } + + void processRun3( + ColEvSels const& cols, + FullTracksIU const& tracks, + aod::BCsWithTimestamps const& bcs) + { + const uint64_t outTable_size = tracks.size(); + dEdxCorrected.reserve(outTable_size); + + for (auto const& trk : tracks) { + double hadronicRate; + int multTPC; + int occupancy; + if (trk.has_collision()) { + auto collision = cols.iteratorAt(trk.collisionId()); + auto bc = collision.bc_as(); + const int runnumber = bc.runNumber(); + hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), runnumber, "ZNC hadronic") * 1.e-3; // kHz + multTPC = collision.multTPC(); + occupancy = collision.trackOccupancyInTimeRange(); + } else { + auto bc = bcs.begin(); + const int runnumber = bc.runNumber(); + hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), runnumber, "ZNC hadronic") * 1.e-3; // kHz + multTPC = 0; + occupancy = 0; + } + + float signedP = trk.sign() * trk.tpcInnerParam(); + float fTPCSignal = trk.tpcSignal(); + float fNormMultTPC = multTPC / 11000.; + + float fTrackOccN = occupancy / 1000.; + float fOccTPCN = fNormMultTPC * 10; //(fNormMultTPC*10).clip(0,12) + if (fOccTPCN > 12) + fOccTPCN = 12; + else if (fOccTPCN < 0) + fOccTPCN = 0; + + float fTrackOccMeanN = hadronicRate / 5; + float side = trk.tgl() > 0 ? 1 : 0; + float a1pt = std::abs(trk.signed1Pt()); + float a1pt2 = a1pt * a1pt; + float atgl = std::abs(trk.tgl()); + float mbb0R = 50 / fTPCSignal; + if (mbb0R > 1.05) + mbb0R = 1.05; + else if (mbb0R < 0.05) + mbb0R = 0.05; + // float mbb0R = max(0.05, min(50 / fTPCSignal, 1.05)); + float a1ptmbb0R = a1pt * mbb0R; + float atglmbb0R = atgl * mbb0R; + + std::vector vec_occu = {fTrackOccN, fOccTPCN, fTrackOccMeanN}; + std::vector vec_track = {mbb0R, a1pt, atgl, atglmbb0R, a1ptmbb0R, side, a1pt2}; + + float fTPCSignalN_CR0 = str_dedx_correction.fReal_fTPCSignalN(vec_occu, vec_track); + + float mbb0R1 = 50 / (fTPCSignal / fTPCSignalN_CR0); + if (mbb0R1 > 1.05) + mbb0R1 = 1.05; + else if (mbb0R1 < 0.05) + mbb0R1 = 0.05; + + std::vector vec_track1 = {mbb0R1, a1pt, atgl, atgl * mbb0R1, a1pt * mbb0R1, side, a1pt2}; + float fTPCSignalN_CR1 = str_dedx_correction.fReal_fTPCSignalN(vec_occu, vec_track1); + + float corrected_dEdx = fTPCSignal / fTPCSignalN_CR1; + dEdxCorrected(corrected_dEdx); + } + } + PROCESS_SWITCH(DeDxCorrection, processRun3, "dEdx correction process", false); + + void processDummy(ColEvSels const&) {} + + PROCESS_SWITCH(DeDxCorrection, processDummy, "Do nothing", true); +}; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; } diff --git a/Common/TableProducer/PID/pidTPCBase.h b/Common/TableProducer/PID/pidTPCBase.h index 937eea21c48..f6017758519 100644 --- a/Common/TableProducer/PID/pidTPCBase.h +++ b/Common/TableProducer/PID/pidTPCBase.h @@ -21,13 +21,22 @@ #include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponseTPC.h" +#include "TMatrixD.h" + namespace o2::aod { +namespace pid +{ +DECLARE_SOA_COLUMN(TpcSignalCorrected, tpcSignalCorrected, float); //! +}; // namespace pid + DECLARE_SOA_TABLE(PIDMults, "AOD", "PIDMults", //! TPC auxiliary table for the PID o2::soa::Marker<1>, mult::MultTPC); +DECLARE_SOA_TABLE_FULL(DEdxsCorrected, "DEdxsCorrected", "AOD", "DEDXCORR", pid::TpcSignalCorrected); //! using PIDMult = PIDMults::iterator; +using DEdxCorrected = DEdxsCorrected::iterator; //! } // namespace o2::aod @@ -57,4 +66,41 @@ int getPIDIndex(const int pdgCode) // Get O2 PID index corresponding to MC PDG c } } +typedef struct Str_dEdx_correction { + TMatrixD fMatrix; + bool warning = true; + + // void init(std::vector& params) + void init() + { + double elements[32] = {0.99091, -0.015053, 0.0018912, -0.012305, + 0.081387, 0.003205, -0.0087404, -0.0028608, + 0.013066, 0.017012, -0.0018469, -0.0052177, + -0.0035655, 0.0017846, 0.0019127, -0.00012964, + 0.0049428, 0.0055592, -0.0010618, -0.0016134, + -0.0059098, 0.0013335, 0.00052133, 3.1119e-05, + -0.004882, 0.00077317, -0.0013827, 0.003249, + -0.00063689, 0.0016218, -0.00045215, -1.5815e-05}; + fMatrix.ResizeTo(4, 8); + fMatrix.SetMatrixArray(elements); + } + + float fReal_fTPCSignalN(std::vector& vec1, std::vector& vec2) + { + float result = 0.f; + // push 1. + vec1.insert(vec1.begin(), 1.0); + vec2.insert(vec2.begin(), 1.0); + for (int i = 0; i < fMatrix.GetNrows(); i++) { + for (int j = 0; j < fMatrix.GetNcols(); j++) { + double param = fMatrix(i, j); + double value1 = i > vec1.size() ? 0 : vec1[i]; + double value2 = j > vec2.size() ? 0 : vec2[j]; + result += param * value1 * value2; + } + } + return result; + } +} Str_dEdx_correction; + #endif // COMMON_TABLEPRODUCER_PID_PIDTPCBASE_H_ diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index 6c3d1304cfe..3c3747bd31e 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -17,29 +17,33 @@ /// \author Jeremy Wilkinson #include "tpcSkimsTableCreator.h" + #include + #include -#include #include +#include /// ROOT #include "TRandom3.h" /// O2 -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" #include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "Framework/runDataProcessing.h" /// O2Physics +#include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/trackUtilities.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/OccupancyTables.h" #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/PIDResponseITS.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "PWGLF/DataModel/LFStrangenessTables.h" -#include "PWGDQ/DataModel/ReducedInfoTables.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/CCDB/ctpRateFetcher.h" -#include "Common/DataModel/OccupancyTables.h" +#include "Common/TableProducer/PID/pidTPCBase.h" using namespace o2; using namespace o2::framework; @@ -52,6 +56,7 @@ struct TreeWriterTpcV0 { Service ccdb; using Trks = soa::Join; + using TrksWithDEdxCorrection = soa::Join; using Colls = soa::Join; using MyBCTable = soa::Join; using V0sWithID = soa::Join; @@ -89,7 +94,7 @@ struct TreeWriterTpcV0 { ctpRateFetcher mRateFetcher; /// Funktion to fill skimmed tables - template + template void fillSkimmedV0Table(V0 const& v0, T const& track, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, float hadronicRate) { @@ -111,7 +116,13 @@ struct TreeWriterTpcV0 { const double pseudoRndm = track.pt() * 1000. - static_cast(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { - rowTPCTree(track.tpcSignal(), + float usedDedx; + if constexpr (doUseCorreceddEdx) { + usedDedx = track.tpcSignalCorrected(); + } else { + usedDedx = track.tpcSignal(); + } + rowTPCTree(usedDedx, 1. / dEdxExp, track.tpcInnerParam(), track.tgl(), @@ -140,7 +151,7 @@ struct TreeWriterTpcV0 { } }; - template + template void fillSkimmedV0TableWithdEdxTrQA(V0 const& v0, T const& track, TQA const& trackQA, bool existTrkQA, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, float hadronicRate) { @@ -162,7 +173,13 @@ struct TreeWriterTpcV0 { const double pseudoRndm = track.pt() * 1000. - static_cast(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { - rowTPCTreeWithdEdxTrkQA(track.tpcSignal(), + float usedDedx; + if constexpr (doUseCorreceddEdx) { + usedDedx = track.tpcSignalCorrected(); + } else { + usedDedx = track.tpcSignal(); + } + rowTPCTreeWithdEdxTrkQA(usedDedx, 1. / dEdxExp, track.tpcInnerParam(), track.tgl(), @@ -193,7 +210,7 @@ struct TreeWriterTpcV0 { }; /// Function to fill skimmed tables - template + template void fillSkimmedV0TableWithTrQA(V0 const& v0, T const& track, TQA const& trackQA, bool existTrkQA, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, float hadronicRate, int bcGlobalIndex, int bcTimeFrameId, int bcBcInTimeFrame) { @@ -215,7 +232,13 @@ struct TreeWriterTpcV0 { const double pseudoRndm = track.pt() * 1000. - static_cast(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { - rowTPCTreeWithTrkQA(track.tpcSignal(), + float usedDedx; + if constexpr (doUseCorreceddEdx) { + usedDedx = track.tpcSignalCorrected(); + } else { + usedDedx = track.tpcSignal(); + } + rowTPCTreeWithTrkQA(usedDedx, 1. / dEdxExp, track.tpcInnerParam(), track.tgl(), @@ -376,6 +399,69 @@ struct TreeWriterTpcV0 { } /// process Standard PROCESS_SWITCH(TreeWriterTpcV0, processStandard, "Standard V0 Samples for PID", true); + void processStandardWithCorrecteddEdx(Colls::iterator const& collision, soa::Filtered const& tracks, V0sWithID const& v0s, aod::BCsWithTimestamps const&) + { + /// Check event slection + if (!isEventSelected(collision, tracks)) { + return; + } + auto bc = collision.bc_as(); + const int runnumber = bc.runNumber(); + float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), runnumber, irSource) * 1.e-3; + + rowTPCTree.reserve(tracks.size()); + + /// Loop over v0 candidates + for (const auto& v0 : v0s) { + auto posTrack = v0.posTrack_as>(); + auto negTrack = v0.negTrack_as>(); + if (v0.v0addid() == -1) { + continue; + } + // gamma + if (static_cast(posTrack.pidbit() & (1 << 0)) && static_cast(negTrack.pidbit() & (1 << 0))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisElectrons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Electron], maxPt4dwnsmplTsalisElectrons)) { + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaEl(), posTrack.tofNSigmaEl(), posTrack.tpcExpSignalEl(posTrack.tpcSignalCorrected()), o2::track::PID::Electron, runnumber, dwnSmplFactor_El, hadronicRate); + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisElectrons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Electron], maxPt4dwnsmplTsalisElectrons)) { + fillSkimmedV0Table(v0, negTrack, collision, negTrack.tpcNSigmaEl(), negTrack.tofNSigmaEl(), negTrack.tpcExpSignalEl(negTrack.tpcSignalCorrected()), o2::track::PID::Electron, runnumber, dwnSmplFactor_El, hadronicRate); + } + } + // Ks0 + if (static_cast(posTrack.pidbit() & (1 << 1)) && static_cast(negTrack.pidbit() & (1 << 1))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0Table(v0, negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(negTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + } + // Lambda + if (static_cast(posTrack.pidbit() & (1 << 2)) && static_cast(negTrack.pidbit() & (1 << 2))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], maxPt4dwnsmplTsalisProtons)) { + if (std::abs(posTrack.tofNSigmaPr()) <= nSigmaTOFdautrack) { + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaPr(), posTrack.tofNSigmaPr(), posTrack.tpcExpSignalPr(posTrack.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate); + } + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0Table(v0, negTrack, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(negTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + } + // Antilambda + if (static_cast(posTrack.pidbit() & (1 << 3)) && static_cast(negTrack.pidbit() & (1 << 3))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0Table(v0, posTrack, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], maxPt4dwnsmplTsalisProtons)) { + if (std::abs(negTrack.tofNSigmaPr()) <= nSigmaTOFdautrack) { + fillSkimmedV0Table(v0, negTrack, collision, negTrack.tpcNSigmaPr(), negTrack.tofNSigmaPr(), negTrack.tpcExpSignalPr(negTrack.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate); + } + } + } + } + } /// process Standard + PROCESS_SWITCH(TreeWriterTpcV0, processStandardWithCorrecteddEdx, "Standard V0 Samples for PID with corrected dEdx", false); + Preslice perCollisionTracks = aod::track::collisionId; Preslice perCollisionV0s = aod::v0data::collisionId; void processWithdEdxTrQA(Colls const& collisions, Trks const& myTracks, V0sWithID const& myV0s, aod::BCsWithTimestamps const&, aod::TracksQAVersion const& tracksQA) @@ -470,6 +556,99 @@ struct TreeWriterTpcV0 { } /// process with dEdx from TrackQA PROCESS_SWITCH(TreeWriterTpcV0, processWithdEdxTrQA, "Standard V0 Samples with dEdx from Track QA for PID", false); + Preslice perCollisionTracksWithNewDEdx = aod::track::collisionId; + void processWithdEdxTrQAWithCorrecteddEdx(Colls const& collisions, TrksWithDEdxCorrection const& myTracks, V0sWithID const& myV0s, aod::BCsWithTimestamps const&, aod::TracksQAVersion const& tracksQA) + { + std::vector labelTrack2TrackQA; + labelTrack2TrackQA.clear(); + labelTrack2TrackQA.resize(myTracks.size(), -1); + for (const auto& trackQA : tracksQA) { + int64_t trackId = trackQA.trackId(); + int64_t trackQAIndex = trackQA.globalIndex(); + labelTrack2TrackQA[trackId] = trackQAIndex; + } + for (const auto& collision : collisions) { + auto tracks = myTracks.sliceBy(perCollisionTracksWithNewDEdx, collision.globalIndex()); + auto v0s = myV0s.sliceBy(perCollisionV0s, collision.globalIndex()); + /// Check event slection + if (!isEventSelected(collision, tracks)) { + continue; + } + auto bc = collision.bc_as(); + const int runnumber = bc.runNumber(); + float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), runnumber, irSource) * 1.e-3; + rowTPCTreeWithTrkQA.reserve(tracks.size()); + /// Loop over v0 candidates + for (const auto& v0 : v0s) { + auto posTrack = v0.posTrack_as(); + auto negTrack = v0.negTrack_as(); + if (v0.v0addid() == -1) { + continue; + } + aod::TracksQA posTrackQA; + aod::TracksQA negTrackQA; + bool existPosTrkQA; + bool existNegTrkQA; + if (labelTrack2TrackQA[posTrack.globalIndex()] != -1) { + posTrackQA = tracksQA.iteratorAt(labelTrack2TrackQA[posTrack.globalIndex()]); + existPosTrkQA = true; + } else { + posTrackQA = tracksQA.iteratorAt(0); + existPosTrkQA = false; + } + if (labelTrack2TrackQA[negTrack.globalIndex()] != -1) { + negTrackQA = tracksQA.iteratorAt(labelTrack2TrackQA[negTrack.globalIndex()]); + existNegTrkQA = true; + } else { + negTrackQA = tracksQA.iteratorAt(0); + existNegTrkQA = false; + } + + // gamma + if (static_cast(posTrack.pidbit() & (1 << 0)) && static_cast(negTrack.pidbit() & (1 << 0))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisElectrons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Electron], maxPt4dwnsmplTsalisElectrons)) { + fillSkimmedV0TableWithdEdxTrQA(v0, posTrack, posTrackQA, existPosTrkQA, collision, posTrack.tpcNSigmaEl(), posTrack.tofNSigmaEl(), posTrack.tpcExpSignalEl(posTrack.tpcSignalCorrected()), o2::track::PID::Electron, runnumber, dwnSmplFactor_El, hadronicRate); + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisElectrons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Electron], maxPt4dwnsmplTsalisElectrons)) { + fillSkimmedV0TableWithdEdxTrQA(v0, negTrack, negTrackQA, existNegTrkQA, collision, negTrack.tpcNSigmaEl(), negTrack.tofNSigmaEl(), negTrack.tpcExpSignalEl(negTrack.tpcSignalCorrected()), o2::track::PID::Electron, runnumber, dwnSmplFactor_El, hadronicRate); + } + } + // Ks0 + if (static_cast(posTrack.pidbit() & (1 << 1)) && static_cast(negTrack.pidbit() & (1 << 1))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0TableWithdEdxTrQA(v0, posTrack, posTrackQA, existPosTrkQA, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0TableWithdEdxTrQA(v0, negTrack, negTrackQA, existNegTrkQA, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(negTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + } + // Lambda + if (static_cast(posTrack.pidbit() & (1 << 2)) && static_cast(negTrack.pidbit() & (1 << 2))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], maxPt4dwnsmplTsalisProtons)) { + if (std::abs(posTrack.tofNSigmaPr()) <= nSigmaTOFdautrack) { + fillSkimmedV0TableWithdEdxTrQA(v0, posTrack, posTrackQA, existPosTrkQA, collision, posTrack.tpcNSigmaPr(), posTrack.tofNSigmaPr(), posTrack.tpcExpSignalPr(posTrack.tpcSignal()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate); + } + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0TableWithdEdxTrQA(v0, negTrack, negTrackQA, existNegTrkQA, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(negTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + } + // Antilambda + if (static_cast(posTrack.pidbit() & (1 << 3)) && static_cast(negTrack.pidbit() & (1 << 3))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0TableWithdEdxTrQA(v0, posTrack, posTrackQA, existPosTrkQA, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], maxPt4dwnsmplTsalisProtons)) { + if (std::abs(negTrack.tofNSigmaPr()) <= nSigmaTOFdautrack) { + fillSkimmedV0TableWithdEdxTrQA(v0, negTrack, negTrackQA, existNegTrkQA, collision, negTrack.tpcNSigmaPr(), negTrack.tofNSigmaPr(), negTrack.tpcExpSignalPr(negTrack.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate); + } + } + } + } + } + } /// process with dEdx from TrackQA + PROCESS_SWITCH(TreeWriterTpcV0, processWithdEdxTrQAWithCorrecteddEdx, "Standard V0 Samples with dEdx from Track QA for PID with corrected dEdx", false); + void processWithTrQA(Colls const& collisions, Trks const& myTracks, V0sWithID const& myV0s, MyBCTable const&, aod::TracksQAVersion const& tracksQA) { std::vector labelTrack2TrackQA; @@ -565,10 +744,105 @@ struct TreeWriterTpcV0 { } /// process with TrackQA PROCESS_SWITCH(TreeWriterTpcV0, processWithTrQA, "Standard V0 Samples with Track QA for PID", false); + void processWithTrQAWithCorrecteddEdx(Colls const& collisions, TrksWithDEdxCorrection const& myTracks, V0sWithID const& myV0s, MyBCTable const&, aod::TracksQAVersion const& tracksQA) + { + std::vector labelTrack2TrackQA; + labelTrack2TrackQA.clear(); + labelTrack2TrackQA.resize(myTracks.size(), -1); + for (const auto& trackQA : tracksQA) { + int64_t trackId = trackQA.trackId(); + int64_t trackQAIndex = trackQA.globalIndex(); + labelTrack2TrackQA[trackId] = trackQAIndex; + } + for (const auto& collision : collisions) { + auto tracks = myTracks.sliceBy(perCollisionTracksWithNewDEdx, collision.globalIndex()); + auto v0s = myV0s.sliceBy(perCollisionV0s, collision.globalIndex()); + /// Check event slection + if (!isEventSelected(collision, tracks)) { + continue; + } + auto bc = collision.bc_as(); + const int runnumber = bc.runNumber(); + float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), runnumber, irSource) * 1.e-3; + const int bcGlobalIndex = bc.globalIndex(); + const int bcTimeFrameId = bc.tfId(); + const int bcBcInTimeFrame = bc.bcInTF(); + rowTPCTreeWithTrkQA.reserve(tracks.size()); + /// Loop over v0 candidates + for (const auto& v0 : v0s) { + auto posTrack = v0.posTrack_as(); + auto negTrack = v0.negTrack_as(); + if (v0.v0addid() == -1) { + continue; + } + aod::TracksQA posTrackQA; + aod::TracksQA negTrackQA; + bool existPosTrkQA; + bool existNegTrkQA; + if (labelTrack2TrackQA[posTrack.globalIndex()] != -1) { + posTrackQA = tracksQA.iteratorAt(labelTrack2TrackQA[posTrack.globalIndex()]); + existPosTrkQA = true; + } else { + posTrackQA = tracksQA.iteratorAt(0); + existPosTrkQA = false; + } + if (labelTrack2TrackQA[negTrack.globalIndex()] != -1) { + negTrackQA = tracksQA.iteratorAt(labelTrack2TrackQA[negTrack.globalIndex()]); + existNegTrkQA = true; + } else { + negTrackQA = tracksQA.iteratorAt(0); + existNegTrkQA = false; + } + + // gamma + if (static_cast(posTrack.pidbit() & (1 << 0)) && static_cast(negTrack.pidbit() & (1 << 0))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisElectrons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Electron], maxPt4dwnsmplTsalisElectrons)) { + fillSkimmedV0TableWithTrQA(v0, posTrack, posTrackQA, existPosTrkQA, collision, posTrack.tpcNSigmaEl(), posTrack.tofNSigmaEl(), posTrack.tpcExpSignalEl(posTrack.tpcSignalCorrected()), o2::track::PID::Electron, runnumber, dwnSmplFactor_El, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisElectrons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Electron], maxPt4dwnsmplTsalisElectrons)) { + fillSkimmedV0TableWithTrQA(v0, negTrack, negTrackQA, existNegTrkQA, collision, negTrack.tpcNSigmaEl(), negTrack.tofNSigmaEl(), negTrack.tpcExpSignalEl(negTrack.tpcSignalCorrected()), o2::track::PID::Electron, runnumber, dwnSmplFactor_El, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + } + // Ks0 + if (static_cast(posTrack.pidbit() & (1 << 1)) && static_cast(negTrack.pidbit() & (1 << 1))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0TableWithTrQA(v0, posTrack, posTrackQA, existPosTrkQA, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0TableWithTrQA(v0, negTrack, negTrackQA, existNegTrkQA, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(negTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + } + // Lambda + if (static_cast(posTrack.pidbit() & (1 << 2)) && static_cast(negTrack.pidbit() & (1 << 2))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], maxPt4dwnsmplTsalisProtons)) { + if (std::abs(posTrack.tofNSigmaPr()) <= nSigmaTOFdautrack) { + fillSkimmedV0TableWithTrQA(v0, posTrack, posTrackQA, existPosTrkQA, collision, posTrack.tpcNSigmaPr(), posTrack.tofNSigmaPr(), posTrack.tpcExpSignalPr(posTrack.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0TableWithTrQA(v0, negTrack, negTrackQA, existNegTrkQA, collision, negTrack.tpcNSigmaPi(), negTrack.tofNSigmaPi(), negTrack.tpcExpSignalPi(negTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + } + // Antilambda + if (static_cast(posTrack.pidbit() & (1 << 3)) && static_cast(negTrack.pidbit() & (1 << 3))) { + if (downsampleTsalisCharged(posTrack.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion], maxPt4dwnsmplTsalisPions)) { + fillSkimmedV0TableWithTrQA(v0, posTrack, posTrackQA, existPosTrkQA, collision, posTrack.tpcNSigmaPi(), posTrack.tofNSigmaPi(), posTrack.tpcExpSignalPi(posTrack.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + if (downsampleTsalisCharged(negTrack.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton], maxPt4dwnsmplTsalisProtons)) { + if (std::abs(negTrack.tofNSigmaPr()) <= nSigmaTOFdautrack) { + fillSkimmedV0TableWithTrQA(v0, negTrack, negTrackQA, existNegTrkQA, collision, negTrack.tpcNSigmaPr(), negTrack.tofNSigmaPr(), negTrack.tpcExpSignalPr(negTrack.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + } + } + } + } + } /// process with TrackQA + PROCESS_SWITCH(TreeWriterTpcV0, processWithTrQAWithCorrecteddEdx, "Standard V0 Samples with Track QA for PID with corrected dEdx", false); + void processDummy(Colls const&) {} PROCESS_SWITCH(TreeWriterTpcV0, processDummy, "Dummy function", false); -}; /// struct TreeWriterTpcV0 +}; /// struct TreeWriterTpcV0 struct TreeWriterTPCTOF { @@ -580,6 +854,12 @@ struct TreeWriterTPCTOF { aod::pidTOFFullEl, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr, aod::pidTOFFullDe, aod::pidTOFFullTr, aod::TrackSelection>; + using TrksWithDEdxCorrection = soa::Join; using Colls = soa::Join; using MyBCTable = soa::Join; @@ -674,7 +954,7 @@ struct TreeWriterTPCTOF { }; /// Function to fill trees - template + template void fillSkimmedTPCTOFTable(T const& track, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, double hadronicRate) { @@ -689,7 +969,13 @@ struct TreeWriterTPCTOF { const double pseudoRndm = track.pt() * 1000. - static_cast(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { - rowTPCTOFTree(track.tpcSignal(), + float usedEdx; + if constexpr (doCorrectdEdx) { + usedEdx = track.tpcSignalCorrected(); + } else { + usedEdx = track.tpcSignal(); + } + rowTPCTOFTree(usedEdx, 1. / dEdxExp, track.tpcInnerParam(), track.tgl(), @@ -711,7 +997,7 @@ struct TreeWriterTPCTOF { hadronicRate); } }; - template + template void fillSkimmedTPCTOFTableWithdEdxTrkQA(T const& track, TQA const& trackQA, bool existTrkQA, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float nSigmaITS, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, double hadronicRate) { @@ -726,7 +1012,13 @@ struct TreeWriterTPCTOF { const double pseudoRndm = track.pt() * 1000. - static_cast(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { - rowTPCTOFTreeWithdEdxTrkQA(track.tpcSignal(), + float usedEdx; + if constexpr (doCorrectdEdx) { + usedEdx = track.tpcSignalCorrected(); + } else { + usedEdx = track.tpcSignal(); + } + rowTPCTOFTreeWithdEdxTrkQA(usedEdx, 1. / dEdxExp, track.tpcInnerParam(), track.tgl(), @@ -751,7 +1043,7 @@ struct TreeWriterTPCTOF { } }; /// Function to fill trees - template + template void fillSkimmedTPCTOFTableWithTrkQA(T const& track, TQA const& trackQA, bool existTrkQA, C const& collision, const float nSigmaTPC, const float nSigmaTOF, const float nSigmaITS, const float dEdxExp, const o2::track::PID::ID id, int runnumber, double dwnSmplFactor, double hadronicRate, int bcGlobalIndex, int bcTimeFrameId, int bcBcInTimeFrame) { @@ -766,7 +1058,13 @@ struct TreeWriterTPCTOF { const double pseudoRndm = track.pt() * 1000. - static_cast(track.pt() * 1000); if (pseudoRndm < dwnSmplFactor) { - rowTPCTOFTreeWithTrkQA(track.tpcSignal(), + float usedEdx; + if constexpr (doCorrectdEdx) { + usedEdx = track.tpcSignalCorrected(); + } else { + usedEdx = track.tpcSignal(); + } + rowTPCTOFTreeWithTrkQA(usedEdx, 1. / dEdxExp, track.tpcInnerParam(), track.tgl(), @@ -869,8 +1167,55 @@ struct TreeWriterTPCTOF { fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignal()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); } } /// Loop tracks - } /// process + } /// process PROCESS_SWITCH(TreeWriterTPCTOF, processStandard, "Standard Samples for PID", true); + + void processStandard2(Colls::iterator const& collision, soa::Filtered const& tracks, aod::BCsWithTimestamps const&) + { + /// Check event selection + if (!isEventSelected(collision, tracks)) { + return; + } + auto bc = collision.bc_as(); + const int runnumber = bc.runNumber(); + float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), runnumber, irSource) * 1.e-3; + + rowTPCTOFTree.reserve(tracks.size()); + for (auto const& trk : tracks) { + /// Fill tree for tritons/* */ + if (trk.tpcInnerParam() < maxMomHardCutOnlyTr && trk.tpcInnerParam() <= maxMomTPCOnlyTr && std::abs(trk.tpcNSigmaTr()) < nSigmaTPCOnlyTr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Triton])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaTr(), trk.tofNSigmaTr(), trk.tpcExpSignalTr(trk.tpcSignalCorrected()), o2::track::PID::Triton, runnumber, dwnSmplFactor_Tr, hadronicRate); + } else if (trk.tpcInnerParam() < maxMomHardCutOnlyTr && trk.tpcInnerParam() > maxMomTPCOnlyTr && std::abs(trk.tofNSigmaTr()) < nSigmaTOF_TPCTOF_Tr && std::abs(trk.tpcNSigmaTr()) < nSigmaTPC_TPCTOF_Tr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Triton])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaTr(), trk.tofNSigmaTr(), trk.tpcExpSignalTr(trk.tpcSignalCorrected()), o2::track::PID::Triton, runnumber, dwnSmplFactor_Tr, hadronicRate); + } + /// Fill tree for deuterons + if (trk.tpcInnerParam() < maxMomHardCutOnlyDe && trk.tpcInnerParam() <= maxMomTPCOnlyDe && std::abs(trk.tpcNSigmaDe()) < nSigmaTPCOnlyDe && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Deuteron])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaDe(), trk.tofNSigmaDe(), trk.tpcExpSignalDe(trk.tpcSignalCorrected()), o2::track::PID::Deuteron, runnumber, dwnSmplFactor_De, hadronicRate); + } else if (trk.tpcInnerParam() < maxMomHardCutOnlyDe && trk.tpcInnerParam() > maxMomTPCOnlyDe && std::abs(trk.tofNSigmaDe()) < nSigmaTOF_TPCTOF_De && std::abs(trk.tpcNSigmaDe()) < nSigmaTPC_TPCTOF_De && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Deuteron])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaDe(), trk.tofNSigmaDe(), trk.tpcExpSignalDe(trk.tpcSignalCorrected()), o2::track::PID::Deuteron, runnumber, dwnSmplFactor_De, hadronicRate); + } + /// Fill tree for protons + if (trk.tpcInnerParam() <= maxMomTPCOnlyPr && std::abs(trk.tpcNSigmaPr()) < nSigmaTPCOnlyPr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.tpcExpSignalPr(trk.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && std::abs(trk.tofNSigmaPr()) < nSigmaTOF_TPCTOF_Pr && std::abs(trk.tpcNSigmaPr()) < nSigmaTPC_TPCTOF_Pr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.tpcExpSignalPr(trk.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate); + } + /// Fill tree for kaons + if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() <= maxMomTPCOnlyKa && std::abs(trk.tpcNSigmaKa()) < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignalCorrected()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + } else if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() > maxMomTPCOnlyKa && std::abs(trk.tofNSigmaKa()) < nSigmaTOF_TPCTOF_Ka && std::abs(trk.tpcNSigmaKa()) < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignalCorrected()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + } + /// Fill tree pions + if (trk.tpcInnerParam() <= maxMomTPCOnlyPi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && std::abs(trk.tofNSigmaPi()) < nSigmaTOF_TPCTOF_Pi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPC_TPCTOF_Pi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + fillSkimmedTPCTOFTable(trk, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + } /// Loop tracks + } /// process + PROCESS_SWITCH(TreeWriterTPCTOF, processStandard2, "Standard Samples for PID with corrected dEdx", false); + Preslice perCollisionTracks = aod::track::collisionId; void processWithdEdxTrQA(Colls const& collisions, Trks const& myTracks, aod::BCsWithTimestamps const&, aod::TracksQAVersion const& tracksQA) { @@ -937,7 +1282,7 @@ struct TreeWriterTPCTOF { if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() <= maxMomTPCOnlyKa && std::abs(trk.tpcNSigmaKa()) < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); } else if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() > maxMomTPCOnlyKa && std::abs(trk.tofNSigmaKa()) < nSigmaTOF_TPCTOF_Ka && std::abs(trk.tpcNSigmaKa()) < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { - fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + (trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); } /// Fill tree pions if (trk.tpcInnerParam() <= maxMomTPCOnlyPi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { @@ -949,6 +1294,86 @@ struct TreeWriterTPCTOF { } } /// process PROCESS_SWITCH(TreeWriterTPCTOF, processWithdEdxTrQA, "Samples for PID with TrackQA info", false); + + Preslice perCollisionTracksWithCorrecteddEdx = aod::track::collisionId; + void processWithdEdxTrQAWithCorrecteddEdx(Colls const& collisions, TrksWithDEdxCorrection const& myTracks, aod::BCsWithTimestamps const&, aod::TracksQAVersion const& tracksQA) + { + std::vector labelTrack2TrackQA; + labelTrack2TrackQA.clear(); + labelTrack2TrackQA.resize(myTracks.size(), -1); + for (const auto& trackQA : tracksQA) { + int64_t trackId = trackQA.trackId(); + int64_t trackQAIndex = trackQA.globalIndex(); + labelTrack2TrackQA[trackId] = trackQAIndex; + } + for (const auto& collision : collisions) { + auto tracks = myTracks.sliceBy(perCollisionTracksWithCorrecteddEdx, collision.globalIndex()); + auto tracksWithITSPid = soa::Attach(tracks); + /// Check event selection + if (!isEventSelected(collision, tracks)) { + continue; + } + auto bc = collision.bc_as(); + const int runnumber = bc.runNumber(); + float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), runnumber, irSource) * 1.e-3; + rowTPCTOFTreeWithTrkQA.reserve(tracks.size()); + for (auto const& trk : tracksWithITSPid) { + if (!((trackSelection == 0) || + ((trackSelection == 1) && trk.isGlobalTrack()) || + ((trackSelection == 2) && trk.isGlobalTrackWoPtEta()) || + ((trackSelection == 3) && trk.isGlobalTrackWoDCA()) || + ((trackSelection == 4) && trk.isQualityTrack()) || + ((trackSelection == 5) && trk.isInAcceptanceTrack()))) { + continue; + } + // get the corresponding trackQA using labelTracks2TracKQA and get variables of interest + aod::TracksQA trackQA; + bool existTrkQA; + if (labelTrack2TrackQA[trk.globalIndex()] != -1) { + trackQA = tracksQA.iteratorAt(labelTrack2TrackQA[trk.globalIndex()]); + existTrkQA = true; + } else { + trackQA = tracksQA.iteratorAt(0); + existTrkQA = false; + } + /// Fill tree for tritons + if (trk.tpcInnerParam() < maxMomHardCutOnlyTr && trk.tpcInnerParam() <= maxMomTPCOnlyTr && std::abs(trk.tpcNSigmaTr()) < nSigmaTPCOnlyTr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Triton])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaTr(), trk.tofNSigmaTr(), trk.itsNSigmaTr(), trk.tpcExpSignalTr(trk.tpcSignalCorrected()), o2::track::PID::Triton, runnumber, dwnSmplFactor_Tr, hadronicRate); + } else if (trk.tpcInnerParam() < maxMomHardCutOnlyTr && trk.tpcInnerParam() > maxMomTPCOnlyTr && std::abs(trk.tofNSigmaTr()) < nSigmaTOF_TPCTOF_Tr && std::abs(trk.tpcNSigmaTr()) < nSigmaTPC_TPCTOF_Tr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Triton])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaTr(), trk.tofNSigmaTr(), trk.itsNSigmaTr(), trk.tpcExpSignalTr(trk.tpcSignalCorrected()), o2::track::PID::Triton, runnumber, dwnSmplFactor_Tr, hadronicRate); + } + /// Fill tree for deuterons + if (trk.tpcInnerParam() < maxMomHardCutOnlyDe && trk.tpcInnerParam() <= maxMomTPCOnlyDe && std::abs(trk.tpcNSigmaDe()) < nSigmaTPCOnlyDe && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Deuteron])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaDe(), trk.tofNSigmaDe(), trk.itsNSigmaDe(), trk.tpcExpSignalDe(trk.tpcSignalCorrected()), o2::track::PID::Deuteron, runnumber, dwnSmplFactor_De, hadronicRate); + } else if (trk.tpcInnerParam() < maxMomHardCutOnlyDe && trk.tpcInnerParam() > maxMomTPCOnlyDe && std::abs(trk.tofNSigmaDe()) < nSigmaTOF_TPCTOF_De && std::abs(trk.tpcNSigmaDe()) < nSigmaTPC_TPCTOF_De && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Deuteron])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaDe(), trk.tofNSigmaDe(), trk.itsNSigmaDe(), trk.tpcExpSignalDe(trk.tpcSignalCorrected()), o2::track::PID::Deuteron, runnumber, dwnSmplFactor_De, hadronicRate); + } + /// Fill tree for protons + if (trk.tpcInnerParam() <= maxMomTPCOnlyPr && std::abs(trk.tpcNSigmaPr()) < nSigmaTPCOnlyPr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.itsNSigmaPr(), trk.tpcExpSignalPr(trk.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && std::abs(trk.tofNSigmaPr()) < nSigmaTOF_TPCTOF_Pr && std::abs(trk.tpcNSigmaPr()) < nSigmaTPC_TPCTOF_Pr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.itsNSigmaPr(), trk.tpcExpSignalPr(trk.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate); + } + /// Fill tree for kaons + if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() <= maxMomTPCOnlyKa && std::abs(trk.tpcNSigmaKa()) < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignalCorrected()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + } else if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() > maxMomTPCOnlyKa && std::abs(trk.tofNSigmaKa()) < nSigmaTOF_TPCTOF_Ka && std::abs(trk.tpcNSigmaKa()) < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignalCorrected()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + } + /// Fill tree pions + if (trk.tpcInnerParam() <= maxMomTPCOnlyPi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.itsNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && std::abs(trk.tofNSigmaPi()) < nSigmaTOF_TPCTOF_Pi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPC_TPCTOF_Pi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.itsNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate); + } + } /// Loop tracks + } + } /// process + PROCESS_SWITCH(TreeWriterTPCTOF, processWithdEdxTrQAWithCorrecteddEdx, "Samples for PID with TrackQA info with corrected dEdx", false); + void processWithTrQA(Colls const& collisions, Trks const& myTracks, MyBCTable const&, aod::TracksQAVersion const& tracksQA) { std::vector labelTrack2TrackQA; @@ -1029,6 +1454,88 @@ struct TreeWriterTPCTOF { } } /// process PROCESS_SWITCH(TreeWriterTPCTOF, processWithTrQA, "Samples for PID with TrackQA info", false); + + void processWithTrQAWithCorrecteddEdx(Colls const& collisions, TrksWithDEdxCorrection const& myTracks, MyBCTable const&, aod::TracksQAVersion const& tracksQA) + { + std::vector labelTrack2TrackQA; + labelTrack2TrackQA.clear(); + labelTrack2TrackQA.resize(myTracks.size(), -1); + for (const auto& trackQA : tracksQA) { + int64_t trackId = trackQA.trackId(); + int64_t trackQAIndex = trackQA.globalIndex(); + labelTrack2TrackQA[trackId] = trackQAIndex; + } + for (const auto& collision : collisions) { + auto tracks = myTracks.sliceBy(perCollisionTracks, collision.globalIndex()); + auto tracksWithITSPid = soa::Attach(tracks); + /// Check event selection + if (!isEventSelected(collision, tracks)) { + continue; + } + auto bc = collision.bc_as(); + const int runnumber = bc.runNumber(); + float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), runnumber, irSource) * 1.e-3; + const int bcGlobalIndex = bc.globalIndex(); + const int bcTimeFrameId = bc.tfId(); + const int bcBcInTimeFrame = bc.bcInTF(); + rowTPCTOFTreeWithTrkQA.reserve(tracks.size()); + for (auto const& trk : tracksWithITSPid) { + if (!((trackSelection == 0) || + ((trackSelection == 1) && trk.isGlobalTrack()) || + ((trackSelection == 2) && trk.isGlobalTrackWoPtEta()) || + ((trackSelection == 3) && trk.isGlobalTrackWoDCA()) || + ((trackSelection == 4) && trk.isQualityTrack()) || + ((trackSelection == 5) && trk.isInAcceptanceTrack()))) { + continue; + } + // get the corresponding trackQA using labelTracks2TracKQA and get variables of interest + aod::TracksQA trackQA; + bool existTrkQA; + if (labelTrack2TrackQA[trk.globalIndex()] != -1) { + trackQA = tracksQA.iteratorAt(labelTrack2TrackQA[trk.globalIndex()]); + existTrkQA = true; + } else { + trackQA = tracksQA.iteratorAt(0); + existTrkQA = false; + } + /// Fill tree for tritons + if (trk.tpcInnerParam() < maxMomHardCutOnlyTr && trk.tpcInnerParam() <= maxMomTPCOnlyTr && std::abs(trk.tpcNSigmaTr()) < nSigmaTPCOnlyTr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Triton])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaTr(), trk.tofNSigmaTr(), trk.itsNSigmaTr(), trk.tpcExpSignalTr(trk.tpcSignalCorrected()), o2::track::PID::Triton, runnumber, dwnSmplFactor_Tr, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } else if (trk.tpcInnerParam() < maxMomHardCutOnlyTr && trk.tpcInnerParam() > maxMomTPCOnlyTr && std::abs(trk.tofNSigmaTr()) < nSigmaTOF_TPCTOF_Tr && std::abs(trk.tpcNSigmaTr()) < nSigmaTPC_TPCTOF_Tr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Triton])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaTr(), trk.tofNSigmaTr(), trk.itsNSigmaTr(), trk.tpcExpSignalTr(trk.tpcSignalCorrected()), o2::track::PID::Triton, runnumber, dwnSmplFactor_Tr, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + /// Fill tree for deuterons + if (trk.tpcInnerParam() < maxMomHardCutOnlyDe && trk.tpcInnerParam() <= maxMomTPCOnlyDe && std::abs(trk.tpcNSigmaDe()) < nSigmaTPCOnlyDe && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Deuteron])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaDe(), trk.tofNSigmaDe(), trk.itsNSigmaDe(), trk.tpcExpSignalDe(trk.tpcSignalCorrected()), o2::track::PID::Deuteron, runnumber, dwnSmplFactor_De, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } else if (trk.tpcInnerParam() < maxMomHardCutOnlyDe && trk.tpcInnerParam() > maxMomTPCOnlyDe && std::abs(trk.tofNSigmaDe()) < nSigmaTOF_TPCTOF_De && std::abs(trk.tpcNSigmaDe()) < nSigmaTPC_TPCTOF_De && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Deuteron])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaDe(), trk.tofNSigmaDe(), trk.itsNSigmaDe(), trk.tpcExpSignalDe(trk.tpcSignalCorrected()), o2::track::PID::Deuteron, runnumber, dwnSmplFactor_De, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + /// Fill tree for protons + if (trk.tpcInnerParam() <= maxMomTPCOnlyPr && std::abs(trk.tpcNSigmaPr()) < nSigmaTPCOnlyPr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.itsNSigmaPr(), trk.tpcExpSignalPr(trk.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPr && std::abs(trk.tofNSigmaPr()) < nSigmaTOF_TPCTOF_Pr && std::abs(trk.tpcNSigmaPr()) < nSigmaTPC_TPCTOF_Pr && downsampleTsalisCharged(trk.pt(), downsamplingTsalisProtons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Proton])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaPr(), trk.tofNSigmaPr(), trk.itsNSigmaPr(), trk.tpcExpSignalPr(trk.tpcSignalCorrected()), o2::track::PID::Proton, runnumber, dwnSmplFactor_Pr, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + /// Fill tree for kaons + if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() <= maxMomTPCOnlyKa && std::abs(trk.tpcNSigmaKa()) < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignalCorrected()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } else if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() > maxMomTPCOnlyKa && std::abs(trk.tofNSigmaKa()) < nSigmaTOF_TPCTOF_Ka && std::abs(trk.tpcNSigmaKa()) < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignalCorrected()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + /// Fill tree pions + if (trk.tpcInnerParam() <= maxMomTPCOnlyPi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.itsNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } else if (trk.tpcInnerParam() > maxMomTPCOnlyPi && std::abs(trk.tofNSigmaPi()) < nSigmaTOF_TPCTOF_Pi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPC_TPCTOF_Pi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) { + fillSkimmedTPCTOFTableWithTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaPi(), trk.tofNSigmaPi(), trk.itsNSigmaPi(), trk.tpcExpSignalPi(trk.tpcSignalCorrected()), o2::track::PID::Pion, runnumber, dwnSmplFactor_Pi, hadronicRate, bcGlobalIndex, bcTimeFrameId, bcBcInTimeFrame); + } + } /// Loop tracks + } + } /// process + PROCESS_SWITCH(TreeWriterTPCTOF, processWithTrQAWithCorrecteddEdx, "Samples for PID with TrackQA info with correced dEdx", false); + void processDummy(Colls const&) {} PROCESS_SWITCH(TreeWriterTPCTOF, processDummy, "Dummy function", false); From c450b27f7f969789e1949e4625e3049bbbb95b87 Mon Sep 17 00:00:00 2001 From: Ana Marin Date: Tue, 24 Jun 2025 15:46:26 +0200 Subject: [PATCH 2/4] fixing macos errors --- Common/TableProducer/PID/pidTPCBase.cxx | 1 - Common/TableProducer/PID/pidTPCBase.h | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/Common/TableProducer/PID/pidTPCBase.cxx b/Common/TableProducer/PID/pidTPCBase.cxx index 3ba7e881eff..2d7525fbc3e 100644 --- a/Common/TableProducer/PID/pidTPCBase.cxx +++ b/Common/TableProducer/PID/pidTPCBase.cxx @@ -127,7 +127,6 @@ struct DeDxCorrection { occupancy = 0; } - float signedP = trk.sign() * trk.tpcInnerParam(); float fTPCSignal = trk.tpcSignal(); float fNormMultTPC = multTPC / 11000.; diff --git a/Common/TableProducer/PID/pidTPCBase.h b/Common/TableProducer/PID/pidTPCBase.h index f6017758519..5fe8bb4726b 100644 --- a/Common/TableProducer/PID/pidTPCBase.h +++ b/Common/TableProducer/PID/pidTPCBase.h @@ -94,8 +94,8 @@ typedef struct Str_dEdx_correction { for (int i = 0; i < fMatrix.GetNrows(); i++) { for (int j = 0; j < fMatrix.GetNcols(); j++) { double param = fMatrix(i, j); - double value1 = i > vec1.size() ? 0 : vec1[i]; - double value2 = j > vec2.size() ? 0 : vec2[j]; + double value1 = i > static_cast(vec1.size()) ? 0 : vec1[i]; + double value2 = j > static_cast(vec2.size()) ? 0 : vec2[j]; result += param * value1 * value2; } } From b61a718324e7aea749667f711e70fd4b923df80c Mon Sep 17 00:00:00 2001 From: Ana Marin Date: Tue, 24 Jun 2025 15:56:10 +0200 Subject: [PATCH 3/4] fixing some errors --- Common/TableProducer/PID/pidTPCBase.cxx | 2 +- DPG/Tasks/TPC/tpcSkimsTableCreator.cxx | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Common/TableProducer/PID/pidTPCBase.cxx b/Common/TableProducer/PID/pidTPCBase.cxx index 2d7525fbc3e..e28f4a55b4d 100644 --- a/Common/TableProducer/PID/pidTPCBase.cxx +++ b/Common/TableProducer/PID/pidTPCBase.cxx @@ -23,8 +23,8 @@ #include "TableHelper.h" #include "pidTPCBase.h" -#include "Common/DataModel/FT0Corrected.h" #include "Common/CCDB/ctpRateFetcher.h" +#include "Common/DataModel/FT0Corrected.h" #include "CCDB/BasicCCDBManager.h" #include "Framework/AnalysisTask.h" diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index 3c3747bd31e..592cc24e2b1 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -855,11 +855,11 @@ struct TreeWriterTPCTOF { aod::pidTOFFullPr, aod::pidTOFFullDe, aod::pidTOFFullTr, aod::TrackSelection>; using TrksWithDEdxCorrection = soa::Join; + aod::pidTPCFullEl, aod::pidTPCFullPi, aod::pidTPCFullKa, + aod::pidTPCFullPr, aod::pidTPCFullDe, aod::pidTPCFullTr, + aod::pidTOFFullEl, aod::pidTOFFullPi, aod::pidTOFFullKa, + aod::pidTOFFullPr, aod::pidTOFFullDe, aod::pidTOFFullTr, + aod::TrackSelection, aod::DEdxsCorrected>; using Colls = soa::Join; using MyBCTable = soa::Join; From 528c429c8114860e8ccaaa053012e273a9075e20 Mon Sep 17 00:00:00 2001 From: Ana Marin Date: Tue, 24 Jun 2025 18:24:19 +0200 Subject: [PATCH 4/4] fixing an error --- DPG/Tasks/TPC/tpcSkimsTableCreator.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx index 592cc24e2b1..2a2d644e7c4 100644 --- a/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx +++ b/DPG/Tasks/TPC/tpcSkimsTableCreator.cxx @@ -1282,7 +1282,7 @@ struct TreeWriterTPCTOF { if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() <= maxMomTPCOnlyKa && std::abs(trk.tpcNSigmaKa()) < nSigmaTPCOnlyKa && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); } else if (trk.tpcInnerParam() < maxMomHardCutOnlyKa && trk.tpcInnerParam() > maxMomTPCOnlyKa && std::abs(trk.tofNSigmaKa()) < nSigmaTOF_TPCTOF_Ka && std::abs(trk.tpcNSigmaKa()) < nSigmaTPC_TPCTOF_Ka && downsampleTsalisCharged(trk.pt(), downsamplingTsalisKaons, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Kaon])) { - (trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); + fillSkimmedTPCTOFTableWithdEdxTrkQA(trk, trackQA, existTrkQA, collision, trk.tpcNSigmaKa(), trk.tofNSigmaKa(), trk.itsNSigmaKa(), trk.tpcExpSignalKa(trk.tpcSignal()), o2::track::PID::Kaon, runnumber, dwnSmplFactor_Ka, hadronicRate); } /// Fill tree pions if (trk.tpcInnerParam() <= maxMomTPCOnlyPi && std::abs(trk.tpcNSigmaPi()) < nSigmaTPCOnlyPi && downsampleTsalisCharged(trk.pt(), downsamplingTsalisPions, sqrtSNN, o2::track::pid_constants::sMasses[o2::track::PID::Pion])) {