|
| 1 | +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. |
| 2 | +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. |
| 3 | +// All rights not expressly granted are reserved. |
| 4 | +// |
| 5 | +// This software is distributed under the terms of the GNU General Public |
| 6 | +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". |
| 7 | +// |
| 8 | +// In applying this license CERN does not waive the privileges and immunities |
| 9 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 10 | +// or submit itself to any jurisdiction. |
| 11 | +// |
| 12 | +/// \brief optimization of particle identification for femtoscopic analysis. |
| 13 | +/// \author Sofia Tomassini |
| 14 | +/// \since July 2025 |
| 15 | + |
| 16 | +#include "Common/CCDB/ctpRateFetcher.h" |
| 17 | +#include "Common/Core/TrackSelection.h" |
| 18 | +#include "Common/DataModel/Centrality.h" |
| 19 | +#include "Common/DataModel/EventSelection.h" |
| 20 | +#include "Common/DataModel/Multiplicity.h" |
| 21 | +#include "Common/DataModel/PIDResponse.h" |
| 22 | +#include "Common/DataModel/PIDResponseITS.h" |
| 23 | +#include "Common/DataModel/TrackSelectionTables.h" |
| 24 | +#include "EventFiltering/Zorro.h" |
| 25 | +#include "EventFiltering/ZorroSummary.h" |
| 26 | + |
| 27 | +#include "CCDB/BasicCCDBManager.h" |
| 28 | +#include "CommonConstants/MathConstants.h" |
| 29 | +#include "Framework/ASoAHelpers.h" |
| 30 | +#include "Framework/AnalysisDataModel.h" |
| 31 | +#include "Framework/AnalysisTask.h" |
| 32 | +#include "Framework/Configurable.h" |
| 33 | +#include "Framework/HistogramRegistry.h" |
| 34 | +#include "Framework/StaticFor.h" |
| 35 | +#include "Framework/runDataProcessing.h" |
| 36 | + |
| 37 | +#include "Math/GenVector/Boost.h" |
| 38 | +#include "Math/Vector4D.h" |
| 39 | +#include "TMath.h" |
| 40 | +#include "TRandom3.h" |
| 41 | +#include <TParameter.h> |
| 42 | + |
| 43 | +#include "fairlogger/Logger.h" |
| 44 | + |
| 45 | +#include <iostream> |
| 46 | +#include <iterator> |
| 47 | +#include <memory> |
| 48 | +#include <string> |
| 49 | +#include <utility> |
| 50 | +#include <vector> |
| 51 | + |
| 52 | +using namespace o2; |
| 53 | +using namespace o2::framework; |
| 54 | +using namespace o2::framework::expressions; |
| 55 | +namespace o2::aod |
| 56 | +{ |
| 57 | +using SelectedTracks = soa::Join<aod::FullTracks, aod::TrackSelection, aod::TracksDCA, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTPCFullDe, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr, aod::pidTOFFullDe>; |
| 58 | +} |
| 59 | +struct PidOptimization { |
| 60 | + |
| 61 | + HistogramRegistry histos{"Histos"}; |
| 62 | + |
| 63 | + Configurable<bool> _removeSameBunchPileup{"removeSameBunchPileup", false, ""}; |
| 64 | + Configurable<bool> _requestGoodZvtxFT0vsPV{"requestGoodZvtxFT0vsPV", false, ""}; |
| 65 | + Configurable<bool> _requestVertexITSTPC{"requestVertexITSTPC", false, ""}; |
| 66 | + Configurable<int> _requestVertexTOForTRDmatched{"requestVertexTOFmatched", 0, "0 -> no selectio; 1 -> vertex is matched to TOF or TRD; 2 -> matched to both;"}; |
| 67 | + Configurable<bool> _requestNoCollInTimeRangeStandard{"requestNoCollInTimeRangeStandard", false, ""}; |
| 68 | + Configurable<std::pair<float, float>> _IRcut{"IRcut", std::pair<float, float>{0.f, 100.f}, "[min., max.] IR range to keep events within"}; |
| 69 | + Configurable<std::pair<int, int>> _OccupancyCut{"OccupancyCut", std::pair<int, int>{0, 10000}, "[min., max.] occupancy range to keep events within"}; |
| 70 | + |
| 71 | + Configurable<int> _sign{"sign", 1, "sign of a track"}; |
| 72 | + Configurable<float> _vertexZ{"VertexZ", 20.0, "abs vertexZ value limit"}; |
| 73 | + Configurable<float> _min_P{"min_P", 0.0, "lower mometum limit"}; |
| 74 | + Configurable<float> _max_P{"max_P", 100.0, "upper mometum limit"}; |
| 75 | + Configurable<float> _eta{"eta", 100.0, "abs eta value limit"}; |
| 76 | + |
| 77 | + Configurable<std::vector<float>> _dcaXY{"dcaXY", std::vector<float>{0.3f, 0.0f, 0.0f}, "abs dcaXY value limit; formula: [0] + [1]*pT^[2]"}; |
| 78 | + Configurable<std::vector<float>> _dcaZ{"dcaZ", std::vector<float>{0.3f, 0.0f, 0.0f}, "abs dcaZ value limit; formula: [0] + [1]*pT^[2]"}; |
| 79 | + Configurable<int16_t> _tpcNClsFound{"minTpcNClsFound", 0, "minimum allowed number of TPC clasters"}; |
| 80 | + Configurable<float> _tpcChi2NCl{"tpcChi2NCl", 100.0, "upper limit for chi2 value of a fit over TPC clasters"}; |
| 81 | + Configurable<float> _tpcCrossedRowsOverFindableCls{"tpcCrossedRowsOverFindableCls", 0, "lower limit of TPC CrossedRows/FindableCls value"}; |
| 82 | + Configurable<float> _tpcFractionSharedCls{"maxTpcFractionSharedCls", 0.4, "maximum fraction of TPC shared clasters"}; |
| 83 | + Configurable<int> _itsNCls{"minItsNCls", 0, "minimum allowed number of ITS clasters for a track"}; |
| 84 | + Configurable<float> _itsChi2NCl{"itsChi2NCl", 100.0, "upper limit for chi2 value of a fit over ITS clasters for a track"}; |
| 85 | + Configurable<int> _particlePDG{"particlePDG", 2212, "PDG code of a particle to perform PID for (only pion, kaon, proton and deurton are supported now)"}; |
| 86 | + Configurable<std::pair<float, float>> _tpcNSigma{"tpcNSigma", std::pair<float, float>{-100, 100}, "Nsigma range in TPC before the TOF is used"}; |
| 87 | + Configurable<std::pair<float, float>> _itsNSigma{"itsNSigma", std::pair<float, float>{-100, 100}, "Nsigma range in ITS to use along with TPC"}; |
| 88 | + Configurable<float> _PIDtrshld{"PIDtrshld", 10.0, "value of momentum from which the PID is done with TOF (before that only TPC is used)"}; |
| 89 | + Configurable<std::pair<float, float>> _tofNSigma{"tofNSigma", std::pair<float, float>{-100, 100}, "Nsigma range in TOF"}; |
| 90 | + Configurable<std::pair<float, float>> _tpcNSigmaResidual{"tpcNSigmaResidual", std::pair<float, float>{-10, 10}, "residual TPC Nsigma cut to use with the TOF"}; |
| 91 | + Configurable<int> _particlePDGtoReject{"particlePDGtoReject", 211, "PDG codes of perticles that will be rejected with TOF (only pion, kaon, proton and deurton are supported now)"}; |
| 92 | + Configurable<std::pair<float, float>> _rejectWithinNsigmaTOF{"rejectWithinNsigmaTOF", std::pair<float, float>{-10, 10}, "TOF rejection Nsigma range for particles specified with PDG to be rejected"}; |
| 93 | + |
| 94 | + std::shared_ptr<TH2> ITShisto; |
| 95 | + std::shared_ptr<TH2> TPChisto; |
| 96 | + std::shared_ptr<TH2> TOFhisto; |
| 97 | + std::shared_ptr<TH2> ITSvsTPChisto; |
| 98 | + std::shared_ptr<TH2> dcaxy_p; |
| 99 | + std::shared_ptr<TH2> dcaxy_pt; |
| 100 | + std::shared_ptr<TH2> dcaz_p; |
| 101 | + std::shared_ptr<TH2> dcaz_pt; |
| 102 | + |
| 103 | + void init(o2::framework::InitContext& context) |
| 104 | + { |
| 105 | + o2::aod::ITSResponse::setParameters(context); |
| 106 | + histos.add("vtz", "vtz", kTH1F, {{100, -20., 20., "vtxz"}}); |
| 107 | + histos.add("eta", "eta", kTH1F, {{200, -2.5, 2.5, "eta"}}); |
| 108 | + histos.add("phi", "phi", kTH1F, {{200, 0., 2. * M_PI, "phi"}}); |
| 109 | + histos.add("px", "px", kTH1F, {{100, 0., 5., "px"}}); |
| 110 | + histos.add("py", "py", kTH1F, {{100, 0., 5., "py"}}); |
| 111 | + histos.add("pz", "pz", kTH1F, {{100, 0., 5., "pz"}}); |
| 112 | + histos.add("p", "p", kTH1F, {{100, 0., 5., "p"}}); |
| 113 | + histos.add("pt", "pt", kTH1F, {{100, 0., 5., "pt"}}); |
| 114 | + histos.add("sign", "sign", kTH1F, {{3, -1.5, 1.5, "sign"}}); |
| 115 | + histos.add("TPCClusters", "TPCClusters", kTH1F, {{163, -0.5, 162.5, "NTPCClust"}}); |
| 116 | + histos.add("TPCCrossedRowsOverFindableCls", "TPCCrossedRowsOverFindableCls", kTH1F, {{100, 0.0, 10.0, "NcrossedRowsOverFindable"}}); |
| 117 | + histos.add("TPCFractionSharedCls", "TPCFractionSharedCls", kTH1F, {{100, 0.0, 1.0, "TPCsharedFraction"}}); |
| 118 | + histos.add("ITSClusters", "ITSClusters", kTH1F, {{10, -0.5, 9.5, "NITSClust"}}); |
| 119 | + histos.add("ITSchi2", "ITSchi2", kTH1F, {{100, 0.0, 40., "ITSchi2"}}); |
| 120 | + histos.add("TPCchi2", "TPCchi2", kTH1F, {{100, 0.0, 6., "TPCchi2"}}); |
| 121 | + |
| 122 | + dcaxy_p = histos.add<TH2>("dcaxy_p", "dcaxy_p", kTH2F, {{100, 0., 5.0, "p"}, {500, -0.5, 0.5, "dcaxy"}}); |
| 123 | + dcaxy_pt = histos.add<TH2>("dcaxy_pt", "dcaxy_pt", kTH2F, {{100, 0., 5.0, "pt"}, {500, -0.5, 0.5, "dcaxy"}}); |
| 124 | + dcaz_p = histos.add<TH2>("dcaz_p", "dcaz_p", kTH2F, {{100, 0., 5.0, "p"}, {500, -0.5, 0.5, "dcaz"}}); |
| 125 | + dcaz_pt = histos.add<TH2>("dcaz_pt", "dcaz_pt", kTH2F, {{100, 0., 5.0, "pt"}, {500, -0.5, 0.5, "dcaxy"}}); |
| 126 | + |
| 127 | + ITShisto = histos.add<TH2>(Form("nsigmaITS_PDG%i", _particlePDG.value), Form("nsigmaITS_PDG%i", _particlePDG.value), kTH2F, {{100, 0., 10.}, {1000, -50., 50.}}); |
| 128 | + TPChisto = histos.add<TH2>(Form("nsigmaTPC_PDG%i", _particlePDG.value), Form("nsigmaTPC_PDG%i", _particlePDG.value), kTH2F, {{100, 0., 10.}, {1000, -50., 50.}}); |
| 129 | + TOFhisto = histos.add<TH2>(Form("nsigmaTOF_PDG%i", _particlePDG.value), Form("nsigmaTOF_PDG%i", _particlePDG.value), kTH2F, {{100, 0., 10.}, {2000, -100., 100.}}); |
| 130 | + |
| 131 | + ITSvsTPChisto = histos.add<TH2>(Form("nsigmaITSvsTPC_PDG%i", _particlePDG.value), Form("nsigmaITSvsTPC_PDG%i", _particlePDG.value), kTH2F, {{1000, -50., 50.}, {1000, -50., 50.}}); |
| 132 | + } |
| 133 | + |
| 134 | + template <typename TrackType> |
| 135 | + inline float getITSNsigma(TrackType const& track, int const& PDG) |
| 136 | + { |
| 137 | + switch (PDG) { |
| 138 | + case 211: |
| 139 | + return track.itsNSigmaPi(); |
| 140 | + case 321: |
| 141 | + return track.itsNSigmaKa(); |
| 142 | + case 2212: |
| 143 | + return track.itsNSigmaPr(); |
| 144 | + case 1000010020: |
| 145 | + return track.itsNSigmaDe(); |
| 146 | + case 0: |
| 147 | + return -1000.0; |
| 148 | + default: |
| 149 | + LOG(fatal) << "Cannot interpret PDG for ITS selection: " << PDG; |
| 150 | + return -1000.0; |
| 151 | + } |
| 152 | + } |
| 153 | + template <typename TrackType> |
| 154 | + inline float getTPCNsigma(TrackType const& track, int const& PDG) |
| 155 | + { |
| 156 | + switch (PDG) { |
| 157 | + case 211: |
| 158 | + return track.tpcNSigmaPi(); |
| 159 | + case 321: |
| 160 | + return track.tpcNSigmaKa(); |
| 161 | + case 2212: |
| 162 | + return track.tpcNSigmaPr(); |
| 163 | + case 1000010020: |
| 164 | + return track.tpcNSigmaDe(); |
| 165 | + case 0: |
| 166 | + return -1000.0; |
| 167 | + default: |
| 168 | + LOG(fatal) << "Cannot interpret PDG for TPC selection: " << PDG; |
| 169 | + return -1000.0; |
| 170 | + } |
| 171 | + } |
| 172 | + template <typename TrackType> |
| 173 | + inline float getTOFNsigma(TrackType const& track, int const& PDG) |
| 174 | + { |
| 175 | + switch (PDG) { |
| 176 | + case 211: |
| 177 | + return track.tofNSigmaPi(); |
| 178 | + case 321: |
| 179 | + return track.tofNSigmaKa(); |
| 180 | + case 2212: |
| 181 | + return track.tofNSigmaPr(); |
| 182 | + case 1000010020: |
| 183 | + return track.tofNSigmaDe(); |
| 184 | + case 0: |
| 185 | + return -1000.0; |
| 186 | + default: |
| 187 | + LOG(fatal) << "Cannot interpret PDG for TOF selection: " << PDG; |
| 188 | + return -1000.0; |
| 189 | + } |
| 190 | + } |
| 191 | + |
| 192 | + bool isInRange(float value, std::pair<float, float> range) |
| 193 | + { |
| 194 | + return value > range.first && value < range.second; |
| 195 | + } |
| 196 | + |
| 197 | + void process(soa::Join<aod::Collisions, aod::EvSels> const& collisions, aod::SelectedTracks const& tracks) |
| 198 | + { |
| 199 | + auto tracksWithItsPid = soa::Attach<aod::SelectedTracks, aod::pidits::ITSNSigmaPi, aod::pidits::ITSNSigmaKa, aod::pidits::ITSNSigmaPr, aod::pidits::ITSNSigmaDe>(tracks); |
| 200 | + |
| 201 | + for (auto& collision : collisions) { |
| 202 | + if (!collision.sel8() || abs(collision.posZ()) > _vertexZ) |
| 203 | + continue; |
| 204 | + if (_requestGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsTriggerTVX)) |
| 205 | + continue; |
| 206 | + if (_removeSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) |
| 207 | + continue; |
| 208 | + if (_requestGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) |
| 209 | + continue; |
| 210 | + if (_requestVertexITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) |
| 211 | + continue; |
| 212 | + // if (_requestVertexTOForTRDmatched > collision.selection_bit(o2::aod::evsel::kisVertexTOForTRDmatched())) |
| 213 | + // continue; |
| 214 | + if (_requestNoCollInTimeRangeStandard && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) |
| 215 | + continue; |
| 216 | + // if (collision.multPerc() < _centCut.value.first || collision.multPerc() >= _centCut.value.second) |
| 217 | + // continue; |
| 218 | + // if (collision.hadronicRate() < _IRcut.value.first || collision.hadronicRate() >= _IRcut.value.second) |
| 219 | + // continue; |
| 220 | + if (collision.trackOccupancyInTimeRange() < _OccupancyCut.value.first || collision.trackOccupancyInTimeRange() >= _OccupancyCut.value.second) |
| 221 | + continue; |
| 222 | + |
| 223 | + histos.fill(HIST("vtz"), collision.posZ()); |
| 224 | + } |
| 225 | + |
| 226 | + for (auto& track : tracksWithItsPid) { |
| 227 | + if (track.sign() != _sign) |
| 228 | + continue; |
| 229 | + if (track.p() < _min_P || track.p() > _max_P || abs(track.eta()) > _eta) |
| 230 | + continue; |
| 231 | + if ((track.itsChi2NCl() >= _itsChi2NCl) || (track.itsChi2NCl() <= 0.f) || (track.tpcChi2NCl() <= 0.f) || (track.tpcChi2NCl() >= _tpcChi2NCl)) |
| 232 | + continue; |
| 233 | + if ((track.tpcFractionSharedCls()) > _tpcFractionSharedCls || (track.tpcNClsFound()) < _tpcNClsFound || (track.tpcCrossedRowsOverFindableCls()) < _tpcCrossedRowsOverFindableCls || (track.itsNCls()) < _itsNCls) |
| 234 | + continue; |
| 235 | + if (std::fabs(track.dcaXY()) > _dcaXY.value[0] + _dcaXY.value[1] * std::pow(track.pt(), _dcaXY.value[2]) || std::fabs(track.dcaZ()) > _dcaZ.value[0] + _dcaZ.value[1] * std::pow(track.pt(), _dcaZ.value[2])) |
| 236 | + continue; |
| 237 | + |
| 238 | + bool belowThreshold = track.p() < _PIDtrshld; |
| 239 | + float tpcSigma = getTPCNsigma(track, _particlePDG); |
| 240 | + float itsSigma = getITSNsigma(track, _particlePDG); |
| 241 | + float tofSigma = getTOFNsigma(track, _particlePDG); |
| 242 | + float tofRejection = getTOFNsigma(track, _particlePDGtoReject); |
| 243 | + |
| 244 | + bool passTPC = belowThreshold ? isInRange(tpcSigma, _tpcNSigma.value) : isInRange(tpcSigma, _tpcNSigmaResidual.value); |
| 245 | + |
| 246 | + bool passITS = belowThreshold ? isInRange(itsSigma, _itsNSigma.value) : true; |
| 247 | + |
| 248 | + bool passTOF = belowThreshold ? true : (isInRange(tofSigma, _tofNSigma.value) && !isInRange(tofRejection, _rejectWithinNsigmaTOF.value)); |
| 249 | + |
| 250 | + if (passTPC && passITS && passTOF) { |
| 251 | + histos.fill(HIST("eta"), track.eta()); |
| 252 | + histos.fill(HIST("phi"), track.phi()); |
| 253 | + histos.fill(HIST("px"), track.px()); |
| 254 | + histos.fill(HIST("py"), track.py()); |
| 255 | + histos.fill(HIST("pz"), track.pz()); |
| 256 | + histos.fill(HIST("p"), track.p()); |
| 257 | + histos.fill(HIST("pt"), track.pt()); |
| 258 | + histos.fill(HIST("sign"), track.sign()); |
| 259 | + histos.fill(HIST("dcaxy_p"), track.p(), track.dcaXY()); |
| 260 | + histos.fill(HIST("dcaxy_pt"), track.pt(), track.dcaXY()); |
| 261 | + histos.fill(HIST("dcaz_p"), track.p(), track.dcaZ()); |
| 262 | + histos.fill(HIST("dcaz_pt"), track.pt(), track.dcaZ()); |
| 263 | + histos.fill(HIST("TPCClusters"), track.tpcNClsFound()); |
| 264 | + histos.fill(HIST("TPCCrossedRowsOverFindableCls"), track.tpcCrossedRowsOverFindableCls()); |
| 265 | + histos.fill(HIST("TPCFractionSharedCls"), track.tpcFractionSharedCls()); |
| 266 | + histos.fill(HIST("ITSClusters"), track.itsNCls()); |
| 267 | + histos.fill(HIST("ITSchi2"), track.itsChi2NCl()); |
| 268 | + histos.fill(HIST("TPCchi2"), track.tpcChi2NCl()); |
| 269 | + |
| 270 | + switch (_particlePDG) { |
| 271 | + case 211: |
| 272 | + ITShisto->Fill(track.p(), track.itsNSigmaPi()); |
| 273 | + TPChisto->Fill(track.p(), track.tpcNSigmaPi()); |
| 274 | + TOFhisto->Fill(track.p(), track.tofNSigmaPi()); |
| 275 | + ITSvsTPChisto->Fill(track.itsNSigmaPi(), track.tpcNSigmaPi()); |
| 276 | + break; |
| 277 | + case 321: |
| 278 | + ITShisto->Fill(track.p(), track.itsNSigmaKa()); |
| 279 | + TPChisto->Fill(track.p(), track.tpcNSigmaKa()); |
| 280 | + TOFhisto->Fill(track.p(), track.tofNSigmaKa()); |
| 281 | + ITSvsTPChisto->Fill(track.itsNSigmaKa(), track.tpcNSigmaKa()); |
| 282 | + break; |
| 283 | + case 2212: |
| 284 | + ITShisto->Fill(track.p(), track.itsNSigmaPr()); |
| 285 | + TPChisto->Fill(track.p(), track.tpcNSigmaPr()); |
| 286 | + TOFhisto->Fill(track.p(), track.tofNSigmaPr()); |
| 287 | + ITSvsTPChisto->Fill(track.itsNSigmaPr(), track.tpcNSigmaPr()); |
| 288 | + break; |
| 289 | + case 1000010020: |
| 290 | + ITShisto->Fill(track.p(), track.itsNSigmaDe()); |
| 291 | + TPChisto->Fill(track.p(), track.tpcNSigmaDe()); |
| 292 | + TOFhisto->Fill(track.p(), track.tofNSigmaDe()); |
| 293 | + ITSvsTPChisto->Fill(track.itsNSigmaDe(), track.tpcNSigmaDe()); |
| 294 | + break; |
| 295 | + } |
| 296 | + } |
| 297 | + } |
| 298 | + } |
| 299 | +}; |
| 300 | + |
| 301 | +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) |
| 302 | +{ |
| 303 | + return WorkflowSpec{adaptAnalysisTask<PidOptimization>(cfgc)}; |
| 304 | +} |
0 commit comments