Skip to content

Commit 5c8a318

Browse files
njacazionjacazio
authored
[ALICE3] Implement proto PID with tracking layers (#11278)
1 parent 2d44451 commit 5c8a318

File tree

3 files changed

+270
-0
lines changed

3 files changed

+270
-0
lines changed

‎ALICE3/DataModel/OTFPIDTrk.h‎

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
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+
///
13+
/// \file OTFPIDTrk.h
14+
/// \author Berkin Ulukutlu TUM
15+
/// \author Henrik Fribert TUM
16+
/// \author Nicolò Jacazio Università del Piemonte Orientale
17+
/// \since May 22, 2025
18+
/// \brief Set of tables for the ALICE3 Trk PID information
19+
///
20+
21+
#ifndef ALICE3_DATAMODEL_OTFPIDTRK_H_
22+
#define ALICE3_DATAMODEL_OTFPIDTRK_H_
23+
24+
// O2 includes
25+
#include "Framework/AnalysisDataModel.h"
26+
27+
namespace o2::aod
28+
{
29+
namespace upgrade::trk
30+
{
31+
32+
DECLARE_SOA_COLUMN(TimeOverThresholdBarrel, timeOverThresholdBarrel, float); //! Time over threshold for the barrel layers
33+
DECLARE_SOA_COLUMN(ClusterSizeBarrel, clusterSizeBarrel, float); //! Cluster size for the barrel layers
34+
DECLARE_SOA_COLUMN(TimeOverThresholdForward, timeOverThresholdForward, float); //! Time over threshold for the Forward layers
35+
DECLARE_SOA_COLUMN(ClusterSizeForward, clusterSizeForward, float); //! Cluster size for the barrel layers
36+
37+
DECLARE_SOA_COLUMN(NSigmaTrkEl, nSigmaEl, float); //! NSigma electron from the tracker layers
38+
DECLARE_SOA_COLUMN(NSigmaTrkMu, nSigmaMu, float); //! NSigma muon from the tracker layers
39+
DECLARE_SOA_COLUMN(NSigmaTrkPi, nSigmaPi, float); //! NSigma pion from the tracker layers
40+
DECLARE_SOA_COLUMN(NSigmaTrkKa, nSigmaKa, float); //! NSigma kaon from the tracker layers
41+
DECLARE_SOA_COLUMN(NSigmaTrkPr, nSigmaPr, float); //! NSigma proton from the tracker layers
42+
43+
DECLARE_SOA_DYNAMIC_COLUMN(NSigmaTrk, nSigmaTrk, //! General function to get the nSigma for the tracker layers
44+
[](const float el,
45+
const float mu,
46+
const float pi,
47+
const float ka,
48+
const float pr,
49+
const int id) -> float {
50+
switch (std::abs(id)) {
51+
case 0:
52+
return el;
53+
case 1:
54+
return mu;
55+
case 2:
56+
return pi;
57+
case 3:
58+
return ka;
59+
case 4:
60+
return pr;
61+
default:
62+
LOG(fatal) << "Unrecognized PDG code for InnerTOF";
63+
return 999.f;
64+
}
65+
});
66+
67+
} // namespace upgrade::trk
68+
69+
DECLARE_SOA_TABLE(UpgradeTrkPidSignals, "AOD", "UPGRADETRKSIG",
70+
upgrade::trk::TimeOverThresholdBarrel,
71+
upgrade::trk::ClusterSizeBarrel);
72+
73+
DECLARE_SOA_TABLE(UpgradeTrkPids, "AOD", "UPGRADETRKPID",
74+
upgrade::trk::NSigmaTrkEl,
75+
upgrade::trk::NSigmaTrkMu,
76+
upgrade::trk::NSigmaTrkPi,
77+
upgrade::trk::NSigmaTrkKa,
78+
upgrade::trk::NSigmaTrkPr,
79+
upgrade::trk::NSigmaTrk<upgrade::trk::NSigmaTrkEl,
80+
upgrade::trk::NSigmaTrkMu,
81+
upgrade::trk::NSigmaTrkPi,
82+
upgrade::trk::NSigmaTrkKa,
83+
upgrade::trk::NSigmaTrkPr>);
84+
85+
using UpgradeTrkPidSignal = UpgradeTrkPidSignals::iterator;
86+
using UpgradeTrkPid = UpgradeTrkPids::iterator;
87+
88+
} // namespace o2::aod
89+
90+
#endif // ALICE3_DATAMODEL_OTFPIDTRK_H_

‎ALICE3/TableProducer/OTF/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -23,3 +23,8 @@ o2physics_add_dpl_workflow(onthefly-richpid
2323
SOURCES onTheFlyRichPid.cxx
2424
PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core
2525
COMPONENT_NAME Analysis)
26+
27+
o2physics_add_dpl_workflow(onthefly-trkpid
28+
SOURCES onTheFlyTrackerPid.cxx
29+
PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core
30+
COMPONENT_NAME Analysis)

‎ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx‎

Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
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+
/// \file onTheFlyTrackerPid.cxx
13+
///
14+
/// \brief This task produces the PID information that can be obtained from the tracker layers (i.e. cluster size and ToT).
15+
/// It currently contemplates 5 particle types: electrons, muons, pions, kaons and protons.
16+
///
17+
/// \author Berkin Ulukutlu TUM
18+
/// \author Henrik Fribert TUM
19+
/// \author Nicolò Jacazio Università del Piemonte Orientale
20+
/// \since May 22, 2025
21+
///
22+
23+
#include <utility>
24+
#include <map>
25+
#include <string>
26+
#include <algorithm>
27+
#include <vector>
28+
29+
#include "Framework/AnalysisDataModel.h"
30+
#include "Framework/AnalysisTask.h"
31+
#include "Framework/runDataProcessing.h"
32+
#include "Framework/RunningWorkflowInfo.h"
33+
#include "Framework/HistogramRegistry.h"
34+
#include "Framework/O2DatabasePDGPlugin.h"
35+
#include "Framework/ASoAHelpers.h"
36+
#include "Common/DataModel/TrackSelectionTables.h"
37+
#include "Common/Core/trackUtilities.h"
38+
#include "ALICE3/Core/TrackUtilities.h"
39+
#include "ReconstructionDataFormats/DCA.h"
40+
#include "DetectorsBase/Propagator.h"
41+
#include "DetectorsBase/GeometryManager.h"
42+
#include "CommonUtils/NameConf.h"
43+
#include "CCDB/CcdbApi.h"
44+
#include "CCDB/BasicCCDBManager.h"
45+
#include "DataFormatsParameters/GRPMagField.h"
46+
#include "DataFormatsCalibration/MeanVertexObject.h"
47+
#include "CommonConstants/GeomConstants.h"
48+
#include "CommonConstants/PhysicsConstants.h"
49+
#include "TRandom3.h"
50+
#include "TF1.h"
51+
#include "TH2F.h"
52+
#include "TVector3.h"
53+
#include "TString.h"
54+
#include "ALICE3/DataModel/OTFRICH.h"
55+
#include "DetectorsVertexing/HelixHelper.h"
56+
#include "TableHelper.h"
57+
#include "ALICE3/Core/DelphesO2TrackSmearer.h"
58+
#include "ALICE3/DataModel/OTFPIDTrk.h"
59+
60+
using namespace o2;
61+
using namespace o2::framework;
62+
63+
struct OnTheFlyTrackerPid {
64+
Produces<aod::UpgradeTrkPidSignals> tableUpgradeTrkPidSignals;
65+
Produces<aod::UpgradeTrkPids> tableUpgradeTrkPids;
66+
67+
// necessary for particle charges
68+
Service<o2::framework::O2DatabasePDG> pdg;
69+
70+
static constexpr int kMaxBarrelLayers = 8;
71+
static constexpr int kMaxForwardLayers = 9;
72+
73+
struct : ConfigurableGroup {
74+
Configurable<std::string> efficiencyFormula{"efficiencyFormula", "1.0/(1.0+exp(-(x-0.01)/0.2))", "ROOT TF1 formula for efficiency"};
75+
Configurable<std::string> landauFormula{"landauFormula", "TMath::Landau(x, 1, 1, true)", "ROOT TF1 formula for Landau distribution (e.g. ToT response)"};
76+
Configurable<int> averageMethod{"averageMethod", 0, "Method to average the ToT and cluster size. 0: truncated mean"};
77+
} simConfig;
78+
79+
TF1* mEfficiency = nullptr;
80+
static constexpr int kEtaBins = 50;
81+
static constexpr float kEtaMin = -2.5;
82+
static constexpr float kEtaMax = 2.5;
83+
static constexpr int kPtBins = 200;
84+
static constexpr float kPtMin = 0.0;
85+
static constexpr float kPtMax = 20.0;
86+
87+
std::array<std::array<TF1*, kPtBins>, kEtaBins> mElossPi;
88+
89+
void init(o2::framework::InitContext&)
90+
{
91+
92+
for (int i = 0; i < kEtaBins; i++) {
93+
for (int j = 0; j < kPtBins; j++) {
94+
mElossPi[i][j] = new TF1(Form("mElossPi_%d_%d", i, j), simConfig.landauFormula.value.c_str(), 0, 20);
95+
}
96+
}
97+
mEfficiency = new TF1("mEfficiency", simConfig.efficiencyFormula.value.c_str(), 0, 20);
98+
}
99+
100+
void process(soa::Join<aod::Collisions, aod::McCollisionLabels>::iterator const&,
101+
soa::Join<aod::Tracks, aod::TracksCov, aod::McTrackLabels> const& tracks,
102+
aod::McParticles const&,
103+
aod::McCollisions const&)
104+
{
105+
std::array<float, kMaxBarrelLayers> timeOverThresholdBarrel;
106+
std::array<float, kMaxBarrelLayers> clusterSizeBarrel;
107+
// std::array<float, kMaxForwardLayers> timeOverThresholdForward;
108+
// std::array<float, kMaxForwardLayers> clusterSizeForward;
109+
110+
auto noSignalTrack = [&]() {
111+
tableUpgradeTrkPidSignals(0.f, 0.f); // no PID information
112+
tableUpgradeTrkPids(0.f, 0.f, 0.f, 0.f, 0.f); // no PID information
113+
};
114+
115+
for (const auto& track : tracks) {
116+
if (!track.has_mcParticle()) {
117+
noSignalTrack();
118+
continue;
119+
}
120+
const auto& mcParticle = track.mcParticle();
121+
const auto& pdgInfo = pdg->GetParticle(mcParticle.pdgCode());
122+
if (!pdgInfo) {
123+
LOG(warning) << "PDG code " << mcParticle.pdgCode() << " not found in the database";
124+
noSignalTrack();
125+
continue;
126+
}
127+
const float pt = mcParticle.pt();
128+
const float eta = mcParticle.eta();
129+
130+
const int binnedPt = static_cast<int>((pt - kPtMin) / kPtBins);
131+
const int binnedEta = static_cast<int>((eta - kEtaMin) / kEtaBins);
132+
if (binnedPt < 0 || binnedPt >= kPtBins || binnedEta < 0 || binnedEta >= kEtaBins) {
133+
noSignalTrack();
134+
continue;
135+
}
136+
for (int i = 0; i < kMaxBarrelLayers; i++) {
137+
timeOverThresholdBarrel[i] = -1;
138+
clusterSizeBarrel[i] = -1;
139+
140+
// Check if layer is efficient
141+
if (mEfficiency->Eval(pt) > gRandom->Uniform(0, 1)) {
142+
timeOverThresholdBarrel[i] = mElossPi[binnedEta][binnedPt]->GetRandom(); // Simulate ToT
143+
clusterSizeBarrel[i] = mElossPi[binnedEta][binnedPt]->GetRandom(); // Simulate cluster size
144+
}
145+
}
146+
147+
// Now we do the average
148+
switch (simConfig.averageMethod) {
149+
case 0: { // truncated mean
150+
float meanToT = 0;
151+
float meanClusterSize = 0;
152+
// Order them by ToT
153+
std::sort(timeOverThresholdBarrel.begin(), timeOverThresholdBarrel.end());
154+
std::sort(clusterSizeBarrel.begin(), clusterSizeBarrel.end());
155+
static constexpr int kTruncatedMean = 5;
156+
// Take the mean of the first 5 values
157+
for (int i = 0; i < kTruncatedMean; i++) {
158+
meanToT += timeOverThresholdBarrel[i];
159+
meanClusterSize += clusterSizeBarrel[i];
160+
}
161+
meanToT /= kTruncatedMean;
162+
meanClusterSize /= kTruncatedMean;
163+
// Fill the table
164+
tableUpgradeTrkPidSignals(meanToT, meanClusterSize);
165+
} break;
166+
167+
default:
168+
LOG(fatal) << "Unknown average method " << simConfig.averageMethod;
169+
break;
170+
}
171+
}
172+
}
173+
};
174+
175+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<OnTheFlyTrackerPid>(cfgc)}; }

0 commit comments

Comments
 (0)