Project
Loading...
Searching...
No Matches
TrackMCStudy.cxx
Go to the documentation of this file.
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#include <vector>
13#include <TStopwatch.h>
53#include "GPUO2InterfaceRefit.h"
54#include "GPUParam.h"
55#include "GPUParam.inc"
56#include "MathUtils/fit.h"
57#include <TRandom.h>
58#include <map>
59#include <unordered_map>
60#include <array>
61#include <utility>
62#include <gsl/span>
63
64// workflow to study relation of reco tracks to MCTruth
65// o2-trackmc-study-workflow --device-verbosity 3 -b --run
66
67namespace o2::trackstudy
68{
69
70using namespace o2::framework;
73
77using VTIndexV = std::pair<int, o2::dataformats::VtxTrackIndex>;
80
82
83class TrackMCStudy : public Task
84{
85 public:
86 TrackMCStudy(std::shared_ptr<DataRequest> dr, std::shared_ptr<o2::base::GRPGeomRequest> gr, GTrackID::mask_t src, const o2::tpc::CorrectionMapsLoaderGloOpts& sclOpts, bool checkSV)
87 : mDataRequest(dr), mGGCCDBRequest(gr), mTracksSrc(src), mCheckSV(checkSV)
88 {
89 mTPCCorrMapsLoader.setLumiScaleType(sclOpts.lumiType);
90 mTPCCorrMapsLoader.setLumiScaleMode(sclOpts.lumiMode);
91 mTPCCorrMapsLoader.setCheckCTPIDCConsistency(sclOpts.checkCTPIDCconsistency);
92 }
93 ~TrackMCStudy() final = default;
94 void init(InitContext& ic) final;
95 void run(ProcessingContext& pc) final;
96 void endOfStream(EndOfStreamContext& ec) final;
97 void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final;
98 void process(const o2::globaltracking::RecoContainer& recoData);
99
100 private:
101 void processTPCTrackRefs();
102 void loadTPCOccMap(const o2::globaltracking::RecoContainer& recoData);
103 void fillMCClusterInfo(const o2::globaltracking::RecoContainer& recoData);
104 void prepareITSData(const o2::globaltracking::RecoContainer& recoData);
105 bool processMCParticle(int src, int ev, int trid);
106 bool addMCParticle(const MCTrack& mctr, const o2::MCCompLabel& lb, TParticlePDG* pPDG = nullptr);
107 bool acceptMCCharged(const MCTrack& tr, const o2::MCCompLabel& lb, int followDec = -1);
108 bool propagateToRefX(o2::track::TrackParCov& trcTPC, o2::track::TrackParCov& trcITS);
109 bool refitV0(int i, o2::dataformats::V0& v0, const o2::globaltracking::RecoContainer& recoData);
110 void updateTimeDependentParams(ProcessingContext& pc);
111 float getDCAYCut(float pt) const;
112
113 const std::vector<o2::MCTrack>* mCurrMCTracks = nullptr;
114 TVector3 mCurrMCVertex;
115 o2::tpc::VDriftHelper mTPCVDriftHelper{};
116 o2::tpc::CorrectionMapsLoader mTPCCorrMapsLoader{};
117 std::shared_ptr<DataRequest> mDataRequest;
118 std::shared_ptr<o2::base::GRPGeomRequest> mGGCCDBRequest;
119 std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOut;
120 std::vector<float> mTBinClOcc;
121 std::vector<float> mTBinClOccHist; //< original occupancy
122 std::vector<long> mIntBC;
123 std::vector<float> mTPCOcc;
124 std::vector<int> mITSOcc; //< N ITS clusters in the ROF containing collision
125 bool mCheckSV = false; //< check SV binding (apart from prongs availability)
126 bool mRecProcStage = false; //< flag that the MC particle was added only at the stage of reco tracks processing
127 int mNTPCOccBinLength = 0;
128 float mNTPCOccBinLengthInv = -1.f;
129 int mVerbose = 0;
130 float mITSTimeBiasMUS = 0.f;
131 float mITSROFrameLengthMUS = 0.f;
132 float mTPCTBinMUS = 0.;
133
134 int mNCheckDecays = 0;
135
136 GTrackID::mask_t mTracksSrc{};
137 o2::steer::MCKinematicsReader mcReader; // reader of MC information
138 std::vector<int> mITSROF;
139 std::vector<TBracket> mITSROFBracket;
140 std::vector<o2::MCCompLabel> mDecProdLblPool; // labels of decay products to watch, added to MC map
141 std::vector<MCVertex> mMCVtVec{};
142
143 struct DecayRef {
144 o2::MCCompLabel mother{};
145 o2::track::TrackPar parent{};
146 int pdg = 0;
147 int daughterFirst = -1;
148 int daughterLast = -1;
149 int foundSVID = -1;
150 };
151 std::vector<std::vector<DecayRef>> mDecaysMaps; // for every parent particle to watch, store its label and entries of 1st/last decay product labels in mDecProdLblPool
152 std::unordered_map<o2::MCCompLabel, TrackFamily> mSelMCTracks;
153 std::unordered_map<o2::MCCompLabel, std::pair<int, int>> mSelTRefIdx;
154 std::vector<o2::track::TrackPar> mSelTRefs;
156 static constexpr float MaxSnp = 0.9; // max snp of ITS or TPC track at xRef to be matched
157};
158
160{
162 mcReader.initFromDigitContext("collisioncontext.root");
163
164 mDBGOut = std::make_unique<o2::utils::TreeStreamRedirector>("trackMCStudy.root", "recreate");
165 mVerbose = ic.options().get<int>("device-verbosity");
166
168 for (int id = 0; id < sizeof(params.decayPDG) / sizeof(int); id++) {
169 if (params.decayPDG[id] < 0) {
170 break;
171 }
172 mNCheckDecays++;
173 }
174 mDecaysMaps.resize(mNCheckDecays);
175 mTPCCorrMapsLoader.init(ic);
176}
177
179{
181 for (int i = 0; i < mNCheckDecays; i++) {
182 mDecaysMaps[i].clear();
183 }
184 mDecProdLblPool.clear();
185 mMCVtVec.clear();
186 mCurrMCTracks = nullptr;
187
188 recoData.collectData(pc, *mDataRequest.get()); // select tracks of needed type, with minimal cuts, the real selected will be done in the vertexer
189 updateTimeDependentParams(pc); // Make sure this is called after recoData.collectData, which may load some conditions
190 mRecProcStage = false;
191 process(recoData);
192}
193
194void TrackMCStudy::updateTimeDependentParams(ProcessingContext& pc)
195{
197 mTPCVDriftHelper.extractCCDBInputs(pc);
198 mTPCCorrMapsLoader.extractCCDBInputs(pc);
199 bool updateMaps = false;
200 if (mTPCCorrMapsLoader.isUpdated()) {
201 mTPCCorrMapsLoader.acknowledgeUpdate();
202 updateMaps = true;
203 }
204 if (mTPCVDriftHelper.isUpdated()) {
205 LOGP(info, "Updating TPC fast transform map with new VDrift factor of {} wrt reference {} and DriftTimeOffset correction {} wrt {} from source {}",
206 mTPCVDriftHelper.getVDriftObject().corrFact, mTPCVDriftHelper.getVDriftObject().refVDrift,
207 mTPCVDriftHelper.getVDriftObject().timeOffsetCorr, mTPCVDriftHelper.getVDriftObject().refTimeOffset,
208 mTPCVDriftHelper.getSourceName());
209 mTPCVDriftHelper.acknowledgeUpdate();
210 updateMaps = true;
211 }
212 if (updateMaps) {
213 mTPCCorrMapsLoader.updateVDrift(mTPCVDriftHelper.getVDriftObject().corrFact, mTPCVDriftHelper.getVDriftObject().refVDrift, mTPCVDriftHelper.getVDriftObject().getTimeOffset());
214 }
215 static bool initOnceDone = false;
216 if (!initOnceDone) { // this params need to be queried only once
217 initOnceDone = true;
219 mITSROFrameLengthMUS = o2::base::GRPGeomHelper::instance().getGRPECS()->isDetContinuousReadOut(o2::detectors::DetID::ITS) ? alpParamsITS.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingMUS : alpParamsITS.roFrameLengthTrig * 1.e-3;
220 LOGP(info, "VertexTrackMatcher ITSROFrameLengthMUS:{}", mITSROFrameLengthMUS);
221
223 mTPCTBinMUS = elParam.ZbinWidth;
224
225 if (mCheckSV) {
226 const auto& svparam = o2::vertexing::SVertexerParams::Instance();
227 mFitterV0.setBz(o2::base::Propagator::Instance()->getNominalBz());
228 mFitterV0.setUseAbsDCA(svparam.useAbsDCA);
229 mFitterV0.setPropagateToPCA(false);
230 mFitterV0.setMaxR(svparam.maxRIni);
231 mFitterV0.setMinParamChange(svparam.minParamChange);
232 mFitterV0.setMinRelChi2Change(svparam.minRelChi2Change);
233 mFitterV0.setMaxDZIni(svparam.maxDZIni);
234 mFitterV0.setMaxDXYIni(svparam.maxDXYIni);
235 mFitterV0.setMaxChi2(svparam.maxChi2);
236 mFitterV0.setMatCorrType(o2::base::Propagator::MatCorrType(svparam.matCorr));
237 mFitterV0.setUsePropagator(svparam.usePropagator);
238 mFitterV0.setRefitWithMatCorr(svparam.refitWithMatCorr);
239 mFitterV0.setMaxStep(svparam.maxStep);
240 mFitterV0.setMaxSnp(svparam.maxSnp);
241 mFitterV0.setMinXSeed(svparam.minXSeed);
242 }
243 }
244}
245
247{
248 constexpr float SQRT12Inv = 0.288675f;
250 auto pvvec = recoData.getPrimaryVertices();
251 auto pvvecLbl = recoData.getPrimaryVertexMCLabels();
252 auto trackIndex = recoData.getPrimaryVertexMatchedTracks(); // Global ID's for associated tracks
253 auto vtxRefs = recoData.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs
254 auto prop = o2::base::Propagator::Instance();
255 int nv = vtxRefs.size();
256 float vdriftTB = mTPCVDriftHelper.getVDriftObject().getVDrift() * o2::tpc::ParameterElectronics::Instance().ZbinWidth; // VDrift expressed in cm/TimeBin
257 float itsBias = 0.5 * mITSROFrameLengthMUS + o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameBiasInBC * o2::constants::lhc::LHCBunchSpacingMUS; // ITS time is supplied in \mus as beginning of ROF
258
259 prepareITSData(recoData);
260 loadTPCOccMap(recoData);
261 auto getITSPatt = [&](GTrackID gid, uint8_t& ncl) {
262 int8_t patt = 0;
263 if (gid.getSource() == VTIndex::ITSAB) {
264 const auto& itsTrf = recoData.getITSABRefs()[gid];
265 ncl = itsTrf.getNClusters();
266 for (int il = 0; il < 7; il++) {
267 if (itsTrf.hasHitOnLayer(il)) {
268 patt |= 0x1 << il;
269 }
270 }
271 patt |= 0x1 << 7;
272 } else {
273 const auto& itsTr = recoData.getITSTrack(gid);
274 for (int il = 0; il < 7; il++) {
275 if (itsTr.hasHitOnLayer(il)) {
276 patt |= 0x1 << il;
277 ncl++;
278 }
279 }
280 }
281 return patt;
282 };
283
284 auto fillTPCClusterInfo = [&recoData](const o2::tpc::TrackTPC& trc, RecTrack& tref) {
285 if (recoData.inputsTPCclusters) {
286 uint8_t clSect = 0, clRow = 0, lowestR = -1;
287 uint32_t clIdx = 0;
288 const auto clRefs = recoData.getTPCTracksClusterRefs();
289 const auto tpcClusAcc = recoData.getTPCClusters();
290 const auto shMap = recoData.clusterShMapTPC;
291 for (int ic = 0; ic < trc.getNClusterReferences(); ic++) { // outside -> inside ordering, but on the sector boundaries backward jumps are possible
292 trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
293 if (clRow < lowestR) {
294 tref.rowCountTPC++;
295 lowestR = clRow;
296 }
297 unsigned int absoluteIndex = tpcClusAcc.clusterOffset[clSect][clRow] + clIdx;
298 if (shMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) {
299 tref.nClTPCShared++;
300 }
301 }
302 tref.lowestPadRow = lowestR;
303 const auto& clus = tpcClusAcc.clusters[clSect][clRow][clIdx];
304 int padFromEdge = int(clus.getPad()), npads = o2::gpu::GPUTPCGeometry::NPads(clRow);
305 if (padFromEdge > npads / 2) {
306 padFromEdge = npads - 1 - padFromEdge;
307 }
308 tref.padFromEdge = uint8_t(padFromEdge);
309 trc.getClusterReference(clRefs, 0, clSect, clRow, clIdx);
310 tref.rowMaxTPC = clRow;
311 }
312 };
313
314 auto flagTPCClusters = [&recoData](const o2::tpc::TrackTPC& trc, o2::MCCompLabel lbTrc) {
315 if (recoData.inputsTPCclusters) {
316 const auto clRefs = recoData.getTPCTracksClusterRefs();
317 const auto* TPCClMClab = recoData.inputsTPCclusters->clusterIndex.clustersMCTruth;
318 const auto& TPCClusterIdxStruct = recoData.inputsTPCclusters->clusterIndex;
319 for (int ic = 0; ic < trc.getNClusterReferences(); ic++) {
320 uint8_t clSect = 0, clRow = 0;
321 uint32_t clIdx = 0;
322 trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
323 auto labels = TPCClMClab->getLabels(clIdx + TPCClusterIdxStruct.clusterOffset[clSect][clRow]);
324 for (auto& lbl : labels) {
325 if (lbl == lbTrc) {
326 const_cast<o2::MCCompLabel&>(lbl).setFakeFlag(true); // actually, in this way we are flagging that this cluster was correctly attached
327 break;
328 }
329 }
330 }
331 }
332 };
333
334 {
335 const auto* digconst = mcReader.getDigitizationContext();
336 const auto& mcEvRecords = digconst->getEventRecords(false);
337 int ITSTimeBias = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameBiasInBC;
338 int ITSROFLen = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameLengthInBC;
339 unsigned int rofCount = 0;
340 const auto ITSClusROFRec = recoData.getITSClustersROFRecords();
341 for (const auto& mcIR : mcEvRecords) {
342 long tbc = mcIR.differenceInBC(recoData.startIR);
343 auto& mcVtx = mMCVtVec.emplace_back();
344 mcVtx.ts = tbc * o2::constants::lhc::LHCBunchSpacingMUS + mcIR.getTimeOffsetWrtBC() * 1e-3;
345 mcVtx.ID = mIntBC.size();
346 mIntBC.push_back(tbc);
347 int occBin = tbc / 8 * mNTPCOccBinLengthInv;
348 mTPCOcc.push_back(occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]));
349 // fill ITS occupancy
350 long gbc = mcIR.toLong();
351 while (rofCount < ITSClusROFRec.size()) {
352 long rofbcMin = ITSClusROFRec[rofCount].getBCData().toLong() + ITSTimeBias, rofbcMax = rofbcMin + ITSROFLen;
353 if (gbc < rofbcMin) { // IRs and ROFs are sorted, so this IR is prior of all ROFs
354 mITSOcc.push_back(0);
355 } else if (gbc < rofbcMax) {
356 mITSOcc.push_back(ITSClusROFRec[rofCount].getNEntries());
357 } else {
358 rofCount++; // test next ROF
359 continue;
360 }
361 break;
362 }
363 if (mNTPCOccBinLengthInv > 0.f) {
364 mcVtx.occTPCV.resize(params.nOccBinsDrift);
365 int grp = TMath::Max(1, TMath::Nint(params.nTBPerOccBin * mNTPCOccBinLengthInv));
366 for (int ib = 0; ib < params.nOccBinsDrift; ib++) {
367 float smb = 0;
368 int tbs = occBin + TMath::Nint(ib * params.nTBPerOccBin * mNTPCOccBinLengthInv);
369 for (int ig = 0; ig < grp; ig++) {
370 if (tbs >= 0 && tbs < int(mTBinClOccHist.size())) {
371 smb += mTBinClOccHist[tbs];
372 }
373 tbs++;
374 }
375 mcVtx.occTPCV[ib] = smb;
376 }
377 }
378 if (rofCount >= ITSClusROFRec.size()) {
379 mITSOcc.push_back(0); // IR after the last ROF
380 }
381 }
382 }
383 // collect interesting MC particle (tracks and parents)
384 int curSrcMC = 0, curEvMC = 0;
385 for (curSrcMC = 0; curSrcMC < (int)mcReader.getNSources(); curSrcMC++) {
386 if (mVerbose > 1) {
387 LOGP(info, "Source {}", curSrcMC);
388 }
389 int nev = mcReader.getNEvents(curSrcMC);
390 bool okAccVtx = true;
391 if (nev != (int)mMCVtVec.size()) {
392 LOGP(debug, "source {} has {} events while {} MC vertices were booked", curSrcMC, nev, mMCVtVec.size());
393 okAccVtx = false;
394 if (nev > (int)mMCVtVec.size()) { // QED
395 continue;
396 }
397 }
398 for (curEvMC = 0; curEvMC < nev; curEvMC++) {
399 if (mVerbose > 1) {
400 LOGP(info, "Event {}", curEvMC);
401 }
402 mCurrMCTracks = &mcReader.getTracks(curSrcMC, curEvMC);
403 const_cast<o2::dataformats::MCEventHeader&>(mcReader.getMCEventHeader(curSrcMC, curEvMC)).GetVertex(mCurrMCVertex);
404 if (okAccVtx) {
405 auto& pos = mMCVtVec[curEvMC].pos;
406 if (pos[2] < -999) {
407 pos[0] = mCurrMCVertex.X();
408 pos[1] = mCurrMCVertex.Y();
409 pos[2] = mCurrMCVertex.Z();
410 }
411 }
412 for (int itr = 0; itr < mCurrMCTracks->size(); itr++) {
413 processMCParticle(curSrcMC, curEvMC, itr);
414 }
415 }
416 }
417 if (mVerbose > 0) {
418 for (int id = 0; id < mNCheckDecays; id++) {
419 LOGP(info, "Decay PDG={} : {} entries", params.decayPDG[id], mDecaysMaps[id].size());
420 }
421 }
422
423 // add reconstruction info to MC particles. If MC particle was not selected before but was reconstrected, account MC info
424 mRecProcStage = true; // MC particles accepted only at this stage will be flagged
425 for (int iv = 0; iv < nv; iv++) {
426 if (mVerbose > 1) {
427 LOGP(info, "processing PV {} of {}", iv, nv);
428 }
429 o2::MCEventLabel pvLbl;
430 int pvID = -1;
431 if (iv < (int)pvvecLbl.size()) {
432 pvLbl = pvvecLbl[iv];
433 pvID = iv;
434 if (pvLbl.isSet() && pvLbl.getEventID() < mMCVtVec.size()) {
435 mMCVtVec[pvLbl.getEventID()].recVtx.emplace_back(RecPV{pvvec[iv], pvLbl});
436 }
437 }
438 const auto& vtref = vtxRefs[iv];
439 for (int is = GTrackID::NSources; is--;) {
440 DetID::mask_t dm = GTrackID::getSourceDetectorsMask(is);
441 if (!mTracksSrc[is] || !recoData.isTrackSourceLoaded(is) || !(dm[DetID::ITS] || dm[DetID::TPC])) {
442 continue;
443 }
444 int idMin = vtref.getFirstEntryOfSource(is), idMax = idMin + vtref.getEntriesOfSource(is);
445 for (int i = idMin; i < idMax; i++) {
446 auto vid = trackIndex[i];
447 const auto& trc = recoData.getTrackParam(vid);
448 if (trc.getPt() < params.minPt || std::abs(trc.getTgl()) > params.maxTgl) {
449 continue;
450 }
451 auto lbl = recoData.getTrackMCLabel(vid);
452 if (lbl.isValid()) {
453 lbl.setFakeFlag(false);
454 auto entry = mSelMCTracks.find(lbl);
455 if (entry == mSelMCTracks.end()) { // add the track which was not added during MC scan
456 if (lbl.getSourceID() != curSrcMC || lbl.getEventID() != curEvMC) {
457 curSrcMC = lbl.getSourceID();
458 curEvMC = lbl.getEventID();
459 mCurrMCTracks = &mcReader.getTracks(curSrcMC, curEvMC);
460 const_cast<o2::dataformats::MCEventHeader&>(mcReader.getMCEventHeader(curSrcMC, curEvMC)).GetVertex(mCurrMCVertex);
461 }
462 if (!acceptMCCharged((*mCurrMCTracks)[lbl.getTrackID()], lbl)) {
463 continue;
464 }
465 entry = mSelMCTracks.find(lbl);
466 }
467 auto& trackFamily = entry->second;
468 if (vid.isAmbiguous()) { // do not repeat ambiguous tracks
469 if (trackFamily.contains(vid)) {
470 continue;
471 }
472 }
473 auto& trf = trackFamily.recTracks.emplace_back();
474 trf.gid = vid; // account(iv, vid);
475 trf.pvID = pvID;
476 trf.pvLabel = pvLbl;
477 while (dm[DetID::ITS] && dm[DetID::TPC]) { // this track should have both ITS and TPC parts, if ITS was mismatched, fill it to its proper MC track slot
478 auto gidSet = recoData.getSingleDetectorRefs(vid);
479 if (!gidSet[GTrackID::ITS].isSourceSet()) {
480 break; // AB track, nothing to check
481 }
482 auto lblITS = recoData.getTrackMCLabel(gidSet[GTrackID::ITS]);
483 if (lblITS == trackFamily.mcTrackInfo.label) {
484 break; // correct match, no need for special treatment
485 }
486 const auto& trcITSF = recoData.getTrackParam(gidSet[GTrackID::ITS]);
487 if (trcITSF.getPt() < params.minPt || std::abs(trcITSF.getTgl()) > params.maxTgl) {
488 break; // ignore this track
489 }
490 auto entryOfFake = mSelMCTracks.find(lblITS);
491 if (entryOfFake == mSelMCTracks.end()) { // this MC track was not selected
492 break;
493 }
494 auto& trackFamilyOfFake = entryOfFake->second;
495 auto& trfOfFake = trackFamilyOfFake.recTracks.emplace_back();
496 trfOfFake.gid = gidSet[GTrackID::ITS]; // account(iv, vid);
497 break;
498 }
499 if (mVerbose > 1) {
500 LOGP(info, "Matched rec track {} to MC track {}", vid.asString(), entry->first.asString());
501 }
502 } else {
503 continue;
504 }
505 }
506 }
507 }
508
509 LOGP(info, "collected {} MC tracks", mSelMCTracks.size());
510 if (params.minTPCRefsToExtractClRes > 0 || params.storeTPCTrackRefs) { // prepare MC trackrefs for TPC
511 processTPCTrackRefs();
512 }
513
514 int mcnt = 0;
515 for (auto& entry : mSelMCTracks) {
516 auto& trackFam = entry.second;
517 auto& tracks = trackFam.recTracks;
518 mcnt++;
519 if (tracks.empty()) {
520 continue;
521 }
522 if (mVerbose > 1) {
523 LOGP(info, "Processing MC track#{} {} -> {} reconstructed tracks", mcnt - 1, entry.first.asString(), tracks.size());
524 }
525 // sort according to the gid complexity (in principle, should be already sorted due to the backwards loop over NSources above
526 std::sort(tracks.begin(), tracks.end(), [](const RecTrack& lhs, const RecTrack& rhs) {
527 const auto mskL = lhs.gid.getSourceDetectorsMask();
528 const auto mskR = rhs.gid.getSourceDetectorsMask();
529 bool itstpcL = mskL[DetID::ITS] && mskL[DetID::TPC], itstpcR = mskR[DetID::ITS] && mskR[DetID::TPC];
530 if (itstpcL && !itstpcR) { // to avoid TPC/TRD or TPC/TOF shadowing ITS/TPC
531 return true;
532 }
533 return lhs.gid.getSource() > rhs.gid.getSource();
534 });
535 if (params.storeTPCTrackRefs) {
536 auto rft = mSelTRefIdx.find(entry.first);
537 if (rft != mSelTRefIdx.end()) {
538 auto rfent = rft->second;
539 for (int irf = rfent.first; irf < rfent.second; irf++) {
540 trackFam.mcTrackInfo.trackRefsTPC.push_back(mSelTRefs[irf]);
541 }
542 }
543 }
544 // fill track params
545 int tcnt = 0;
546 for (auto& tref : tracks) {
547 if (tref.gid.isSourceSet()) {
548 auto gidSet = recoData.getSingleDetectorRefs(tref.gid);
549 tref.track = recoData.getTrackParam(tref.gid);
550 if (recoData.getTrackMCLabel(tref.gid).isFake()) {
551 tref.flags |= RecTrack::FakeGLO;
552 }
553 auto msk = tref.gid.getSourceDetectorsMask();
554 if (msk[DetID::ITS]) {
555 if (gidSet[GTrackID::ITS].isSourceSet()) { // has ITS track rather than AB tracklet
556 tref.pattITS = getITSPatt(gidSet[GTrackID::ITS], tref.nClITS);
557 if (trackFam.entITS < 0) {
558 trackFam.entITS = tcnt;
559 }
560 auto lblITS = recoData.getTrackMCLabel(gidSet[GTrackID::ITS]);
561 if (lblITS.isFake()) {
562 tref.flags |= RecTrack::FakeITS;
563 }
564 if (lblITS == trackFam.mcTrackInfo.label) {
565 trackFam.entITSFound = tcnt;
566 }
567 } else { // AB ITS tracklet
568 tref.pattITS = getITSPatt(gidSet[GTrackID::ITSAB], tref.nClITS);
569 if (recoData.getTrackMCLabel(gidSet[GTrackID::ITSAB]).isFake()) {
570 tref.flags |= RecTrack::FakeITS;
571 }
572 }
573 if (msk[DetID::TPC] && trackFam.entITSTPC < 0) { // has both ITS and TPC contribution
574 trackFam.entITSTPC = tcnt;
575 if (recoData.getTrackMCLabel(gidSet[GTrackID::ITSTPC]).isFake()) {
576 tref.flags |= RecTrack::FakeITSTPC;
577 }
578 }
579 }
580 if (msk[DetID::TPC]) {
581 const auto& trtpc = recoData.getTPCTrack(gidSet[GTrackID::TPC]);
582 tref.nClTPC = trtpc.getNClusters();
583 if (trtpc.hasBothSidesClusters()) {
584 tref.flags |= RecTrack::HASACSides;
585 }
586 fillTPCClusterInfo(trtpc, tref);
587 flagTPCClusters(trtpc, entry.first);
588 if (trackFam.entTPC < 0) {
589 trackFam.entTPC = tcnt;
590 trackFam.tpcT0 = trtpc.getTime0();
591 }
592 if (recoData.getTrackMCLabel(gidSet[GTrackID::TPC]).isFake()) {
593 tref.flags |= RecTrack::FakeTPC;
594 }
595 }
596 float ts = 0, terr = 0;
597 if (tref.gid.getSource() != GTrackID::ITS) {
598 recoData.getTrackTime(tref.gid, ts, terr);
599 tref.ts = timeEst{ts, terr};
600 } else {
601 const auto& itsBra = mITSROFBracket[mITSROF[tref.gid.getIndex()]];
602 tref.ts = timeEst{itsBra.mean(), itsBra.delta() * SQRT12Inv};
603 }
604 } else {
605 LOGP(info, "Invalid entry {} of {} getTrackMCLabel {}", tcnt, tracks.size(), tref.gid.asString());
606 }
607 tcnt++;
608 }
609 if (trackFam.entITS > -1 && trackFam.entTPC > -1) { // ITS and TPC were found but matching failed
610 auto vidITS = recoData.getITSContributorGID(tracks[trackFam.entITS].gid);
611 auto vidTPC = recoData.getTPCContributorGID(tracks[trackFam.entTPC].gid);
612 auto trcTPC = recoData.getTrackParam(vidTPC);
613 auto trcITS = recoData.getTrackParamOut(vidITS);
614 if (propagateToRefX(trcTPC, trcITS)) {
615 trackFam.trackITSProp = trcITS;
616 trackFam.trackTPCProp = trcTPC;
617 } else {
618 trackFam.trackITSProp.invalidate();
619 trackFam.trackTPCProp.invalidate();
620 }
621 } else {
622 trackFam.trackITSProp.invalidate();
623 trackFam.trackTPCProp.invalidate();
624 }
625 }
626
627 // SVertices (V0s)
628 if (mCheckSV) {
629 auto v0s = recoData.getV0sIdx();
630 auto prpr = [](o2::trackstudy::TrackFamily& f) {
631 std::string s;
632 s += fmt::format(" par {} Ntpccl={} Nitscl={} ", f.mcTrackInfo.pdgParent, f.mcTrackInfo.nTPCCl, f.mcTrackInfo.nITSCl);
633 for (auto& t : f.recTracks) {
634 s += t.gid.asString();
635 s += " ";
636 }
637 return s;
638 };
639 for (int svID; svID < (int)v0s.size(); svID++) {
640 const auto& v0idx = v0s[svID];
641 int nOKProngs = 0, realMCSVID = -1;
642 int8_t decTypeID = -1;
643 for (int ipr = 0; ipr < v0idx.getNProngs(); ipr++) {
644 auto mcl = recoData.getTrackMCLabel(v0idx.getProngID(ipr)); // was this MC particle selected?
645 auto itl = mSelMCTracks.find(mcl);
646 if (itl == mSelMCTracks.end()) {
647 nOKProngs = -1; // was not selected as interesting one, ignore
648 break;
649 }
650 auto& trackFamily = itl->second;
651 int decayParentIndex = trackFamily.mcTrackInfo.parentEntry;
652 if (decayParentIndex < 0) { // does not come from decay
653 break;
654 }
655 if (ipr == 0) {
656 realMCSVID = decayParentIndex;
657 decTypeID = trackFamily.mcTrackInfo.parentDecID;
658 nOKProngs = 1;
659 LOGP(debug, "Prong{} {} comes from {}/{}", ipr, prpr(trackFamily), decTypeID, realMCSVID);
660 continue;
661 }
662 if (realMCSVID != decayParentIndex || decTypeID != trackFamily.mcTrackInfo.parentDecID) {
663 break;
664 }
665 LOGP(debug, "Prong{} {} comes from {}/{}", ipr, prpr(trackFamily), decTypeID, realMCSVID);
666 nOKProngs++;
667 }
668 if (nOKProngs == v0idx.getNProngs()) { // all prongs are from the decay of MC parent which deemed to be interesting, flag it
669 LOGP(debug, "Decay {}/{} was found", decTypeID, realMCSVID);
670 mDecaysMaps[decTypeID][realMCSVID].foundSVID = svID;
671 }
672 }
673 }
674
675 // collect ITS/TPC cluster info for selected MC particles
676 fillMCClusterInfo(recoData);
677
678 // single tracks
679 for (auto& entry : mSelMCTracks) {
680 auto& trackFam = entry.second;
681 (*mDBGOut) << "tracks" << "tr=" << trackFam << "\n";
682 }
683
684 // decays
685 std::vector<TrackFamily> decFam;
686 for (int id = 0; id < mNCheckDecays; id++) {
687 std::string decTreeName = fmt::format("dec{}", params.decayPDG[id]);
688 for (const auto& dec : mDecaysMaps[id]) {
689 decFam.clear();
690 bool skip = false;
691 for (int idd = dec.daughterFirst; idd <= dec.daughterLast; idd++) {
692 auto dtLbl = mDecProdLblPool[idd]; // daughter MC label
693 const auto& dtFamily = mSelMCTracks[dtLbl];
694 if (dtFamily.mcTrackInfo.pdgParent != dec.pdg) {
695 LOGP(error, "{}-th decay (pdg={}): {} in {}:{} range refers to MC track with pdgParent = {}", id, params.decayPDG[id], idd, dec.daughterFirst, dec.daughterLast, dtFamily.mcTrackInfo.pdgParent);
696 skip = true;
697 break;
698 }
699 decFam.push_back(dtFamily);
700 }
701 if (!skip) {
703 if (dec.foundSVID >= 0 && !refitV0(dec.foundSVID, v0, recoData)) {
704 v0.invalidate();
705 }
706 (*mDBGOut) << decTreeName.c_str() << "pdgPar=" << dec.pdg << "trPar=" << dec.parent << "prod=" << decFam << "found=" << dec.foundSVID << "sv=" << v0 << "\n";
707 }
708 }
709 }
710
711 for (auto& mcVtx : mMCVtVec) { // sort rec.vertices in mult. order
712 std::sort(mcVtx.recVtx.begin(), mcVtx.recVtx.end(), [](const RecPV& lhs, const RecPV& rhs) {
713 return lhs.pv.getNContributors() > rhs.pv.getNContributors();
714 });
715 (*mDBGOut) << "mcVtxTree" << "mcVtx=" << mcVtx << "\n";
716 }
717}
718
719void TrackMCStudy::processTPCTrackRefs()
720{
721 constexpr float alpsec[18] = {0.174533, 0.523599, 0.872665, 1.221730, 1.570796, 1.919862, 2.268928, 2.617994, 2.967060, 3.316126, 3.665191, 4.014257, 4.363323, 4.712389, 5.061455, 5.410521, 5.759587, 6.108652};
722 constexpr float sinAlp[18] = {0.173648, 0.500000, 0.766044, 0.939693, 1.000000, 0.939693, 0.766044, 0.500000, 0.173648, -0.173648, -0.500000, -0.766044, -0.939693, -1.000000, -0.939693, -0.766044, -0.500000, -0.173648};
723 constexpr float cosAlp[18] = {0.984808, 0.866025, 0.642788, 0.342020, 0.000000, -0.342020, -0.642788, -0.866025, -0.984808, -0.984808, -0.866025, -0.642788, -0.342020, -0.000000, 0.342020, 0.642788, 0.866025, 0.984808};
725 for (auto& entry : mSelMCTracks) {
726 auto lb = entry.first;
727 auto trspan = mcReader.getTrackRefs(lb.getSourceID(), lb.getEventID(), lb.getTrackID());
728 int q = entry.second.mcTrackInfo.track.getCharge();
729 if (q * q != 1) {
730 continue;
731 }
732 int ref0entry = mSelTRefs.size(), nrefsSel = 0;
733 for (const auto& trf : trspan) {
734 if (trf.getDetectorId() != 1) { // process TPC only
735 continue;
736 }
737 float pT = std::sqrt(trf.Px() * trf.Px() + trf.Py() * trf.Py());
738 if (pT < 0.05) {
739 continue;
740 }
741 float secX, secY, phi = std::atan2(trf.Y(), trf.X());
742 int sector = o2::math_utils::angle2Sector(phi);
743 o2::math_utils::rotateZInv(trf.X(), trf.Y(), secX, secY, sinAlp[sector], cosAlp[sector]); // sector coordinates
744 float phiPt = std::atan2(trf.Py(), trf.Px());
745 o2::math_utils::bringTo02Pi(phiPt);
746 auto dphiPt = phiPt - alpsec[sector];
747 if (dphiPt > o2::constants::math::PI) { // account for wraps
749 } else if (dphiPt < -o2::constants::math::PI) {
751 } else if (std::abs(dphiPt) > o2::constants::math::PIHalf * 0.8) {
752 continue; // ignore backward going or parallel to padrows tracks
753 }
754 float tgL = trf.Pz() / pT;
755 std::array<float, 5> pars = {secY, trf.Z(), std::sin(dphiPt), tgL, q / pT};
756 auto& refTrack = mSelTRefs.emplace_back(secX, alpsec[sector], pars);
757 refTrack.setUserField(uint16_t(sector));
758 nrefsSel++;
759 }
760 if (nrefsSel < params.minTPCRefsToExtractClRes) {
761 mSelTRefs.resize(ref0entry); // discard unused tracks
762 continue;
763 } else {
764 mSelTRefIdx[lb] = std::make_pair(ref0entry, ref0entry + nrefsSel);
765 }
766 }
767}
768
769void TrackMCStudy::fillMCClusterInfo(const o2::globaltracking::RecoContainer& recoData)
770{
771 // TPC clusters info
772 const auto& TPCClusterIdxStruct = recoData.inputsTPCclusters->clusterIndex;
773 const auto* TPCClMClab = recoData.inputsTPCclusters->clusterIndex.clustersMCTruth;
775
776 ClResTPC clRes{};
777 for (uint8_t row = 0; row < 152; row++) { // we need to go in increasing row, so this should be the outer loop
778 for (uint8_t sector = 0; sector < 36; sector++) {
779 unsigned int offs = TPCClusterIdxStruct.clusterOffset[sector][row];
780 for (unsigned int icl0 = 0; icl0 < TPCClusterIdxStruct.nClusters[sector][row]; icl0++) {
781 const auto labels = TPCClMClab->getLabels(icl0 + offs);
782 int ncontLb = 0; // number of real contrubutors to this label (w/o noise)
783 for (const auto& lbl : labels) {
784 if (!lbl.isValid()) {
785 continue;
786 }
787 ncontLb++;
788 }
789 const auto& clus = TPCClusterIdxStruct.clusters[sector][row][icl0];
790 int tbinH = int(clus.getTime() * mNTPCOccBinLengthInv); // time bin converted to slot of the occ. histo
791 clRes.contTracks.clear();
792 bool doClusRes = (params.minTPCRefsToExtractClRes > 0) && (params.rejectClustersResStat <= 0. || gRandom->Rndm() < params.rejectClustersResStat);
793 for (auto lbl : labels) {
794 bool corrAttach = lbl.isFake(); // was this flagged in the flagTPCClusters called from process ?
795 lbl.setFakeFlag(false);
796 auto entry = mSelMCTracks.find(lbl);
797 if (entry == mSelMCTracks.end()) { // not selected
798 continue;
799 }
800 auto& mctr = entry->second.mcTrackInfo;
801 mctr.nTPCCl++;
802 if (row > mctr.maxTPCRow) {
803 mctr.maxTPCRow = row;
804 mctr.maxTPCRowSect = sector;
805 mctr.nUsedPadRows++;
806 } else if (row == 0 && mctr.nUsedPadRows == 0) {
807 mctr.nUsedPadRows++;
808 }
809 if (row < mctr.minTPCRow) {
810 mctr.minTPCRow = row;
811 mctr.minTPCRowSect = sector;
812 }
813 if (mctr.minTPCRowSect == sector && row > mctr.maxTPCRowInner) {
814 mctr.maxTPCRowInner = row;
815 }
816 if (ncontLb > 1) {
817 mctr.nTPCClShared++;
818 }
819 // try to extract ideal track position
820 if (doClusRes) {
821 auto entTRefIDsIt = mSelTRefIdx.find(lbl);
822 if (entTRefIDsIt == mSelTRefIdx.end()) {
823 continue;
824 }
825 float xc, yc, zc;
826 mTPCCorrMapsLoader.Transform(sector, row, clus.getPad(), clus.getTime(), xc, yc, zc, mctr.bcInTF / 8.); // nominal time of the track
827
828 const auto& entTRefIDs = entTRefIDsIt->second;
829 // find bracketing TRef params
830 int entIDBelow = -1, entIDAbove = -1;
831 float xBelow = -1e6, xAbove = 1e6;
832
833 for (int entID = entTRefIDs.first; entID < entTRefIDs.second; entID++) {
834 const auto& refTr = mSelTRefs[entID];
835 if (refTr.getUserField() != sector % 18) {
836 continue;
837 }
838 if ((refTr.getX() < xc) && (refTr.getX() > xBelow) && (refTr.getX() > xc - params.maxTPCRefExtrap)) {
839 xBelow = refTr.getX();
840 entIDBelow = entID;
841 }
842 if ((refTr.getX() > xc) && (refTr.getX() < xAbove) && (refTr.getX() < xc + params.maxTPCRefExtrap)) {
843 xAbove = refTr.getX();
844 entIDAbove = entID;
845 }
846 }
847 if ((entIDBelow < 0 && entIDAbove < 0) || (params.requireTopBottomRefs && (entIDBelow < 0 || entIDAbove < 0))) {
848 continue;
849 }
850 auto prop = o2::base::Propagator::Instance();
851 o2::track::TrackPar tparAbove, tparBelow;
852 bool okBelow = entIDBelow >= 0 && prop->PropagateToXBxByBz((tparBelow = mSelTRefs[entIDBelow]), xc, 0.99, 2.);
853 bool okAbove = entIDAbove >= 0 && prop->PropagateToXBxByBz((tparAbove = mSelTRefs[entIDAbove]), xc, 0.99, 2.);
854 if ((!okBelow && !okAbove) || (params.requireTopBottomRefs && (!okBelow || !okAbove))) {
855 continue;
856 }
857
858 int nmeas = 0;
859 auto& clCont = clRes.contTracks.emplace_back();
860 clCont.corrAttach = corrAttach;
861 if (okBelow) {
862 clCont.below = {mSelTRefs[entIDBelow].getX(), tparBelow.getY(), tparBelow.getZ()};
863 clCont.snp += tparBelow.getSnp();
864 clCont.tgl += tparBelow.getTgl();
865 clCont.q2pt += tparBelow.getQ2Pt();
866 nmeas++;
867 }
868 if (okAbove) {
869 clCont.above = {mSelTRefs[entIDAbove].getX(), tparAbove.getY(), tparAbove.getZ()};
870 clCont.snp += tparAbove.getSnp();
871 clCont.tgl += tparAbove.getTgl();
872 clCont.q2pt += tparAbove.getQ2Pt();
873 nmeas++;
874 }
875 if (nmeas) {
876 if (clRes.contTracks.size() == 1) {
877 int occBin = mctr.bcInTF / 8 * mNTPCOccBinLengthInv;
878 clRes.occ = occBin < 0 ? mTBinClOcc[0] : (occBin >= mTBinClOcc.size() ? mTBinClOcc.back() : mTBinClOcc[occBin]);
879 }
880 clCont.xyz = {xc, yc, zc};
881 if (nmeas > 1) {
882 clCont.snp *= 0.5;
883 clCont.tgl *= 0.5;
884 clCont.q2pt *= 0.5;
885 }
886 } else {
887 clRes.contTracks.pop_back();
888 }
889 }
890 }
891 if (clRes.getNCont()) {
892 clRes.sect = sector;
893 clRes.row = row;
894 clRes.qtot = clus.getQtot();
895 clRes.qmax = clus.getQmax();
896 clRes.flags = clus.getFlags();
897 clRes.sigmaTimePacked = clus.sigmaTimePacked;
898 clRes.sigmaPadPacked = clus.sigmaPadPacked;
899 clRes.ncont = ncontLb;
900 clRes.sortCont();
901
902 if (tbinH < 0) {
903 tbinH = 0;
904 } else if (tbinH >= int(mTBinClOccHist.size())) {
905 tbinH = (int)mTBinClOccHist.size() - 1;
906 }
907 clRes.occBin = mTBinClOccHist[tbinH];
908
909 (*mDBGOut) << "clres" << "clr=" << clRes << "\n";
910 }
911 }
912 }
913 }
914 // fill ITS cluster info
915 const auto* mcITSClusters = recoData.getITSClustersMCLabels();
916 const auto& ITSClusters = recoData.getITSClusters();
917 for (unsigned int icl = 0; icl < ITSClusters.size(); icl++) {
918 const auto labels = mcITSClusters->getLabels(icl);
919 for (const auto& lbl : labels) {
920 auto entry = mSelMCTracks.find(lbl);
921 if (entry == mSelMCTracks.end()) { // not selected
922 continue;
923 }
924 auto& mctr = entry->second.mcTrackInfo;
925 mctr.nITSCl++;
926 mctr.pattITSCl |= 0x1 << o2::itsmft::ChipMappingITS::getLayer(ITSClusters[icl].getChipID());
927 }
928 }
929}
930
931bool TrackMCStudy::propagateToRefX(o2::track::TrackParCov& trcTPC, o2::track::TrackParCov& trcITS)
932{
933 bool refReached = false;
934 constexpr float TgHalfSector = 0.17632698f;
936 int trialsLeft = 2;
937 while (o2::base::Propagator::Instance()->PropagateToXBxByBz(trcTPC, par.XMatchingRef, MaxSnp, 2., par.matCorr)) {
938 if (refReached) {
939 break;
940 }
941 // make sure the track is indeed within the sector defined by alpha
942 if (fabs(trcTPC.getY()) < par.XMatchingRef * TgHalfSector) {
943 refReached = true;
944 break; // ok, within
945 }
946 if (!trialsLeft--) {
947 break;
948 }
949 auto alphaNew = o2::math_utils::angle2Alpha(trcTPC.getPhiPos());
950 if (!trcTPC.rotate(alphaNew) != 0) {
951 break; // failed (RS: check effect on matching tracks to neighbouring sector)
952 }
953 }
954 if (!refReached) {
955 return false;
956 }
957 refReached = false;
958 float alp = trcTPC.getAlpha();
959 if (!trcITS.rotate(alp) != 0 || !o2::base::Propagator::Instance()->PropagateToXBxByBz(trcITS, par.XMatchingRef, MaxSnp, 2., par.matCorr)) {
960 return false;
961 }
962 return true;
963}
964
965void TrackMCStudy::endOfStream(EndOfStreamContext& ec)
966{
967 mDBGOut.reset();
968}
969
970void TrackMCStudy::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj)
971{
972 if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) {
973 return;
974 }
975 if (mTPCVDriftHelper.accountCCDBInputs(matcher, obj)) {
976 return;
977 }
978 if (mTPCCorrMapsLoader.accountCCDBInputs(matcher, obj)) {
979 return;
980 }
981 if (matcher == ConcreteDataMatcher("ITS", "ALPIDEPARAM", 0)) {
982 LOG(info) << "ITS Alpide param updated";
984 par.printKeyValues();
985 mITSTimeBiasMUS = par.roFrameBiasInBC * o2::constants::lhc::LHCBunchSpacingNS * 1e-3;
986 mITSROFrameLengthMUS = par.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS * 1e-3;
987 return;
988 }
989}
990
991//_____________________________________________________
992void TrackMCStudy::prepareITSData(const o2::globaltracking::RecoContainer& recoData)
993{
994 const auto ITSTracksArray = recoData.getITSTracks();
995 const auto ITSTrackROFRec = recoData.getITSTracksROFRecords();
996 int nROFs = ITSTrackROFRec.size();
997 mITSROF.clear();
998 mITSROFBracket.clear();
999 mITSROF.reserve(ITSTracksArray.size());
1000 mITSROFBracket.reserve(ITSTracksArray.size());
1001 for (int irof = 0; irof < nROFs; irof++) {
1002 const auto& rofRec = ITSTrackROFRec[irof];
1003 long nBC = rofRec.getBCData().differenceInBC(recoData.startIR);
1004 float tMin = nBC * o2::constants::lhc::LHCBunchSpacingMUS + mITSTimeBiasMUS;
1005 float tMax = tMin + mITSROFrameLengthMUS;
1006 mITSROFBracket.emplace_back(tMin, tMax);
1007 for (int it = 0; it < rofRec.getNEntries(); it++) {
1008 mITSROF.push_back(irof);
1009 }
1010 }
1011}
1012/*
1013float TrackMCStudy::getDCAYCut(float pt) const
1014{
1015 static TF1 fun("dcayvspt", mDCAYFormula.c_str(), 0, 20);
1016 return fun.Eval(pt);
1017}
1018*/
1019
1020bool TrackMCStudy::processMCParticle(int src, int ev, int trid)
1021{
1022 const auto& mcPart = (*mCurrMCTracks)[trid];
1023 int pdg = mcPart.GetPdgCode();
1024 bool res = false;
1025 while (true) {
1026 auto lbl = o2::MCCompLabel(trid, ev, src);
1027 int decay = -1; // is this decay to watch?
1029 if (mcPart.T() < params.decayMotherMaxT) {
1030 for (int id = 0; id < mNCheckDecays; id++) {
1031 if (params.decayPDG[id] == std::abs(pdg)) {
1032 decay = id;
1033 break;
1034 }
1035 }
1036 if (decay >= 0) { // check if decay and kinematics is acceptable
1037 auto& decayPool = mDecaysMaps[decay];
1038 int idd0 = mcPart.getFirstDaughterTrackId(), idd1 = mcPart.getLastDaughterTrackId(); // we want only charged and trackable daughters
1039 int dtStart = mDecProdLblPool.size(), dtEnd = -1;
1040 if (idd0 < 0) {
1041 break;
1042 }
1043 for (int idd = idd0; idd <= idd1; idd++) {
1044 const auto& product = (*mCurrMCTracks)[idd];
1045 auto lbld = o2::MCCompLabel(idd, ev, src);
1046 if (!acceptMCCharged(product, lbld, decay)) {
1047 decay = -1; // discard decay
1048 mDecProdLblPool.resize(dtStart);
1049 break;
1050 }
1051 mDecProdLblPool.push_back(lbld); // register prong entry and label
1052 }
1053 if (decay >= 0) {
1054 // account decay
1055 dtEnd = mDecProdLblPool.size();
1056 for (int dtid = dtStart; dtid < dtEnd; dtid++) { // flag selected decay parent entry in the prongs MCs
1057 mSelMCTracks[mDecProdLblPool[dtid]].mcTrackInfo.parentEntry = decayPool.size();
1058 mSelMCTracks[mDecProdLblPool[dtid]].mcTrackInfo.parentDecID = int8_t(decay);
1059 }
1060 dtEnd--;
1061 std::array<float, 3> xyz{(float)mcPart.GetStartVertexCoordinatesX(), (float)mcPart.GetStartVertexCoordinatesY(), (float)mcPart.GetStartVertexCoordinatesZ()};
1062 std::array<float, 3> pxyz{(float)mcPart.GetStartVertexMomentumX(), (float)mcPart.GetStartVertexMomentumY(), (float)mcPart.GetStartVertexMomentumZ()};
1063 decayPool.emplace_back(DecayRef{lbl,
1064 o2::track::TrackPar(xyz, pxyz, TMath::Nint(O2DatabasePDG::Instance()->GetParticle(mcPart.GetPdgCode())->Charge() / 3), false),
1065 mcPart.GetPdgCode(), dtStart, dtEnd});
1066 if (mVerbose > 1) {
1067 LOGP(info, "Adding MC parent pdg={} {}, with prongs in {}:{} range", pdg, lbl.asString(), dtStart, dtEnd);
1068 }
1069 res = true; // Accept!
1070 }
1071 break;
1072 }
1073 }
1074 // check if this is a charged which should be processed but was not accounted as a decay product
1075 if (mSelMCTracks.find(lbl) == mSelMCTracks.end()) {
1076 res = acceptMCCharged(mcPart, lbl);
1077 }
1078 break;
1079 }
1080 return res;
1081}
1082
1083bool TrackMCStudy::acceptMCCharged(const MCTrack& tr, const o2::MCCompLabel& lb, int followDecay)
1084{
1086 if (tr.GetPt() < params.minPtMC ||
1087 std::abs(tr.GetTgl()) > params.maxTglMC ||
1088 tr.R2() > params.maxRMC * params.maxRMC) {
1089 if (mVerbose > 1 && followDecay > -1) {
1090 LOGP(info, "rejecting decay {} prong : pdg={}, pT={}, tgL={}, r={}", followDecay, tr.GetPdgCode(), tr.GetPt(), tr.GetTgl(), std::sqrt(tr.R2()));
1091 }
1092 return false;
1093 }
1094 float dx = tr.GetStartVertexCoordinatesX() - mCurrMCVertex.X(), dy = tr.GetStartVertexCoordinatesY() - mCurrMCVertex.Y(), dz = tr.GetStartVertexCoordinatesZ() - mCurrMCVertex.Z();
1095 float r2 = dx * dx + dy * dy;
1096 float posTgl2 = r2 > 1 && std::abs(dz) < 20 ? dz * dz / r2 : 0;
1097 if (posTgl2 > params.maxPosTglMC * params.maxPosTglMC) {
1098 if (mVerbose > 1 && followDecay > -1) {
1099 LOGP(info, "rejecting decay {} prong : pdg={}, pT={}, tgL={}, dr={}, dz={} r={}, z={}, posTgl={}", followDecay, tr.GetPdgCode(), tr.GetPt(), tr.GetTgl(), std::sqrt(r2), dz, std::sqrt(tr.R2()), tr.GetStartVertexCoordinatesZ(), std::sqrt(posTgl2));
1100 }
1101 return false;
1102 }
1103 if (params.requireITSorTPCTrackRefs) {
1104 auto trspan = mcReader.getTrackRefs(lb.getSourceID(), lb.getEventID(), lb.getTrackID());
1105 bool ok = false;
1106 for (const auto& trf : trspan) {
1107 if (trf.getDetectorId() == DetID::ITS || trf.getDetectorId() == DetID::TPC) {
1108 ok = true;
1109 break;
1110 }
1111 }
1112 if (!ok) {
1113 return false;
1114 }
1115 }
1116 TParticlePDG* pPDG = O2DatabasePDG::Instance()->GetParticle(tr.GetPdgCode());
1117 if (!pPDG) {
1118 LOGP(debug, "Unknown particle {}", tr.GetPdgCode());
1119 return false;
1120 }
1121 if (pPDG->Charge() == 0.) {
1122 return false;
1123 }
1124 return addMCParticle(tr, lb, pPDG);
1125}
1126
1127bool TrackMCStudy::addMCParticle(const MCTrack& mcPart, const o2::MCCompLabel& lb, TParticlePDG* pPDG)
1128{
1129 std::array<float, 3> xyz{(float)mcPart.GetStartVertexCoordinatesX(), (float)mcPart.GetStartVertexCoordinatesY(), (float)mcPart.GetStartVertexCoordinatesZ()};
1130 std::array<float, 3> pxyz{(float)mcPart.GetStartVertexMomentumX(), (float)mcPart.GetStartVertexMomentumY(), (float)mcPart.GetStartVertexMomentumZ()};
1131 if (!pPDG && !(pPDG = O2DatabasePDG::Instance()->GetParticle(mcPart.GetPdgCode()))) {
1132 LOGP(debug, "Unknown particle {}", mcPart.GetPdgCode());
1133 return false;
1134 }
1135 auto& mcEntry = mSelMCTracks[lb];
1136 mcEntry.mcTrackInfo.pdg = mcPart.GetPdgCode();
1137 mcEntry.mcTrackInfo.track = o2::track::TrackPar(xyz, pxyz, TMath::Nint(pPDG->Charge() / 3), true);
1138 mcEntry.mcTrackInfo.label = lb;
1139 mcEntry.mcTrackInfo.bcInTF = mIntBC[lb.getEventID()];
1140 mcEntry.mcTrackInfo.occTPC = mTPCOcc[lb.getEventID()];
1141 mcEntry.mcTrackInfo.occITS = mITSOcc[lb.getEventID()];
1142 mcEntry.mcTrackInfo.occTPCV = mMCVtVec[lb.getEventID()].occTPCV;
1143 if (mRecProcStage) {
1144 mcEntry.mcTrackInfo.setAddedAtRecStage();
1145 }
1146 if (o2::mcutils::MCTrackNavigator::isPhysicalPrimary(mcPart, *mCurrMCTracks)) {
1147 mcEntry.mcTrackInfo.setPrimary();
1148 }
1149 int moth = -1;
1150 o2::MCCompLabel mclbPar;
1151 if ((moth = mcPart.getMotherTrackId()) >= 0) {
1152 const auto& mcPartPar = (*mCurrMCTracks)[moth];
1153 mcEntry.mcTrackInfo.pdgParent = mcPartPar.GetPdgCode();
1154 }
1155 if (mcPart.isPrimary() && mcReader.getNEvents(lb.getSourceID()) == mMCVtVec.size()) {
1156 mMCVtVec[lb.getEventID()].nTrackSel++;
1157 }
1158 if (mVerbose > 1) {
1159 LOGP(info, "Adding charged MC pdg={} {} ", mcPart.GetPdgCode(), lb.asString());
1160 }
1161 return true;
1162}
1163
1164bool TrackMCStudy::refitV0(int i, o2::dataformats::V0& v0, const o2::globaltracking::RecoContainer& recoData)
1165{
1166 const auto& id = recoData.getV0sIdx()[i];
1167 auto seedP = recoData.getTrackParam(id.getProngID(0));
1168 auto seedN = recoData.getTrackParam(id.getProngID(1));
1169 bool isTPConly = (id.getProngID(0).getSource() == GTrackID::TPC) || (id.getProngID(1).getSource() == GTrackID::TPC);
1170 const auto& svparam = o2::vertexing::SVertexerParams::Instance();
1171 if (svparam.mTPCTrackPhotonTune && isTPConly) {
1172 mFitterV0.setMaxDZIni(svparam.mTPCTrackMaxDZIni);
1173 mFitterV0.setMaxDXYIni(svparam.mTPCTrackMaxDXYIni);
1174 mFitterV0.setMaxChi2(svparam.mTPCTrackMaxChi2);
1175 mFitterV0.setCollinear(true);
1176 }
1177 int nCand = mFitterV0.process(seedP, seedN);
1178 if (svparam.mTPCTrackPhotonTune && isTPConly) { // restore
1179 // Reset immediately to the defaults
1180 mFitterV0.setMaxDZIni(svparam.maxDZIni);
1181 mFitterV0.setMaxDXYIni(svparam.maxDXYIni);
1182 mFitterV0.setMaxChi2(svparam.maxChi2);
1183 mFitterV0.setCollinear(false);
1184 }
1185 if (nCand == 0) { // discard this pair
1186 return false;
1187 }
1188 const int cand = 0;
1189 if (!mFitterV0.isPropagateTracksToVertexDone(cand) && !mFitterV0.propagateTracksToVertex(cand)) {
1190 return false;
1191 }
1192 const auto& trPProp = mFitterV0.getTrack(0, cand);
1193 const auto& trNProp = mFitterV0.getTrack(1, cand);
1194 std::array<float, 3> pP{}, pN{};
1195 trPProp.getPxPyPzGlo(pP);
1196 trNProp.getPxPyPzGlo(pN);
1197 std::array<float, 3> pV0 = {pP[0] + pN[0], pP[1] + pN[1], pP[2] + pN[2]};
1198 auto p2V0 = pV0[0] * pV0[0] + pV0[1] * pV0[1] + pV0[2] * pV0[2];
1199 const auto& pv = recoData.getPrimaryVertex(id.getVertexID());
1200 const auto v0XYZ = mFitterV0.getPCACandidatePos(cand);
1201 float dx = v0XYZ[0] - pv.getX(), dy = v0XYZ[1] - pv.getY(), dz = v0XYZ[2] - pv.getZ(), prodXYZv0 = dx * pV0[0] + dy * pV0[1] + dz * pV0[2];
1202 float cosPA = prodXYZv0 / std::sqrt((dx * dx + dy * dy + dz * dz) * p2V0);
1203 new (&v0) o2::dataformats::V0(v0XYZ, pV0, mFitterV0.calcPCACovMatrixFlat(cand), trPProp, trNProp);
1204 v0.setDCA(mFitterV0.getChi2AtPCACandidate(cand));
1205 v0.setCosPA(cosPA);
1206 return true;
1207}
1208
1209void TrackMCStudy::loadTPCOccMap(const o2::globaltracking::RecoContainer& recoData)
1210{
1211 auto NHBPerTF = o2::base::GRPGeomHelper::instance().getGRPECS()->getNHBFPerTF();
1212 const auto& TPCOccMap = recoData.occupancyMapTPC;
1213 auto prop = o2::base::Propagator::Instance();
1214 auto TPCRefitter = std::make_unique<o2::gpu::GPUO2InterfaceRefit>(&recoData.inputsTPCclusters->clusterIndex, &mTPCCorrMapsLoader, prop->getNominalBz(),
1215 recoData.getTPCTracksClusterRefs().data(), 0, recoData.clusterShMapTPC.data(), TPCOccMap.data(), TPCOccMap.size(), nullptr, prop);
1216 mNTPCOccBinLength = TPCRefitter->getParam()->rec.tpc.occupancyMapTimeBins;
1217 mTBinClOcc.clear();
1218 if (mNTPCOccBinLength > 1 && TPCOccMap.size()) {
1219 mNTPCOccBinLengthInv = 1. / mNTPCOccBinLength;
1220 int nTPCBins = NHBPerTF * o2::constants::lhc::LHCMaxBunches / 8, ninteg = 0;
1221 int nTPCOccBins = nTPCBins * mNTPCOccBinLengthInv, sumBins = std::max(1, int(o2::constants::lhc::LHCMaxBunches / 8 * mNTPCOccBinLengthInv));
1222 mTBinClOcc.resize(nTPCOccBins);
1223 mTBinClOccHist.resize(nTPCOccBins);
1224 float sm = 0., tb = 0.5 * mNTPCOccBinLength;
1225 for (int i = 0; i < nTPCOccBins; i++) {
1226 mTBinClOccHist[i] = TPCRefitter->getParam()->GetUnscaledMult(tb);
1227 tb += mNTPCOccBinLength;
1228 }
1229 for (int i = nTPCOccBins; i--;) {
1230 sm += mTBinClOccHist[i];
1231 if (i + sumBins < nTPCOccBins) {
1232 sm -= mTBinClOccHist[i + sumBins];
1233 }
1234 mTBinClOcc[i] = sm;
1235 }
1236 } else {
1237 mTBinClOcc.resize(1);
1238 mTBinClOccHist.resize(1);
1239 }
1240}
1241
1243{
1244 std::vector<OutputSpec> outputs;
1245 Options opts{
1246 {"device-verbosity", VariantType::Int, 0, {"Verbosity level"}},
1247 {"dcay-vs-pt", VariantType::String, "0.0105 + 0.0350 / pow(x, 1.1)", {"Formula for global tracks DCAy vs pT cut"}},
1248 {"min-tpc-clusters", VariantType::Int, 60, {"Cut on TPC clusters"}},
1249 {"max-tpc-dcay", VariantType::Float, 2.f, {"Cut on TPC dcaY"}},
1250 {"max-tpc-dcaz", VariantType::Float, 2.f, {"Cut on TPC dcaZ"}},
1251 {"min-x-prop", VariantType::Float, 6.f, {"track should be propagated to this X at least"}}};
1252 auto dataRequest = std::make_shared<DataRequest>();
1253 bool useMC = true;
1254 dataRequest->requestTracks(srcTracks, useMC);
1255 dataRequest->requestClusters(srcClusters, useMC);
1256 dataRequest->requestPrimaryVertices(useMC);
1257 if (checkSV) {
1258 dataRequest->requestSecondaryVertices(useMC);
1259 }
1260 o2::tpc::VDriftHelper::requestCCDBInputs(dataRequest->inputs);
1261 o2::tpc::CorrectionMapsLoader::requestCCDBInputs(dataRequest->inputs, opts, sclOpts);
1262 auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(false, // orbitResetTime
1263 true, // GRPECS=true
1264 true, // GRPLHCIF
1265 true, // GRPMagField
1266 true, // askMatLUT
1268 dataRequest->inputs,
1269 true);
1270
1271 return DataProcessorSpec{
1272 "track-mc-study",
1273 dataRequest->inputs,
1274 outputs,
1275 AlgorithmSpec{adaptFromTask<TrackMCStudy>(dataRequest, ggRequest, srcTracks, sclOpts, checkSV)},
1276 opts};
1277}
1278
1279} // namespace o2::trackstudy
Helper class to access load maps from CCDB.
Defintions for N-prongs secondary vertex fit.
Wrapper container for different reconstructed object types.
Definition of the GeometryManager class.
Definition of the FIT RecPoints class.
int32_t i
Helper for geometry and GRP related CCDB requests.
Global index for barrel track: provides provenance (detectors combination), index in respective array...
Utility functions for MC particles.
Configurable params for TPC ITS matching.
Definition of the Names Generator class.
Definition of the parameter class for the detector electronics.
uint16_t pos
Definition RawData.h:3
uint32_t res
Definition RawData.h:0
Wrapper container for different reconstructed object types.
o2::track::TrackParCov TrackParCov
Definition Recon.h:39
Configurable params for secondary vertexer.
Result of refitting TPC-ITS matched track.
Reference on ITS/MFT clusters set.
Helper class to extract VDrift from different sources.
std::ostringstream debug
Referenc on track indices contributing to the vertex, with possibility chose tracks from specific sou...
bool isFake() const
Definition MCCompLabel.h:78
void setFakeFlag(bool v=true)
int getTrackID() const
int getSourceID() const
int getEventID() const
std::string asString() const
bool isSet() const
int getEventID() const
Double_t GetStartVertexMomentumZ() const
Definition MCTrack.h:81
Double_t GetStartVertexMomentumX() const
Definition MCTrack.h:79
bool isPrimary() const
Definition MCTrack.h:75
Double_t GetStartVertexCoordinatesY() const
Definition MCTrack.h:83
Double_t GetPt() const
Definition MCTrack.h:109
Double_t GetStartVertexCoordinatesZ() const
Definition MCTrack.h:84
Double_t R2() const
production radius squared
Definition MCTrack.h:88
Double_t GetStartVertexMomentumY() const
Definition MCTrack.h:80
Double_t GetTgl() const
Definition MCTrack.h:142
Double_t GetStartVertexCoordinatesX() const
Definition MCTrack.h:82
Int_t GetPdgCode() const
Accessors.
Definition MCTrack.h:72
Int_t getMotherTrackId() const
Definition MCTrack.h:73
static TDatabasePDG * Instance()
void checkUpdates(o2::framework::ProcessingContext &pc)
static GRPGeomHelper & instance()
void setRequest(std::shared_ptr< GRPGeomRequest > req)
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:143
void setDCA(float d)
Definition V0.h:47
Static class with identifiers, bitmasks and names for ALICE detectors.
Definition DetID.h:58
static constexpr ID ITS
Definition DetID.h:63
static constexpr ID TPC
Definition DetID.h:64
ConfigParamRegistry const & options()
Definition InitContext.h:33
static constexpr int getLayer(int chipSW)
static bool isPhysicalPrimary(o2::MCTrack const &p, std::vector< o2::MCTrack > const &pcontainer)
Definition MCUtils.cxx:73
std::vector< o2::InteractionTimeRecord > & getEventRecords(bool withQED=false)
bool initFromDigitContext(std::string_view filename)
DigitizationContext const * getDigitizationContext() const
size_t getNEvents(int source) const
Get number of events.
o2::dataformats::MCEventHeader const & getMCEventHeader(int source, int event) const
retrieves the MCEventHeader for a given eventID and sourceID
size_t getNSources() const
Get number of sources.
std::vector< MCTrack > const & getTracks(int source, int event) const
variant returning all tracks for source and event at once
void extractCCDBInputs(o2::framework::ProcessingContext &pc)
void updateVDrift(float vdriftCorr, float vdrifRef, float driftTimeOffset=0)
static void requestCCDBInputs(std::vector< o2::framework::InputSpec > &inputs, std::vector< o2::framework::ConfigParamSpec > &options, const CorrectionMapsLoaderGloOpts &gloOpts)
void init(o2::framework::InitContext &ic)
static void requestCCDBInputs(std::vector< o2::framework::InputSpec > &inputs, bool laser=true, bool itstpcTgl=true)
void extractCCDBInputs(o2::framework::ProcessingContext &pc, bool laser=true, bool itstpcTgl=true)
const VDriftCorrFact & getVDriftObject() const
static std::string_view getSourceName(Source s)
bool isUpdated() const
TrackMCStudy(std::shared_ptr< DataRequest > dr, std::shared_ptr< o2::base::GRPGeomRequest > gr, GTrackID::mask_t src, const o2::tpc::CorrectionMapsLoaderGloOpts &sclOpts, bool checkSV)
void endOfStream(EndOfStreamContext &ec) final
This is invoked whenever we have an EndOfStream event.
void init(InitContext &ic) final
void finaliseCCDB(ConcreteDataMatcher &matcher, void *obj) final
void run(ProcessingContext &pc) final
void process(const o2::globaltracking::RecoContainer &recoData)
~TrackMCStudy() final=default
GLenum src
Definition glcorearb.h:1767
GLuint entry
Definition glcorearb.h:5735
GLdouble f
Definition glcorearb.h:310
GLenum const GLfloat * params
Definition glcorearb.h:272
GLfloat v0
Definition glcorearb.h:811
GLuint id
Definition glcorearb.h:650
@ ITSClusters
constexpr double LHCBunchSpacingMUS
constexpr int LHCMaxBunches
constexpr double LHCBunchSpacingNS
constexpr float TwoPI
constexpr float PI
constexpr float PIHalf
Node par(int index)
Parameters.
Defining PrimaryVertex explicitly as messageable.
Definition TFIDInfo.h:20
std::vector< ConfigParamSpec > Options
detail::Bracket< float > Bracketf_t
Definition Primitive2D.h:40
float angle2Alpha(float phi)
Definition Utils.h:203
int angle2Sector(float phi)
Definition Utils.h:183
std::tuple< float, float > rotateZInv(float xG, float yG, float snAlp, float csAlp)
Definition Utils.h:142
TrackParCovF TrackParCov
Definition Track.h:33
TrackParF TrackPar
Definition Track.h:29
std::pair< int, o2::dataformats::VtxTrackIndex > VTIndexV
o2::dataformats::VtxTrackRef V2TRef
o2::framework::DataProcessorSpec getTrackMCStudySpec(o2::dataformats::GlobalTrackID::mask_t srcTracks, o2::dataformats::GlobalTrackID::mask_t srcClus, const o2::tpc::CorrectionMapsLoaderGloOpts &sclOpts, bool checkSV)
create a processor spec
o2::dataformats::VtxTrackIndex VTIndex
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Defining DataPointCompositeObject explicitly as copiable.
GTrackID getITSContributorGID(GTrackID source) const
GlobalIDSet getSingleDetectorRefs(GTrackID gidx) const
const o2::tpc::TrackTPC & getTPCTrack(GTrackID id) const
const o2::tpc::ClusterNativeAccess & getTPCClusters() const
o2::MCCompLabel getTrackMCLabel(GTrackID id) const
GTrackID getTPCContributorGID(GTrackID source) const
const o2::track::TrackParCov & getTrackParam(GTrackID gidx) const
gsl::span< const unsigned char > clusterShMapTPC
externally set TPC clusters sharing map
void collectData(o2::framework::ProcessingContext &pc, const DataRequest &request)
const o2::track::TrackParCov & getTrackParamOut(GTrackID gidx) const
const o2::dataformats::PrimaryVertex & getPrimaryVertex(int i) const
const o2::its::TrackITS & getITSTrack(GTrackID gid) const
std::unique_ptr< o2::tpc::internal::getWorkflowTPCInput_ret > inputsTPCclusters
void getTrackTime(GTrackID gid, float &t, float &tErr) const
gsl::span< const unsigned int > occupancyMapTPC
externally set TPC clusters occupancy map
int lumiType
what estimator to used for corrections scaling: 0: no scaling, 1: CTP, 2: IDC
int lumiMode
what corrections method to use: 0: classical scaling, 1: Using of the derivative map,...
float refTimeOffset
additive time offset reference (\mus)
float refVDrift
reference vdrift for which factor was extracted
float getTimeOffset() const
float timeOffsetCorr
additive time offset correction (\mus)
float corrFact
drift velocity correction factor (multiplicative)
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
std::vector< int > row