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