Skip to content

Commit 346c9ef

Browse files
nburmasonburmaso
authored
[PWGUD] Introduce compact MCH and MCH-MID efficiency table producer (#18569)
1 parent 8eda4f5 commit 346c9ef

2 files changed

Lines changed: 264 additions & 0 deletions

File tree

‎PWGUD/TableProducer/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -100,3 +100,8 @@ o2physics_add_dpl_workflow(upc-cand-producer-barrel
100100
SOURCES upcCandProducerBarrel.cxx
101101
PUBLIC_LINK_LIBRARIES O2::Framework O2::CCDB O2::DataFormatsParameters O2Physics::AnalysisCore
102102
COMPONENT_NAME Analysis)
103+
104+
o2physics_add_dpl_workflow(muon-efficiency-producer
105+
SOURCES muonEfficiencyProducer.cxx
106+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::GlobalTracking
107+
COMPONENT_NAME Analysis)

‎PWGUD/TableProducer/muonEfficiencyProducer.cxx‎

Lines changed: 259 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,259 @@
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 muonEfficiencyProducer.cxx
13+
/// \brief Compact forward-muon skim for offline MCH and MCH-MID efficiency studies
14+
/// \author Nazar Burmasov (JINR)
15+
16+
#include <CCDB/BasicCCDBManager.h>
17+
#include <DataFormatsParameters/AggregatedRunInfo.h>
18+
#include <DataFormatsParameters/GRPMagField.h>
19+
#include <DetectorsBase/GeometryManager.h>
20+
#include <DetectorsBase/Propagator.h>
21+
#include <Framework/ASoA.h>
22+
#include <Framework/AnalysisDataModel.h>
23+
#include <Framework/AnalysisHelpers.h>
24+
#include <Framework/AnalysisTask.h>
25+
#include <Framework/Configurable.h>
26+
#include <Framework/DataTypes.h>
27+
#include <Framework/InitContext.h>
28+
#include <Framework/Logger.h>
29+
#include <Framework/runDataProcessing.h>
30+
#include <GlobalTracking/MatchGlobalFwd.h>
31+
#include <MCHTracking/TrackExtrap.h>
32+
#include <ReconstructionDataFormats/GlobalFwdTrack.h>
33+
34+
#include <Math/MatrixRepresentationsStatic.h>
35+
#include <Math/SMatrix.h>
36+
#include <TGeoManager.h>
37+
38+
#include <array>
39+
#include <cstdint>
40+
#include <vector>
41+
42+
using namespace o2;
43+
using namespace o2::framework;
44+
45+
namespace o2::aod::mcheffskim
46+
{
47+
DECLARE_SOA_COLUMN(RunNumber, runNumber, int32_t);
48+
DECLARE_SOA_COLUMN(ReferenceBC, referenceBC, int64_t);
49+
DECLARE_SOA_COLUMN(FT0AAmplitude, ft0aAmplitude, float);
50+
DECLARE_SOA_COLUMN(FT0CAmplitude, ft0cAmplitude, float);
51+
DECLARE_SOA_COLUMN(FV0AAmplitude, fv0aAmplitude, float);
52+
DECLARE_SOA_COLUMN(CollisionId, collisionId, int32_t); // -1 for ambiguous tracks
53+
DECLARE_SOA_COLUMN(Pt, pt, float);
54+
DECLARE_SOA_COLUMN(Eta, eta, float);
55+
DECLARE_SOA_COLUMN(Phi, phi, float);
56+
DECLARE_SOA_COLUMN(Sign, sign, int8_t);
57+
DECLARE_SOA_COLUMN(RAbs, rAbs, float);
58+
DECLARE_SOA_COLUMN(PDca, pDca, float);
59+
DECLARE_SOA_COLUMN(Chi2Ndf, chi2Ndf, float); // -1 if propagation failed
60+
DECLARE_SOA_COLUMN(NClusters, nClusters, int8_t);
61+
DECLARE_SOA_COLUMN(MchBitMap, mchBitMap, uint16_t);
62+
DECLARE_SOA_COLUMN(MidBitMap, midBitMap, uint8_t);
63+
DECLARE_SOA_COLUMN(MidBoards, midBoards, uint32_t);
64+
DECLARE_SOA_COLUMN(Chi2MatchMCHMID, chi2MatchMCHMID, float); // -1 for MCH-only tracks
65+
DECLARE_SOA_COLUMN(TrackTime, trackTime, float); // ns relative to referenceBC
66+
DECLARE_SOA_COLUMN(TrackTimeRes, trackTimeRes, float);
67+
} // namespace o2::aod::mcheffskim
68+
69+
namespace o2::aod
70+
{
71+
DECLARE_SOA_TABLE(MchEffBCs, "AOD", "MCHEFFBC",
72+
o2::soa::Index<>,
73+
mcheffskim::RunNumber, mcheffskim::ReferenceBC,
74+
mcheffskim::FT0AAmplitude, mcheffskim::FT0CAmplitude,
75+
mcheffskim::FV0AAmplitude);
76+
} // namespace o2::aod
77+
78+
namespace o2::aod::mcheffskim
79+
{
80+
DECLARE_SOA_INDEX_COLUMN(MchEffBC, mchEffBC); // links a muon to the compact BC table
81+
} // namespace o2::aod::mcheffskim
82+
83+
namespace o2::aod
84+
{
85+
DECLARE_SOA_TABLE(MchEffSkims, "AOD", "MCHEFFSKIM",
86+
o2::soa::Index<>, mcheffskim::MchEffBCId,
87+
mcheffskim::CollisionId,
88+
mcheffskim::Pt, mcheffskim::Eta, mcheffskim::Phi, mcheffskim::Sign,
89+
mcheffskim::RAbs, mcheffskim::PDca, mcheffskim::Chi2Ndf,
90+
mcheffskim::NClusters, mcheffskim::MchBitMap,
91+
mcheffskim::MidBitMap, mcheffskim::MidBoards,
92+
mcheffskim::Chi2MatchMCHMID,
93+
mcheffskim::TrackTime, mcheffskim::TrackTimeRes);
94+
} // namespace o2::aod
95+
96+
struct MuonEfficiencyProducer {
97+
static constexpr float MinEta = -4.f;
98+
static constexpr float MaxEta = -2.5f;
99+
static constexpr float NoAmplitude = -1.f;
100+
static constexpr float MaxFITTime = 30.f;
101+
102+
using Tracks = soa::Join<aod::FwdTracks, aod::FwdTracksCov>;
103+
using Matrix55 = ROOT::Math::SMatrix<double, 5, 5, ROOT::Math::MatRepSym<double, 5>>;
104+
using Vector5 = ROOT::Math::SVector<double, 5>;
105+
106+
Produces<aod::MchEffBCs> outputBCs;
107+
Produces<aod::MchEffSkims> output;
108+
Service<ccdb::BasicCCDBManager> ccdb{};
109+
globaltracking::MatchGlobalFwd matcher;
110+
int cachedFieldRun = -1;
111+
112+
Configurable<float> minPt{"minPt", 0.4f, "Minimum track pT (Gev/c)"};
113+
Configurable<float> minRAbs{"minRAbs", 17.6f, "Minimum absorber-end radius (cm)"};
114+
Configurable<float> maxRAbs{"maxRAbs", 89.5f, "Maximum absorber-end radius (cm)"};
115+
Configurable<float> maxPDca{"maxPDca", 400.f, "Maximum stored pDCA (GeV cm/c)"};
116+
Configurable<float> maxChi2{"maxChi2", 100.f, "Maximum MCH chi2/ndf"};
117+
118+
void init(InitContext&)
119+
{
120+
ccdb->setLocalObjectValidityChecking();
121+
}
122+
123+
void setField(int run)
124+
{
125+
if (cachedFieldRun == run) {
126+
return;
127+
}
128+
const auto runInfo = parameters::AggregatedRunInfo::buildAggregatedRunInfo(ccdb->instance(), run);
129+
const auto ts = (runInfo.sor + runInfo.eor) / 2;
130+
auto* grp = ccdb->getForTimeStamp<parameters::GRPMagField>("GLO/Config/GRPMagField", ts);
131+
if (!grp) {
132+
LOGF(fatal, "GRPMagField not found for run %d", run);
133+
}
134+
base::Propagator::initFieldFromGRP(grp);
135+
if (!base::GeometryManager::isGeometryLoaded()) {
136+
ccdb->get<TGeoManager>("GLO/Config/GeometryAligned");
137+
}
138+
mch::TrackExtrap::setField();
139+
cachedFieldRun = run;
140+
}
141+
142+
struct Kinematics {
143+
float pt, eta, phi;
144+
};
145+
146+
template <typename Track>
147+
bool propagateMCH(const Track& t, const std::array<float, 3>& vtx, const std::array<float, 2>& vtxCov, Kinematics& k)
148+
{
149+
const Vector5 params(t.x(), t.y(), t.phi(), t.tgl(), t.signed1Pt());
150+
const std::vector<double> covElements{t.cXX(), t.cXY(), t.cYY(), t.cPhiX(), t.cPhiY(), t.cPhiPhi(),
151+
t.cTglX(), t.cTglY(), t.cTglPhi(), t.cTglTgl(),
152+
t.c1PtX(), t.c1PtY(), t.c1PtPhi(), t.c1PtTgl(), t.c1Pt21Pt2()};
153+
const Matrix55 cov(covElements.begin(), covElements.end());
154+
dataformats::GlobalFwdTrack source;
155+
source.setParameters(params);
156+
source.setZ(t.z());
157+
source.setCovariances(cov);
158+
auto mchTrack = matcher.FwdtoMCH(source);
159+
if (!mch::TrackExtrap::extrapToVertex(mchTrack, vtx[0], vtx[1], vtx[2], vtxCov[0], vtxCov[1])) {
160+
return false;
161+
}
162+
const auto propagated = matcher.MCHtoFwd(mchTrack);
163+
k = {.pt = static_cast<float>(propagated.getPt()), .eta = static_cast<float>(propagated.getEta()), .phi = static_cast<float>(propagated.getPhi())};
164+
return true;
165+
}
166+
167+
void process(aod::BCs const& bcs, aod::FT0s const& ft0s, aod::FV0As const& fv0as,
168+
Tracks const& tracks, aod::Collisions const& collisions,
169+
aod::AmbiguousFwdTracks const& ambiguous)
170+
{
171+
if (bcs.size() == 0) {
172+
return;
173+
}
174+
const int run = bcs.begin().runNumber();
175+
176+
std::vector<float> ft0aByBC(bcs.size(), NoAmplitude);
177+
std::vector<float> ft0cByBC(bcs.size(), NoAmplitude);
178+
std::vector<float> fv0aByBC(bcs.size(), NoAmplitude);
179+
std::vector<int32_t> skimBCIds(bcs.size(), -1);
180+
for (const auto& ft0 : ft0s) {
181+
if (ft0.timeA() < MaxFITTime) {
182+
ft0aByBC[ft0.bcId()] = ft0.sumAmpA();
183+
}
184+
if (ft0.timeC() < MaxFITTime) {
185+
ft0cByBC[ft0.bcId()] = ft0.sumAmpC();
186+
}
187+
}
188+
for (const auto& fv0 : fv0as) {
189+
if (fv0.time() >= MaxFITTime) {
190+
continue;
191+
}
192+
float amplitude = 0.f;
193+
for (const float a : fv0.amplitude()) {
194+
if (a > 0.f) {
195+
amplitude += a;
196+
}
197+
}
198+
fv0aByBC[fv0.bcId()] = amplitude;
199+
}
200+
201+
std::vector<int32_t> ambiguousBCIds(tracks.size(), -1);
202+
for (const auto& entry : ambiguous) {
203+
const auto& ids = entry.bcIds();
204+
if (ids.empty()) {
205+
continue;
206+
}
207+
ambiguousBCIds[entry.fwdtrackId()] = *ids.begin();
208+
}
209+
210+
for (const auto& t : tracks) {
211+
const auto type = t.trackType();
212+
if (type != aod::fwdtrack::MuonStandaloneTrack && type != aod::fwdtrack::MCHStandaloneTrack) {
213+
continue;
214+
}
215+
const float rAbs = t.rAtAbsorberEnd();
216+
if (rAbs < minRAbs || rAbs > maxRAbs || t.pDca() > maxPDca || t.chi2() > maxChi2) {
217+
continue;
218+
}
219+
Kinematics k{.pt = t.pt(), .eta = t.eta(), .phi = t.phi()};
220+
std::array<float, 3> vtx{0.f, 0.f, 0.f};
221+
std::array<float, 2> vtxCov{0.f, 0.f};
222+
int32_t referenceBCId = ambiguousBCIds[t.globalIndex()];
223+
if (t.has_collision()) {
224+
const auto col = collisions.iteratorAt(t.collisionId());
225+
vtx = {col.posX(), col.posY(), col.posZ()};
226+
vtxCov = {col.covXX(), col.covYY()};
227+
referenceBCId = col.bcId();
228+
}
229+
setField(run);
230+
if (!propagateMCH(t, vtx, vtxCov, k)) {
231+
continue;
232+
}
233+
if (k.eta < MinEta || k.eta > MaxEta) {
234+
continue;
235+
}
236+
if (k.pt < minPt) {
237+
continue;
238+
}
239+
int32_t skimBCId = -1;
240+
if (referenceBCId >= 0) {
241+
skimBCId = skimBCIds[referenceBCId];
242+
if (skimBCId < 0) {
243+
outputBCs(run, bcs.iteratorAt(referenceBCId).globalBC(),
244+
ft0aByBC[referenceBCId], ft0cByBC[referenceBCId], fv0aByBC[referenceBCId]);
245+
skimBCId = outputBCs.lastIndex();
246+
skimBCIds[referenceBCId] = skimBCId;
247+
}
248+
}
249+
output(skimBCId, t.collisionId(), k.pt, k.eta, k.phi, static_cast<int8_t>(t.signed1Pt() > 0 ? 1 : -1),
250+
rAbs, t.pDca(), t.chi2(), t.nClusters(), t.mchBitMap(), t.midBitMap(), t.midBoards(),
251+
type == aod::fwdtrack::MuonStandaloneTrack ? t.chi2MatchMCHMID() : -1.f, t.trackTime(), t.trackTimeRes());
252+
}
253+
}
254+
};
255+
256+
WorkflowSpec defineDataProcessing(ConfigContext const& context)
257+
{
258+
return WorkflowSpec{adaptAnalysisTask<MuonEfficiencyProducer>(context)};
259+
}

0 commit comments

Comments
 (0)