From e2930f1e10b3d001b3d38bf8ffca1c546fc7acdd Mon Sep 17 00:00:00 2001 From: Preet Pati Date: Wed, 24 Jun 2026 09:41:29 +0200 Subject: [PATCH] Addition of PID --- .../TwoParticleCorrelations/Tasks/corrFit.cxx | 434 +++++++++++++++--- 1 file changed, 373 insertions(+), 61 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/corrFit.cxx b/PWGCF/TwoParticleCorrelations/Tasks/corrFit.cxx index dc343428a4c..b78d6f28463 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/corrFit.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/corrFit.cxx @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file corrFit.cxx -/// \brief Ultra long range correlation using forward FIT detectors and TPC, with foxus on multiplicity dependence +/// \brief Ultra long range correlation using forward FIT detectors and TPC, with focus on multiplicity dependence /// \author Thor Jensen (thor.kjaersgaard.jensen@cern.ch) #include "PWGCF/Core/CorrelationContainer.h" @@ -21,6 +21,9 @@ #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" #include @@ -61,9 +64,10 @@ using namespace o2::aod::rctsel; using namespace constants::math; #define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; +static constexpr float LongArrayFloat[3][20] = {{1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}, {2.1, 2.2, 2.3, -2.1, -2.2, -2.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}, {3.1, 3.2, 3.3, -3.1, -3.2, -3.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}}; struct CorrFit { - + o2::aod::ITSResponse itsResponse; Service ccdb; O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, true, "Use additional event cut on mult correlations") @@ -76,6 +80,7 @@ struct CorrFit { O2_DEFINE_CONFIGURABLE(cfgMinMultForCorrelations, int, 0, "minimum multiplicity for correlations") O2_DEFINE_CONFIGURABLE(cfgMaxMultForCorrelations, int, 20, "maximum multiplicity for correlations") O2_DEFINE_CONFIGURABLE(cfgRefMultiplicity, bool, false, "Use multiplicity of reference tracks for multiplicity correlation cut instead of Nch") + Configurable> cfgRunRemoveList{"cfgRunRemoveList", std::vector{-1}, "excluded run numbers"}; struct : ConfigurableGroup{ O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.2f, "minimum accepted track pT") @@ -85,9 +90,7 @@ struct CorrFit { O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 50.0f, "minimum TPC clusters") O2_DEFINE_CONFIGURABLE(cfgCutTPCCrossedRows, float, 70.0f, "minimum TPC crossed rows") O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum ITS clusters") - O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "max DCA to vertex z") - - } cfgTrackCuts; + O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "max DCA to vertex z")} cfgTrackCuts; struct : ConfigurableGroup{ O2_DEFINE_CONFIGURABLE(cfgEvSelkNoSameBunchPileup, bool, false, "rejects collisions which are associated with the same found-by-T0 bunch crossing") @@ -103,17 +106,7 @@ struct CorrFit { O2_DEFINE_CONFIGURABLE(cfgEvSelV0AT0ACut, bool, false, "V0A T0A 5 sigma cut") O2_DEFINE_CONFIGURABLE(cfgEvSelOccupancy, bool, false, "Occupancy cut") O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 2000, "High cut on TPC occupancy") - O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy") - - } cfgEventSelection; - - struct : ConfigurableGroup{ - O2_DEFINE_CONFIGURABLE(cfgSystematicsVariation, bool, false, "Enable systematics variation for track cuts") - O2_DEFINE_CONFIGURABLE(cfgSystematicsCutChi2prTPCcls, float, 3.0f, "max chi2 per TPC clusters for systematics variation") - O2_DEFINE_CONFIGURABLE(cfgSystematicsCutTPCclu, float, 40.0f, "minimum TPC clusters for systematics variation") - O2_DEFINE_CONFIGURABLE(cfgSystematicsCutTPCCrossedRows, float, 60.0f, "minimum TPC crossed rows for systematics variation") - O2_DEFINE_CONFIGURABLE(cfgSystematicsCutITSclu, float, 4.0f, "minimum ITS clusters for systematics variation") - O2_DEFINE_CONFIGURABLE(cfgSystematicsCutDCAz, float, 1.5f, "max DCA to vertex z for systematics variation")} cfgSystematics; + O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy")} cfgEventSelection; O2_DEFINE_CONFIGURABLE(cfgMinMixEventNum, int, 5, "Minimum number of events to mix") O2_DEFINE_CONFIGURABLE(cfgMergingCut, float, 0.02, "Merging cut on track merge") @@ -159,6 +152,16 @@ struct CorrFit { TF1* fT0AV0ASigma = nullptr; } cfgFuncParas; + struct : ConfigurableGroup { + O2_DEFINE_CONFIGURABLE(cfgTofPtCut, float, 0.4f, "Minimum pt to use TOF N-sigma") + O2_DEFINE_CONFIGURABLE(cfgUseItsPID, bool, false, "Use ITS PID for particle identification") + O2_DEFINE_CONFIGURABLE(cfgPIDParticle, int, 0, "1 = pion, 2 = kaon, 3 = proton, 4 = kshort, 5 = lambda, 6 = phi, 0 for no PID") + O2_DEFINE_CONFIGURABLE(cfgGetNsigmaQA, bool, false, "Get QA histograms for selection of pions, kaons, and protons") + O2_DEFINE_CONFIGURABLE(cfgGetdEdx, bool, false, "Get dEdx histograms for pions, kaons, and protons") + O2_DEFINE_CONFIGURABLE(cfgPIDUseRejection, bool, true, "True: use exclusion exclusion criteria for PID determination, false: don't use exclusion") + Configurable> nSigmas{"nSigmas", {LongArrayFloat[0], 6, 3, {"UpCut_pi", "UpCut_ka", "UpCut_pr", "LowCut_pi", "LowCut_ka", "LowCut_pr"}, {"TPC", "TOF", "ITS"}}, "Labeled array for n-sigma values for TPC, TOF, ITS for pions, kaons, protons (positive and negative)"}; + } cfgPIDConfgs; + Configurable cfgCutFV0{"cfgCutFV0", 50., "FV0A threshold"}; Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; @@ -170,7 +173,7 @@ struct CorrFit { ConfigurableAxis axisMult{"axisMult", {10, 0, 100}, "multiplicity axis for histograms"}; ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"}; ConfigurableAxis axisPhi{"axisPhi", {72, 0.0, constants::math::TwoPI}, "phi axis for histograms"}; - ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt axis for histograms"}; + ConfigurableAxis axisPtFiner{"axisPtFiner", {98, 0.2, 10.0}, "pt axis for histograms"}; ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"}; ConfigurableAxis axisDeltaEtaTpcFt0a{"axisDeltaEtaTpcFt0a", {32, -5.8, -2.6}, "delta eta axis, -5.8~-2.6 for TPC-FT0A,"}; ConfigurableAxis axisDeltaEtaTpcFt0c{"axisDeltaEtaTpcFt0c", {32, 1.2, 4.2}, "delta eta axis, 1.2~4.2 for TPC-FT0C"}; @@ -181,6 +184,7 @@ struct CorrFit { ConfigurableAxis axisVtxMix{"axisVtxMix", {VARIABLE_WIDTH, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, "vertex axis for mixed event histograms"}; ConfigurableAxis axisMultMix{"axisMultMix", {VARIABLE_WIDTH, 0, 10, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260}, "multiplicity / centrality axis for mixed event histograms"}; ConfigurableAxis axisSample{"axisSample", {cfgSampleSize, 0, cfgSampleSize}, "sample axis for histograms"}; + ConfigurableAxis axisParticle{"axisParticle", {4, 0, 4}, "particle axis for correlation container"}; ConfigurableAxis axisVertexEfficiency{"axisVertexEfficiency", {10, -10, 10}, "vertex axis for efficiency histograms"}; ConfigurableAxis axisEtaEfficiency{"axisEtaEfficiency", {20, -1.0, 1.0}, "eta axis for efficiency histograms"}; @@ -188,6 +192,12 @@ struct CorrFit { ConfigurableAxis axisAmplitudeFt0a{"axisAmplitudeFt0a", {5000, 0, 1000}, "FT0A amplitude"}; ConfigurableAxis axisChannelFt0aAxis{"axisChannelFt0aAxis", {96, 0.0, 96.0}, "FT0A channel"}; + ConfigurableAxis axisNsigmaTPC{"axisNsigmaTPC", {80, -5, 5}, "nsigmaTPC axis"}; + ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {80, -5, 5}, "nsigmaTOF axis"}; + ConfigurableAxis axisNsigmaITS{"axisNsigmaITS", {80, -5, 5}, "nsigmaITS axis"}; + ConfigurableAxis axisTpcSignal{"axisTpcSignal", {250, 0, 250}, "dEdx axis for TPC"}; + ConfigurableAxis axisSigma{"axisSigma", {200, 0, 20}, "sigma axis for TPC"}; + Configurable cfgGainEqPath{"cfgGainEqPath", "Analysis/EventPlane/GainEq", "CCDB path for gain equalization constants"}; Configurable cfgCorrLevel{"cfgCorrLevel", 0, "calibration step: 0 = no corr, 1 = gain corr"}; ConfigurableAxis cfgaxisFITamp{"cfgaxisFITamp", {1000, 0, 5000}, ""}; @@ -195,10 +205,10 @@ struct CorrFit { AxisSpec axisChID = {220, 0, 220}; Filter collisionFilter = (nabs(aod::collision::posZ) < cfgZVtxCut); - Filter trackFilter = (nabs(aod::track::eta) < cfgTrackCuts.cfgEtaCut) && (cfgTrackCuts.cfgPtCutMin < aod::track::pt) && (cfgTrackCuts.cfgPtCutMax > aod::track::pt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)) && (aod::track::tpcChi2NCl < cfgTrackCuts.cfgCutChi2prTPCcls) && (aod::track::dcaZ < cfgTrackCuts.cfgCutDCAz); + Filter trackFilter = (nabs(aod::track::eta) < cfgTrackCuts.cfgEtaCut) && (cfgTrackCuts.cfgPtCutMin < aod::track::pt) && (cfgTrackCuts.cfgPtCutMax > aod::track::pt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)); using FilteredCollisions = soa::Filtered>; - using FilteredTracks = soa::Filtered>; + using FilteredTracks = soa::Filtered>; // FT0 geometry o2::ft0::Geometry ft0Det; @@ -237,8 +247,35 @@ struct CorrFit { kFT0C = 1 }; + enum PIDIndex { + kCharged = 0, + kPions, + kKaons, + kProtons, + kK0, + kLambda, + kPhi + }; + enum PiKpArrayIndex { + iPionUp = 0, + iKaonUp, + iProtonUp, + iPionLow, + iKaonLow, + iProtonLow + }; + enum DetectorType { + kTPC = 0, + kTOF, + kITS + }; + RCTFlagsChecker rctChecker{"CBT"}; + std::array tofNsigmaCut; + std::array itsNsigmaCut; + std::array tpcNsigmaCut; + void init(InitContext&) { @@ -249,6 +286,27 @@ struct CorrFit { LOGF(info, "Starting init"); + tpcNsigmaCut[iPionUp] = cfgPIDConfgs.nSigmas->getData()[iPionUp][kTPC]; + tpcNsigmaCut[iKaonUp] = cfgPIDConfgs.nSigmas->getData()[iKaonUp][kTPC]; + tpcNsigmaCut[iProtonUp] = cfgPIDConfgs.nSigmas->getData()[iProtonUp][kTPC]; + tpcNsigmaCut[iPionLow] = cfgPIDConfgs.nSigmas->getData()[iPionLow][kTPC]; + tpcNsigmaCut[iKaonLow] = cfgPIDConfgs.nSigmas->getData()[iKaonLow][kTPC]; + tpcNsigmaCut[iProtonLow] = cfgPIDConfgs.nSigmas->getData()[iProtonLow][kTPC]; + + tofNsigmaCut[iPionUp] = cfgPIDConfgs.nSigmas->getData()[iPionUp][kTOF]; + tofNsigmaCut[iKaonUp] = cfgPIDConfgs.nSigmas->getData()[iKaonUp][kTOF]; + tofNsigmaCut[iProtonUp] = cfgPIDConfgs.nSigmas->getData()[iProtonUp][kTOF]; + tofNsigmaCut[iPionLow] = cfgPIDConfgs.nSigmas->getData()[iPionLow][kTOF]; + tofNsigmaCut[iKaonLow] = cfgPIDConfgs.nSigmas->getData()[iKaonLow][kTOF]; + tofNsigmaCut[iProtonLow] = cfgPIDConfgs.nSigmas->getData()[iProtonLow][kTOF]; + + itsNsigmaCut[iPionUp] = cfgPIDConfgs.nSigmas->getData()[iPionUp][kITS]; + itsNsigmaCut[iKaonUp] = cfgPIDConfgs.nSigmas->getData()[iKaonUp][kITS]; + itsNsigmaCut[iProtonUp] = cfgPIDConfgs.nSigmas->getData()[iProtonUp][kITS]; + itsNsigmaCut[iPionLow] = cfgPIDConfgs.nSigmas->getData()[iPionLow][kITS]; + itsNsigmaCut[iKaonLow] = cfgPIDConfgs.nSigmas->getData()[iKaonLow][kITS]; + itsNsigmaCut[iProtonLow] = cfgPIDConfgs.nSigmas->getData()[iProtonLow][kITS]; + if (doprocessSameFt0aFt0c || doprocessSameTpcFt0a || doprocessSameTpcFt0c || doprocessSameTPC) { registry.add("hEventCountRct", "Number of Event;; Count", {HistType::kTH1D, {{2, 0, 2}}}); registry.get(HIST("hEventCountRct"))->GetXaxis()->SetBinLabel(1, "rct fail"); @@ -289,17 +347,41 @@ struct CorrFit { registry.add("Phi", "Phi", {HistType::kTH1D, {axisPhi}}); registry.add("Eta", "Eta", {HistType::kTH1D, {axisEta}}); registry.add("EtaCorrected", "EtaCorrected", {HistType::kTH1D, {axisEta}}); - registry.add("pT", "pT", {HistType::kTH1D, {axisPtTrigger}}); - registry.add("pTCorrected", "pTCorrected", {HistType::kTH1D, {axisPtTrigger}}); + registry.add("hParticleCounts", "hParticleCounts", {HistType::kTH1D, {axisParticle}}); + registry.add("hParticleSelected", "hParticleSelected", {HistType::kTH2D, {axisParticle, axisPtTrigger}}); + registry.add("pTFiner", "pTFiner", {HistType::kTH1D, {axisPtFiner}}); + registry.add("pTFinerCorrected", "pTFinerCorrected", {HistType::kTH1D, {axisPtFiner}}); registry.add("Nch", "N_{ch}", {HistType::kTH1D, {axisMult}}); registry.add("zVtx", "zVtx", {HistType::kTH1D, {axisVertex}}); if (doprocessSameFt0aFt0c || doprocessSameTpcFt0a || doprocessSameTpcFt0c) { registry.add("FT0Amp", "", {HistType::kTH2F, {axisChID, axisFit}}); registry.add("FT0AmpCorrect", "", {HistType::kTH2F, {axisChID, axisFit}}); } - if (cfgQaCheck) { - registry.add("Nch_corrected", "N_{ch} corrected", {HistType::kTH1D, {axisMult}}); + if (cfgPIDConfgs.cfgGetNsigmaQA) { + if (!cfgPIDConfgs.cfgUseItsPID) { + registry.add("TofTpcNsigma_before", "", {HistType::kTHnSparseD, {{axisNsigmaTPC, axisNsigmaTOF, axisPtTrigger}}}); + registry.add("TofTpcNsigma_after", "", {HistType::kTHnSparseD, {{axisNsigmaTPC, axisNsigmaTOF, axisPtTrigger}}}); + } // TPC-TOF PID QA hists + if (cfgPIDConfgs.cfgUseItsPID) { + registry.add("TofItsNsigma_before", "", {HistType::kTHnSparseD, {{axisNsigmaITS, axisNsigmaTOF, axisPtTrigger}}}); + registry.add("TofItsNsigma_after", "", {HistType::kTHnSparseD, {{axisNsigmaITS, axisNsigmaTOF, axisPtTrigger}}}); + } // ITS-TOF PID QA hists + } // end of PID QA hists + + if (cfgPIDConfgs.cfgGetdEdx) { + registry.add("TpcdEdx_ptwise_beforeCut", "", {HistType::kTHnSparseD, {{axisPtTrigger, axisTpcSignal, axisNsigmaTOF}}}); + registry.add("ExpTpcdEdx_ptwise_beforeCut", "", {HistType::kTHnSparseD, {{axisPtTrigger, axisTpcSignal, axisNsigmaTOF}}}); + registry.add("ExpSigma_ptwise_beforeCut", "", {HistType::kTHnSparseD, {{axisPtTrigger, axisSigma, axisNsigmaTOF}}}); + + registry.add("TpcdEdx_ptwise_afterCut", "", {HistType::kTHnSparseD, {{axisPtTrigger, axisTpcSignal, axisNsigmaTOF}}}); + registry.add("ExpTpcdEdx_ptwise_afterCut", "", {HistType::kTHnSparseD, {{axisPtTrigger, axisTpcSignal, axisNsigmaTOF}}}); + registry.add("ExpSigma_ptwise_afterCut", "", {HistType::kTHnSparseD, {{axisPtTrigger, axisSigma, axisNsigmaTOF}}}); } + + } // end of single particle distribution histograms + + if (cfgQaCheck) { + registry.add("Nch_corrected", "N_{ch} corrected", {HistType::kTH1D, {axisMult}}); } if (doprocessSameTpcFt0a) { @@ -593,13 +675,7 @@ struct CorrFit { template bool trackSelected(TTrack const& track) { - return ((track.tpcNClsFound() >= cfgTrackCuts.cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgTrackCuts.cfgCutTPCCrossedRows) && (track.itsNCls() >= cfgTrackCuts.cfgCutITSclu)); - } - - template - bool trackSelectedSystematics(TTrack const& track) - { - return ((track.tpcNClsFound() >= cfgSystematics.cfgSystematicsCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgSystematics.cfgSystematicsCutTPCCrossedRows) && (track.itsNCls() >= cfgSystematics.cfgSystematicsCutITSclu)); + return ((track.tpcNClsFound() >= cfgTrackCuts.cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgTrackCuts.cfgCutTPCCrossedRows) && (track.itsNCls() >= cfgTrackCuts.cfgCutITSclu) && (track.tpcChi2NCl() < cfgTrackCuts.cfgCutChi2prTPCcls) && (track.dcaZ() < cfgTrackCuts.cfgCutDCAz)); } template @@ -679,6 +755,65 @@ struct CorrFit { } } + template + int getNsigmaPID(TTrack track, bool fillYields) + { + // Computing Nsigma arrays for pion, kaon, and protons + std::array nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()}; + std::array nSigmaTOF = {track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; + std::array nSigmaITS = {itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track)}; + int pid = -1; // -1 = not identified, 1 = pion, 2 = kaon, 3 = proton + + std::array nSigmaToUse = cfgPIDConfgs.cfgUseItsPID ? nSigmaITS : nSigmaTPC; // Choose which nSigma to use: TPC or ITS + std::array detectorNsigmaCut = cfgPIDConfgs.cfgUseItsPID ? itsNsigmaCut : tpcNsigmaCut; // Choose which nSigma to use: TPC or ITS + + bool isPion, isKaon, isProton; + bool isDetectedPion = nSigmaToUse[iPionUp] < detectorNsigmaCut[iPionUp] && nSigmaToUse[iPionUp] > detectorNsigmaCut[iPionLow]; + bool isDetectedKaon = nSigmaToUse[iKaonUp] < detectorNsigmaCut[iKaonUp] && nSigmaToUse[iKaonUp] > detectorNsigmaCut[iKaonLow]; + bool isDetectedProton = nSigmaToUse[iProtonUp] < detectorNsigmaCut[iProtonUp] && nSigmaToUse[iProtonUp] > detectorNsigmaCut[iProtonLow]; + + bool isTofPion = nSigmaTOF[iPionUp] < tofNsigmaCut[iPionUp] && nSigmaTOF[iPionUp] > tofNsigmaCut[iPionLow]; + bool isTofKaon = nSigmaTOF[iKaonUp] < tofNsigmaCut[iKaonUp] && nSigmaTOF[iKaonUp] > tofNsigmaCut[iKaonLow]; + bool isTofProton = nSigmaTOF[iProtonUp] < tofNsigmaCut[iProtonUp] && nSigmaTOF[iProtonUp] > tofNsigmaCut[iProtonLow]; + + if (track.pt() > cfgPIDConfgs.cfgTofPtCut && !track.hasTOF()) { + return -1; + } else if (track.pt() > cfgPIDConfgs.cfgTofPtCut && track.hasTOF()) { + isPion = isTofPion && isDetectedPion; + isKaon = isTofKaon && isDetectedKaon; + isProton = isTofProton && isDetectedProton; + } else { + isPion = isDetectedPion; + isKaon = isDetectedKaon; + isProton = isDetectedProton; + } + + if (cfgPIDConfgs.cfgPIDUseRejection && ((isPion && isKaon) || (isPion && isProton) || (isKaon && isProton))) { + return -1; // more than one particle satisfy the criteria + } + + if (fillYields) + registry.fill(HIST("hParticleCounts"), kCharged); + + if (isPion) { + pid = kPions; + if (fillYields) + registry.fill(HIST("hParticleCounts"), kPions); + } else if (isKaon) { + pid = kKaons; + if (fillYields) + registry.fill(HIST("hParticleCounts"), kKaons); + } else if (isProton) { + pid = kProtons; + if (fillYields) + registry.fill(HIST("hParticleCounts"), kProtons); + } else { + return -1; // no particle satisfies the criteria + } + + return pid; // -1 = not identified, 1 = pion, 2 = kaon, 3 = proton + } + void loadCorrection(uint64_t timestamp) { if (correctionsLoaded) { @@ -792,6 +927,138 @@ struct CorrFit { multiplicity = nTracksCorrected; } + template + void fillYield(TCollision collision, TTracks tracks) // function to fill the yield and etaphi histograms. + { + + float weff1 = 1.0; + float zvtx = collision.posZ(); + + for (auto const& track1 : tracks) { + + if (!trackSelected(track1)) { + continue; + } + if (cfgPIDConfgs.cfgPIDParticle && getNsigmaPID(track1, true) != cfgPIDConfgs.cfgPIDParticle) + continue; // if PID is selected, check if the track has the right PID + + if (!getEfficiencyCorrection(weff1, track1.eta(), track1.pt(), zvtx)) { + continue; + } + + registry.fill(HIST("Phi"), RecoDecay::constrainAngle(track1.phi(), 0.0)); + registry.fill(HIST("Eta"), track1.eta()); + registry.fill(HIST("EtaCorrected"), track1.eta(), weff1); + registry.fill(HIST("pTFiner"), track1.pt()); + registry.fill(HIST("pTFinerCorrected"), track1.pt(), weff1); + registry.fill(HIST("hParticleSelected"), cfgPIDConfgs.cfgPIDParticle.value, track1.pt()); + } + } + + template + void fillNsigmaBeforeCut(TTrack track1, Int_t pid) // function to fill the QA before Nsigma selection + { + switch (pid) { + case kPions: // For Pions + if (!cfgPIDConfgs.cfgUseItsPID) { + if (cfgPIDConfgs.cfgGetNsigmaQA) + registry.fill(HIST("TofTpcNsigma_before"), track1.tpcNSigmaPi(), track1.tofNSigmaPi(), track1.pt()); + if (cfgPIDConfgs.cfgGetdEdx) { + double tpcExpSignalPi = track1.tpcSignal() - (track1.tpcNSigmaPi() * track1.tpcExpSigmaPi()); + + registry.fill(HIST("TpcdEdx_ptwise_beforeCut"), track1.pt(), track1.tpcSignal(), track1.tofNSigmaPi()); + registry.fill(HIST("ExpTpcdEdx_ptwise_beforeCut"), track1.pt(), tpcExpSignalPi, track1.tofNSigmaPi()); + registry.fill(HIST("ExpSigma_ptwise_beforeCut"), track1.pt(), track1.tpcExpSigmaPi(), track1.tofNSigmaPi()); + } + } + if (cfgPIDConfgs.cfgGetNsigmaQA && cfgPIDConfgs.cfgUseItsPID) + registry.fill(HIST("TofItsNsigma_before"), itsResponse.nSigmaITS(track1), track1.tofNSigmaPi(), track1.pt()); + break; + case kKaons: // For Kaons + if (!cfgPIDConfgs.cfgUseItsPID) { + if (cfgPIDConfgs.cfgGetNsigmaQA) + registry.fill(HIST("TofTpcNsigma_before"), track1.tpcNSigmaKa(), track1.tofNSigmaKa(), track1.pt()); + if (cfgPIDConfgs.cfgGetdEdx) { + double tpcExpSignalKa = track1.tpcSignal() - (track1.tpcNSigmaKa() * track1.tpcExpSigmaKa()); + + registry.fill(HIST("TpcdEdx_ptwise_beforeCut"), track1.pt(), track1.tpcSignal(), track1.tofNSigmaKa()); + registry.fill(HIST("ExpTpcdEdx_ptwise_beforeCut"), track1.pt(), tpcExpSignalKa, track1.tofNSigmaKa()); + registry.fill(HIST("ExpSigma_ptwise_beforeCut"), track1.pt(), track1.tpcExpSigmaKa(), track1.tofNSigmaKa()); + } + } + if (cfgPIDConfgs.cfgGetNsigmaQA && cfgPIDConfgs.cfgUseItsPID) + registry.fill(HIST("TofItsNsigma_before"), itsResponse.nSigmaITS(track1), track1.tofNSigmaKa(), track1.pt()); + break; + case kProtons: // For Protons + if (!cfgPIDConfgs.cfgUseItsPID) { + if (cfgPIDConfgs.cfgGetNsigmaQA) + registry.fill(HIST("TofTpcNsigma_before"), track1.tpcNSigmaPr(), track1.tofNSigmaPr(), track1.pt()); + if (cfgPIDConfgs.cfgGetdEdx) { + double tpcExpSignalPr = track1.tpcSignal() - (track1.tpcNSigmaPr() * track1.tpcExpSigmaPr()); + + registry.fill(HIST("TpcdEdx_ptwise_beforeCut"), track1.pt(), track1.tpcSignal(), track1.tofNSigmaPr()); + registry.fill(HIST("ExpTpcdEdx_ptwise_beforeCut"), track1.pt(), tpcExpSignalPr, track1.tofNSigmaPr()); + registry.fill(HIST("ExpSigma_ptwise_beforeCut"), track1.pt(), track1.tpcExpSigmaPr(), track1.tofNSigmaPr()); + } + } + if (cfgPIDConfgs.cfgGetNsigmaQA && cfgPIDConfgs.cfgUseItsPID) + registry.fill(HIST("TofItsNsigma_before"), itsResponse.nSigmaITS(track1), track1.tofNSigmaPr(), track1.pt()); + break; + } + } // end of fillNsigmaBeforeCut + + template + void fillNsigmaAfterCut(TTrack track1, Int_t pid) // function to fill the QA after Nsigma selection + { + switch (pid) { + case kPions: // For Pions + if (!cfgPIDConfgs.cfgUseItsPID) { + if (cfgPIDConfgs.cfgGetdEdx) { + double tpcExpSignalPi = track1.tpcSignal() - (track1.tpcNSigmaPi() * track1.tpcExpSigmaPi()); + + registry.fill(HIST("TpcdEdx_ptwise_afterCut"), track1.pt(), track1.tpcSignal(), track1.tofNSigmaPi()); + registry.fill(HIST("ExpTpcdEdx_ptwise_afterCut"), track1.pt(), tpcExpSignalPi, track1.tofNSigmaPi()); + registry.fill(HIST("ExpSigma_ptwise_afterCut"), track1.pt(), track1.tpcExpSigmaPi(), track1.tofNSigmaPi()); + } + if (cfgPIDConfgs.cfgGetNsigmaQA) + registry.fill(HIST("TofTpcNsigma_after"), track1.tpcNSigmaPi(), track1.tofNSigmaPi(), track1.pt()); + } + if (cfgPIDConfgs.cfgUseItsPID) + registry.fill(HIST("TofItsNsigma_after"), itsResponse.nSigmaITS(track1), track1.tofNSigmaPi(), track1.pt()); + break; + case kKaons: // For Kaons + if (!cfgPIDConfgs.cfgUseItsPID) { + if (cfgPIDConfgs.cfgGetdEdx) { + double tpcExpSignalKa = track1.tpcSignal() - (track1.tpcNSigmaKa() * track1.tpcExpSigmaKa()); + + registry.fill(HIST("TpcdEdx_ptwise_afterCut"), track1.pt(), track1.tpcSignal(), track1.tofNSigmaKa()); + registry.fill(HIST("ExpTpcdEdx_ptwise_afterCut"), track1.pt(), tpcExpSignalKa, track1.tofNSigmaKa()); + registry.fill(HIST("ExpSigma_ptwise_afterCut"), track1.pt(), track1.tpcExpSigmaKa(), track1.tofNSigmaKa()); + } + if (cfgPIDConfgs.cfgGetNsigmaQA) + registry.fill(HIST("TofTpcNsigma_after"), track1.tpcNSigmaKa(), track1.tofNSigmaKa(), track1.pt()); + } + if (cfgPIDConfgs.cfgUseItsPID) + registry.fill(HIST("TofItsNsigma_after"), itsResponse.nSigmaITS(track1), track1.tofNSigmaKa(), track1.pt()); + break; + case kProtons: // For Protons + if (!cfgPIDConfgs.cfgUseItsPID) { + if (cfgPIDConfgs.cfgGetdEdx) { + double tpcExpSignalPr = track1.tpcSignal() - (track1.tpcNSigmaPr() * track1.tpcExpSigmaPr()); + + registry.fill(HIST("TpcdEdx_ptwise_afterCut"), track1.pt(), track1.tpcSignal(), track1.tofNSigmaPr()); + registry.fill(HIST("ExpTpcdEdx_ptwise_afterCut"), track1.pt(), tpcExpSignalPr, track1.tofNSigmaPr()); + registry.fill(HIST("ExpSigma_ptwise_afterCut"), track1.pt(), track1.tpcExpSigmaPr(), track1.tofNSigmaPr()); + } + if (cfgPIDConfgs.cfgGetNsigmaQA) + registry.fill(HIST("TofTpcNsigma_after"), track1.tpcNSigmaPr(), track1.tofNSigmaPr(), track1.pt()); + } + if (cfgPIDConfgs.cfgUseItsPID && cfgPIDConfgs.cfgGetNsigmaQA) + registry.fill(HIST("TofItsNsigma_after"), itsResponse.nSigmaITS(track1), track1.tofNSigmaPr(), track1.pt()); + break; + } // end of switch + } + template void fillCorrelationsTPCFT0(TTracks tracks1, TFT0s const& ft0, float posZ, int system, int multiplicity, int corType, float eventWeight) // function to fill the Output functions (sparse) and the delta eta and delta phi histograms { @@ -805,10 +1072,14 @@ struct CorrFit { if (!trackSelected(track1)) continue; - if (cfgSystematics.cfgSystematicsVariation) { - if (!trackSelectedSystematics(track1)) - continue; - } + if (cfgPIDConfgs.cfgPIDParticle > kCharged) + fillNsigmaBeforeCut(track1, cfgPIDConfgs.cfgPIDParticle); + + if (cfgPIDConfgs.cfgPIDParticle && getNsigmaPID(track1, false) != cfgPIDConfgs.cfgPIDParticle) + continue; // if PID is selected, check if the track has the right PID + + if (cfgPIDConfgs.cfgPIDParticle > kCharged) + fillNsigmaAfterCut(track1, cfgPIDConfgs.cfgPIDParticle); if (!getEfficiencyCorrection(triggerWeight, track1.pt(), track1.eta(), posZ)) continue; @@ -875,32 +1146,6 @@ struct CorrFit { } } - template - void fillYield(TCollision collision, TTracks tracks) // function to fill the yield and etaphi histograms. - { - - float weff1 = 1.0; - float zvtx = collision.posZ(); - registry.fill(HIST("zVtx"), zvtx); - registry.fill(HIST("Nch"), tracks.size()); - - for (auto const& track1 : tracks) { - - if (!trackSelected(track1)) { - continue; - } - if (!getEfficiencyCorrection_Nch(weff1, track1.pt())) { - continue; - } - - registry.fill(HIST("Phi"), RecoDecay::constrainAngle(track1.phi(), 0.0)); - registry.fill(HIST("Eta"), track1.eta()); - registry.fill(HIST("EtaCorrected"), track1.eta(), weff1); - registry.fill(HIST("pT"), track1.pt()); - registry.fill(HIST("pTCorrected"), track1.pt(), weff1); - } - } - template void fillCorrelationsFT0AFT0C(TFT0s const& ft0Col1, TFT0s const& ft0Col2, float posZ, int system, int multiplicity, float eventWeight) // function to fill the Output functions (sparse) and the delta eta and delta phi histograms { @@ -966,6 +1211,16 @@ struct CorrFit { if (!trackSelected(track1)) continue; + // Fill Nsigma QA + if (cfgPIDConfgs.cfgPIDParticle > kCharged) + fillNsigmaBeforeCut(track1, cfgPIDConfgs.cfgPIDParticle); + + if (cfgPIDConfgs.cfgPIDParticle && getNsigmaPID(track1, false) != cfgPIDConfgs.cfgPIDParticle) + continue; // if PID is selected, check if the track has the right PID + + if (cfgPIDConfgs.cfgPIDParticle > kCharged) + fillNsigmaAfterCut(track1, cfgPIDConfgs.cfgPIDParticle); + if (!getEfficiencyCorrection_Nch(triggerWeight, track1.pt())) continue; @@ -1033,6 +1288,17 @@ struct CorrFit { } } + bool isGoodRun(int runNumber) + { + for (const auto& ExcludedRun : cfgRunRemoveList.value) { + if (runNumber == ExcludedRun) { + return false; + } + } + + return true; + } + void processSameTpcFt0a(FilteredCollisions::iterator const& collision, FilteredTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&) { if (cfgQaCheck) { @@ -1047,6 +1313,12 @@ struct CorrFit { auto bc = collision.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } + if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) return; @@ -1118,6 +1390,12 @@ struct CorrFit { registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin auto bc = collision1.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } + loadAlignParam(bc.timestamp()); loadCorrection(bc.timestamp()); float eventWeight = 1.0f; @@ -1152,6 +1430,11 @@ struct CorrFit { return; auto bc = collision.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) return; @@ -1215,6 +1498,11 @@ struct CorrFit { registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin auto bc = collision1.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } loadAlignParam(bc.timestamp()); loadCorrection(bc.timestamp()); float eventWeight = 1.0f; @@ -1249,6 +1537,11 @@ struct CorrFit { return; auto bc = collision.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) return; @@ -1315,6 +1608,11 @@ struct CorrFit { continue; auto bc = collision1.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } loadAlignParam(bc.timestamp()); loadCorrection(bc.timestamp()); float eventWeight = 1.0f; @@ -1351,6 +1649,11 @@ struct CorrFit { return; auto bc = collision.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) return; @@ -1404,7 +1707,15 @@ struct CorrFit { continue; registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin - loadCorrection(collision1.bc_as().timestamp()); + + auto bc = collision1.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } + + loadCorrection(bc.timestamp()); double multiplicity = tracks1.size(); @@ -1417,6 +1728,7 @@ struct CorrFit { } PROCESS_SWITCH(CorrFit, processMixedTPC, "Process mixed events for TPC-TPC correlation", false); }; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{