Project
Loading...
Searching...
No Matches
TrackingStudy.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>
46#include "GPUO2InterfaceRefit.h"
47#include "GPUO2Interface.h" // Needed for propper settings in GPUParam.h
48#include "GPUParam.h"
49#include "GPUParam.inc"
50#include "GPUTPCGeometry.h"
52#include "MathUtils/fit.h"
53#include <TF1.h>
54
55namespace o2::trackstudy
56{
57
58using namespace o2::framework;
61
67
69
70class TrackingStudySpec : public Task
71{
72 public:
73 TrackingStudySpec(std::shared_ptr<DataRequest> dr, std::shared_ptr<o2::base::GRPGeomRequest> gr, GTrackID::mask_t src, bool useMC, const o2::tpc::CorrectionMapsLoaderGloOpts& sclOpts)
74 : mDataRequest(dr), mGGCCDBRequest(gr), mTracksSrc(src), mUseMC(useMC)
75 {
76 mTPCCorrMapsLoader.setLumiScaleType(sclOpts.lumiType);
77 mTPCCorrMapsLoader.setLumiScaleMode(sclOpts.lumiMode);
78 }
79 ~TrackingStudySpec() final = default;
80 void init(InitContext& ic) final;
81 void run(ProcessingContext& pc) final;
82 void endOfStream(EndOfStreamContext& ec) final;
83 void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final;
84 void process(o2::globaltracking::RecoContainer& recoData);
85
86 private:
87 void updateTimeDependentParams(ProcessingContext& pc);
88 float getDCAYCut(float pt) const;
89 float getDCAZCut(float pt) const;
90 std::shared_ptr<DataRequest> mDataRequest;
91 std::shared_ptr<o2::base::GRPGeomRequest> mGGCCDBRequest;
92 o2::tpc::VDriftHelper mTPCVDriftHelper{};
93 o2::tpc::CorrectionMapsLoader mTPCCorrMapsLoader{};
94 bool mUseMC{false};
95 std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOut;
96 std::unique_ptr<o2::utils::TreeStreamRedirector> mDBGOutVtx;
97 std::unique_ptr<o2::gpu::GPUO2InterfaceRefit> mTPCRefitter;
98 std::vector<float> mMltHistTB, mTBinClOccAft, mTBinClOccBef, mTBinClOccWgh;
99 std::unique_ptr<TF1> mOccWghFun;
100 float mITSROFrameLengthMUS = 0.f;
101 float mTPCTBinMUS = 0.f; // TPC bin in microseconds
102 float mTPCTBinMUSInv = 0.f;
103 int mMaxNeighbours = 3;
104 float mMaxVTTimeDiff = 80.; // \mus
105 float mTPCDCAYCut = 2.;
106 float mTPCDCAZCut = 2.;
107 float mMinX = 46.;
108 float mMaxEta = 0.8;
109 float mMinPt = 0.1;
110 int mNOccBinsDrift = 10;
111 int mMinTPCClusters = 60;
112 int mNTPCOccBinLength = 0;
113 int mNHBPerTF = 0;
114 float mNTPCOccBinLengthInv;
115 bool mStoreWithITSOnly = false;
116 bool mDoPairsCorr = false;
117 std::string mDCAYFormula = "0.0105 + 0.0350 / pow(x, 1.1)";
118 std::string mDCAZFormula = "0.0105 + 0.0350 / pow(x, 1.1)";
119 GTrackID::mask_t mTracksSrc{};
121 o2::steer::MCKinematicsReader mcReader; // reader of MC information
122};
123
125{
127 mTPCCorrMapsLoader.init(ic);
128 int lane = ic.services().get<const o2::framework::DeviceSpec>().inputTimesliceId;
129 int maxLanes = ic.services().get<const o2::framework::DeviceSpec>().maxInputTimeslices;
130 std::string dbgnm = maxLanes == 1 ? "trackStudy.root" : fmt::format("trackStudy_{}.root", lane);
131 mDBGOut = std::make_unique<o2::utils::TreeStreamRedirector>(dbgnm.c_str(), "recreate");
132 dbgnm = maxLanes == 1 ? "trackStudyVtx.root" : fmt::format("trackStudyVtx_{}.root", lane);
133 mDBGOutVtx = std::make_unique<o2::utils::TreeStreamRedirector>(dbgnm.c_str(), "recreate");
134 mStoreWithITSOnly = ic.options().get<bool>("with-its-only");
135 mMaxVTTimeDiff = ic.options().get<float>("max-vtx-timediff");
136 mMaxNeighbours = ic.options().get<int>("max-vtx-neighbours");
137 mTPCDCAYCut = ic.options().get<float>("max-tpc-dcay");
138 mTPCDCAZCut = ic.options().get<float>("max-tpc-dcaz");
139 mMinX = ic.options().get<float>("min-x-prop");
140 mMaxEta = ic.options().get<float>("max-eta");
141 mMinPt = ic.options().get<float>("min-pt");
142 mMinTPCClusters = ic.options().get<int>("min-tpc-clusters");
143 mDCAYFormula = ic.options().get<std::string>("dcay-vs-pt");
144 mDCAZFormula = ic.options().get<std::string>("dcaz-vs-pt");
145 mDoPairsCorr = ic.options().get<bool>("pair-correlations");
146 mNOccBinsDrift = ic.options().get<int>("noccbins");
147 if (mNOccBinsDrift < 3) {
148 mNOccBinsDrift = 3;
149 }
150 auto str = ic.options().get<std::string>("occ-weight-fun");
151 if (!str.empty()) {
152 mOccWghFun = std::make_unique<TF1>("occFun", str.c_str(), -100., 100.);
153 }
154}
155
157{
159 recoData.collectData(pc, *mDataRequest.get()); // select tracks of needed type, with minimal cuts, the real selected will be done in the vertexer
160 updateTimeDependentParams(pc); // Make sure this is called after recoData.collectData, which may load some conditions
161 if (recoData.inputsTPCclusters) {
162 mTPCRefitter = std::make_unique<o2::gpu::GPUO2InterfaceRefit>(&recoData.inputsTPCclusters->clusterIndex, &mTPCCorrMapsLoader, o2::base::Propagator::Instance()->getNominalBz(),
163 recoData.getTPCTracksClusterRefs().data(), 0, recoData.clusterShMapTPC.data(), recoData.occupancyMapTPC.data(),
164 recoData.occupancyMapTPC.size(), nullptr, o2::base::Propagator::Instance());
165 mTPCRefitter->setTrackReferenceX(900); // disable propagation after refit by setting reference to value > 500
166 mNTPCOccBinLength = mTPCRefitter->getParam()->rec.tpc.occupancyMapTimeBins;
167 mTBinClOccBef.clear();
168 mTBinClOccAft.clear();
169 mTBinClOccWgh.clear();
170 }
171
172 // prepare TPC occupancy data
173 if (mNTPCOccBinLength > 1 && recoData.occupancyMapTPC.size()) {
174 mNTPCOccBinLengthInv = 1. / mNTPCOccBinLength;
175 int nTPCBins = mNHBPerTF * o2::constants::lhc::LHCMaxBunches / 8, ninteg = 0;
176 int nTPCOccBins = nTPCBins * mNTPCOccBinLengthInv, sumBins = std::max(1, int(o2::constants::lhc::LHCMaxBunches / 8 * mNTPCOccBinLengthInv));
177 mTBinClOccAft.resize(nTPCOccBins);
178 mTBinClOccBef.resize(nTPCOccBins);
179 float sm = 0., tb = 0.5 * mNTPCOccBinLength;
180 mMltHistTB.resize(nTPCOccBins);
181 for (int i = 0; i < nTPCOccBins; i++) {
182 mMltHistTB[i] = mTPCRefitter->getParam()->GetUnscaledMult(tb);
183 tb += mNTPCOccBinLength;
184 }
185 for (int i = nTPCOccBins; i--;) {
186 sm += mMltHistTB[i];
187 if (i + sumBins < nTPCOccBins) {
188 sm -= mMltHistTB[i + sumBins];
189 }
190 mTBinClOccAft[i] = sm;
191 }
192 sm = 0;
193 for (int i = 0; i < nTPCOccBins; i++) {
194 sm += mMltHistTB[i];
195 if (i - sumBins > 0) {
196 sm -= mMltHistTB[i - sumBins];
197 }
198 mTBinClOccBef[i] = sm;
199 }
200 } else {
201 mTBinClOccBef.resize(1);
202 mTBinClOccAft.resize(1);
203 }
204
205 process(recoData);
206}
207
208void TrackingStudySpec::updateTimeDependentParams(ProcessingContext& pc)
209{
211 mTPCVDriftHelper.extractCCDBInputs(pc);
212 mTPCCorrMapsLoader.extractCCDBInputs(pc);
213 static bool initOnceDone = false;
214 if (!initOnceDone) { // this params need to be queried only once
215 initOnceDone = true;
216 // Note: reading of the ITS AlpideParam needed for ITS timing is done by the RecoContainer
219 if (!grp->isDetContinuousReadOut(DetID::ITS)) {
220 mITSROFrameLengthMUS = alpParams.roFrameLengthTrig / 1.e3; // ITS ROFrame duration in \mus
221 } else {
222 mITSROFrameLengthMUS = alpParams.roFrameLengthInBC * o2::constants::lhc::LHCBunchSpacingNS * 1e-3; // ITS ROFrame duration in \mus
223 }
225 mNHBPerTF = o2::base::GRPGeomHelper::instance().getGRPECS()->getNHBFPerTF();
227 mTPCTBinMUS = elParam.ZbinWidth; // TPC bin in microseconds
228 mTPCTBinMUSInv = 1. / mTPCTBinMUS;
229 }
230 bool updateMaps = false;
231 if (mTPCCorrMapsLoader.isUpdated()) {
232 mTPCCorrMapsLoader.acknowledgeUpdate();
233 updateMaps = true;
234 }
235 if (mTPCVDriftHelper.isUpdated()) {
236 LOGP(info, "Updating TPC fast transform map with new VDrift factor of {} wrt reference {} and DriftTimeOffset correction {} wrt {} from source {}",
237 mTPCVDriftHelper.getVDriftObject().corrFact, mTPCVDriftHelper.getVDriftObject().refVDrift,
238 mTPCVDriftHelper.getVDriftObject().timeOffsetCorr, mTPCVDriftHelper.getVDriftObject().refTimeOffset,
239 mTPCVDriftHelper.getSourceName());
240 mTPCVDriftHelper.acknowledgeUpdate();
241 updateMaps = true;
242 }
243 if (updateMaps) {
244 mTPCCorrMapsLoader.updateVDrift(mTPCVDriftHelper.getVDriftObject().corrFact, mTPCVDriftHelper.getVDriftObject().refVDrift, mTPCVDriftHelper.getVDriftObject().getTimeOffset());
245 }
246}
247
249{
250 auto pvvec = recoData.getPrimaryVertices();
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 auto FITInfo = recoData.getFT0RecPoints();
255 static int TFCount = 0;
256 int nv = vtxRefs.size();
258 o2::dataformats::PrimaryVertexExt vtxDummy(mMeanVtx.getPos(), {}, {}, 0);
259 std::vector<o2::dataformats::PrimaryVertexExt> pveVec(nv);
260 std::vector<float> tpcOccAftV, tpcOccBefV;
261 pveVec.back() = vtxDummy;
263 float tBiasITS = alpParams.roFrameBiasInBC * o2::constants::lhc::LHCBunchSpacingMUS;
265 std::vector<o2::dataformats::TrackInfoExt> trcExtVec;
266 std::vector<o2::trackstudy::TrackPairInfo> trcPairsVec;
267 auto vdrift = mTPCVDriftHelper.getVDriftObject().getVDrift();
268 float maxDriftTB = 250.f / vdrift / (o2::constants::lhc::LHCBunchSpacingMUS * 8);
269 int groupOcc = std::ceil(maxDriftTB / mNOccBinsDrift / mNTPCOccBinLength);
270
271 bool tpcTrackOK = recoData.isTrackSourceLoaded(GTrackID::TPC);
272
273 auto fillTPCClInfo = [&recoData, this](const o2::tpc::TrackTPC& trc, o2::dataformats::TrackInfoExt& trExt, float timestampTB = -1e9) {
274 const auto clRefs = recoData.getTPCTracksClusterRefs();
275 const auto tpcClusAcc = recoData.getTPCClusters();
276 const auto shMap = recoData.clusterShMapTPC;
277 if (recoData.inputsTPCclusters) {
278 uint8_t clSect = 0, clRow = 0, clRowP = -1;
279 uint32_t clIdx = 0;
280 for (int ic = 0; ic < trc.getNClusterReferences(); ic++) {
281 trc.getClusterReference(clRefs, ic, clSect, clRow, clIdx);
282 if (clRow != clRowP) {
283 trExt.rowCountTPC++;
284 clRowP = clRow;
285 }
286 unsigned int absoluteIndex = tpcClusAcc.clusterOffset[clSect][clRow] + clIdx;
287 if (shMap[absoluteIndex] & o2::gpu::GPUTPCGMMergedTrackHit::flagShared) {
288 trExt.nClTPCShared++;
289 }
290 }
291 trc.getClusterReference(clRefs, trc.getNClusterReferences() - 1, clSect, clRow, clIdx);
292 trExt.rowMinTPC = clRow;
293 const auto& clus = tpcClusAcc.clusters[clSect][clRow][clIdx];
294 trExt.padFromEdge = uint8_t(clus.getPad());
295 int npads = o2::gpu::GPUTPCGeometry::NPads(clRow);
296 if (trExt.padFromEdge > npads / 2) {
297 trExt.padFromEdge = npads - 1 - trExt.padFromEdge;
298 }
299 this->mTPCCorrMapsLoader.Transform(clSect, clRow, clus.getPad(), clus.getTime(), trExt.innerTPCPos0[0], trExt.innerTPCPos0[1], trExt.innerTPCPos0[2], trc.getTime0()); // nominal time of the track
300 if (timestampTB > -1e8) {
301 this->mTPCCorrMapsLoader.Transform(clSect, clRow, clus.getPad(), clus.getTime(), trExt.innerTPCPos[0], trExt.innerTPCPos[1], trExt.innerTPCPos[2], timestampTB); // time assigned from the global track track
302 } else {
303 trExt.innerTPCPos = trExt.innerTPCPos0;
304 }
305 trc.getClusterReference(clRefs, 0, clSect, clRow, clIdx);
306 trExt.rowMaxTPC = clRow;
307 }
308 };
309
310 auto getTPCPairSharing = [&recoData, this](const o2::tpc::TrackTPC& trc0, const o2::tpc::TrackTPC& trc1) {
311 const auto clRefs = recoData.getTPCTracksClusterRefs();
312 uint8_t nsh = 0, nshRows = 0, lastSharedRow = -1;
313 if (recoData.inputsTPCclusters) {
314 uint8_t clSect0 = 0, clRow0 = 0, clSect1 = 0, clRow1 = 0;
315 uint32_t clIdx0 = 0, clIdx1 = 0;
316 int ic1Start = 0;
317 for (int ic0 = 0; ic0 < trc0.getNClusterReferences(); ic0++) { // outside -> inside
318 trc0.getClusterReference(clRefs, ic0, clSect0, clRow0, clIdx0);
319 for (int ic1 = ic1Start; ic1 < trc1.getNClusterReferences(); ic1++) { // outside -> inside
320 trc1.getClusterReference(clRefs, ic1, clSect1, clRow1, clIdx1);
321 if (clRow1 > clRow0) {
322 ic1Start = ic1 + 1;
323 continue; // catch up ic0
324 }
325 if (clRow1 == clRow0) {
326 if (clSect0 == clSect1 && clIdx0 == clIdx1) {
327 nsh++;
328 if (lastSharedRow != clRow0) {
329 lastSharedRow = clRow0;
330 nshRows++;
331 }
332 ic1Start = ic1 + 1;
333 break; // check next ic0
334 }
335 }
336 }
337 }
338 }
339 return std::make_pair(nsh, nshRows);
340 };
341
342 auto assignRecTrack = [&recoData, this](const o2::dataformats::TrackInfoExt& src, o2::trackstudy::RecTrack& dst) {
343 dst.track = src.track;
344 dst.gid = src.gid;
345 dst.ts.setTimeStamp(src.ttime);
346 dst.ts.setTimeStampError(src.ttimeE);
347 dst.nClITS = src.nClITS;
348 dst.nClTPC = src.nClTPC;
349 dst.pattITS = src.pattITS;
350 if (src.q2ptITS == 0. && dst.nClITS > 0) {
351 dst.pattITS |= 0x1 << 7;
352 }
353 dst.lowestPadRow = src.rowMinTPC;
354 if (this->mUseMC) {
355 auto gidSet = recoData.getSingleDetectorRefs(src.gid);
356 if (recoData.getTrackMCLabel(src.gid).isFake()) {
357 dst.flags |= RecTrack::FakeGLO;
358 }
359 auto msk = src.gid.getSourceDetectorsMask();
360 if (msk[DetID::ITS]) {
361 if (gidSet[GTrackID::ITS].isSourceSet()) { // has ITS track rather than AB tracklet
362 auto lblITS = recoData.getTrackMCLabel(gidSet[GTrackID::ITS]);
363 if (lblITS.isFake()) {
364 dst.flags |= RecTrack::FakeITS;
365 }
366 } else { // AB ITS tracklet
367 if (recoData.getTrackMCLabel(gidSet[GTrackID::ITSAB]).isFake()) {
368 dst.flags |= RecTrack::FakeITS;
369 }
370 }
371 if (msk[DetID::TPC]) { // has both ITS and TPC contribution
372 if (recoData.getTrackMCLabel(gidSet[GTrackID::ITSTPC]).isFake()) {
373 dst.flags |= RecTrack::FakeITSTPC;
374 }
375 }
376 }
377 if (msk[DetID::TPC]) {
378 if (recoData.getTrackMCLabel(gidSet[GTrackID::TPC]).isFake()) {
379 dst.flags |= RecTrack::FakeTPC;
380 }
381 }
382 }
383 };
384 tpcOccAftV.resize(mNOccBinsDrift);
385 tpcOccBefV.resize(mNOccBinsDrift);
386
387 for (int iv = 0; iv < nv; iv++) {
388 LOGP(debug, "processing PV {} of {}", iv, nv);
389 const auto& vtref = vtxRefs[iv];
390 if (iv != nv - 1) {
391 auto& pve = pveVec[iv];
392 static_cast<o2::dataformats::PrimaryVertex&>(pve) = pvvec[iv];
393 // find best matching FT0 signal
394 float bestTimeDiff = 1000, bestTime = -999;
395 int bestFTID = -1;
396 if (mTracksSrc[GTrackID::FT0]) {
397 for (int ift0 = vtref.getFirstEntryOfSource(GTrackID::FT0); ift0 < vtref.getFirstEntryOfSource(GTrackID::FT0) + vtref.getEntriesOfSource(GTrackID::FT0); ift0++) {
398 const auto& ft0 = FITInfo[trackIndex[ift0]];
399 if (ft0Params.isSelected(ft0)) {
400 auto fitTime = ft0.getInteractionRecord().differenceInBCMUS(recoData.startIR);
401 if (std::abs(fitTime - pve.getTimeStamp().getTimeStamp()) < bestTimeDiff) {
402 bestTimeDiff = fitTime - pve.getTimeStamp().getTimeStamp();
403 bestFTID = trackIndex[ift0];
404 }
405 }
406 }
407 } else {
408 LOGP(warn, "FT0 is not requested, cannot set complete vertex info");
409 }
410 if (bestFTID >= 0) {
411 pve.FT0A = FITInfo[bestFTID].getTrigger().getAmplA();
412 pve.FT0C = FITInfo[bestFTID].getTrigger().getAmplC();
413 pve.FT0Time = double(FITInfo[bestFTID].getInteractionRecord().differenceInBCMUS(recoData.startIR)) + FITInfo[bestFTID].getCollisionTimeMean() * 1e-6; // time in \mus
414 }
415 pve.VtxID = iv;
416 }
417 trcExtVec.clear();
418 trcPairsVec.clear();
419 float q2ptITS, q2ptTPC, q2ptITSTPC, q2ptITSTPCTRD;
420 for (int is = 0; is < GTrackID::NSources; is++) {
421 DetID::mask_t dm = GTrackID::getSourceDetectorsMask(is);
422 bool skipTracks = !mTracksSrc[is] || !recoData.isTrackSourceLoaded(is) || !(dm[DetID::ITS] || dm[DetID::TPC]);
423 int idMin = vtref.getFirstEntryOfSource(is), idMax = idMin + vtref.getEntriesOfSource(is);
424 for (int i = idMin; i < idMax; i++) {
425 auto vid = trackIndex[i];
426 bool pvCont = vid.isPVContributor();
427 if (pvCont) {
428 pveVec[iv].nSrc[is]++;
429 }
430 if (skipTracks) {
431 continue;
432 }
433 GTrackID tpcTrID;
434 const o2::tpc::TrackTPC* tpcTr = nullptr;
435 int nclTPC = 0;
436 if (dm[DetID::TPC] && tpcTrackOK) {
437 tpcTrID = recoData.getTPCContributorGID(vid);
438 tpcTr = &recoData.getTPCTrack(tpcTrID);
439 nclTPC = tpcTr->getNClusters();
440 if (nclTPC < mMinTPCClusters) {
441 continue;
442 }
443 }
444 bool ambig = vid.isAmbiguous();
445 auto trc = recoData.getTrackParam(vid);
446 if (abs(trc.getEta()) > mMaxEta) {
447 continue;
448 }
449 if (iv < nv - 1 && is == GTrackID::TPC && tpcTr && !tpcTr->hasBothSidesClusters()) { // for unconstrained TPC tracks correct track Z
450 float corz = vdrift * (tpcTr->getTime0() * mTPCTBinMUS - pvvec[iv].getTimeStamp().getTimeStamp());
451 if (tpcTr->hasASideClustersOnly()) {
452 corz = -corz; // A-side
453 }
454 trc.setZ(trc.getZ() + corz);
455 }
456 float xmin = trc.getX();
458 if (!prop->propagateToDCA(iv == nv - 1 ? vtxDummy : pvvec[iv], trc, prop->getNominalBz(), 2., o2::base::PropagatorF::MatCorrType::USEMatCorrLUT, &dca)) {
459 continue;
460 }
461 bool hasITS = GTrackID::getSourceDetectorsMask(is)[GTrackID::ITS];
462 if (std::abs(dca.getY()) > (hasITS ? getDCAYCut(trc.getPt()) : mTPCDCAYCut) ||
463 std::abs(dca.getZ()) > (hasITS ? getDCAZCut(trc.getPt()) : mTPCDCAZCut)) {
464 continue;
465 }
466 if (trc.getPt() < mMinPt) {
467 continue;
468 }
469 if (iv != nv - 1) {
470 pveVec[iv].nSrcA[is]++;
471 if (ambig) {
472 pveVec[iv].nSrcAU[is]++;
473 }
474 }
475 if (!hasITS && mStoreWithITSOnly) {
476 continue;
477 }
478 {
479 auto& trcExt = trcExtVec.emplace_back();
480 recoData.getTrackTime(vid, trcExt.ttime, trcExt.ttimeE);
481 trcExt.track = trc;
482 trcExt.dca = dca;
483 trcExt.gid = vid;
484 trcExt.xmin = xmin;
485 trcExt.dcaTPC.set(-999.f, -999.f);
486
487 if (tpcTr) {
488 float tsuse = trcExt.ttime / (8 * o2::constants::lhc::LHCBunchSpacingMUS);
489 if (tpcTr->hasASideClusters()) {
490 trcExt.setTPCA();
491 }
492 if (tpcTr->hasCSideClusters()) {
493 trcExt.setTPCC();
494 }
495 if (is == GTrackID::TPC) {
496 trcExt.dcaTPC = dca;
497 tsuse = -1e9;
498 } else {
499 o2::track::TrackParCov tmpTPC(*tpcTr);
500 if (iv < nv - 1 && is == GTrackID::TPC && tpcTr && !tpcTr->hasBothSidesClusters()) { // for unconstrained TPC tracks correct track Z
501 float corz = vdrift * (tpcTr->getTime0() * mTPCTBinMUS - pvvec[iv].getTimeStamp().getTimeStamp());
502 if (tpcTr->hasASideClustersOnly()) {
503 corz = -corz; // A-side
504 }
505 tmpTPC.setZ(tmpTPC.getZ() + corz);
506 }
507 if (!prop->propagateToDCA(iv == nv - 1 ? vtxDummy : pvvec[iv], tmpTPC, prop->getNominalBz(), 2., o2::base::PropagatorF::MatCorrType::USEMatCorrLUT, &trcExt.dcaTPC)) {
508 trcExt.dcaTPC.set(-999.f, -999.f);
509 }
510 }
511 fillTPCClInfo(*tpcTr, trcExt, tsuse);
512 trcExt.chi2TPC = tpcTr->getChi2();
513 }
514 auto gidRefs = recoData.getSingleDetectorRefs(vid);
515 if (gidRefs[GTrackID::ITS].isIndexSet()) {
516 const auto& itsTr = recoData.getITSTrack(gidRefs[GTrackID::ITS]);
517 trcExt.q2ptITS = itsTr.getQ2Pt();
518 trcExt.nClITS = itsTr.getNClusters();
519 for (int il = 0; il < 7; il++) {
520 if (itsTr.hasHitOnLayer(il)) {
521 trcExt.pattITS |= 0x1 << il;
522 }
523 }
524 } else if (gidRefs[GTrackID::ITSAB].isIndexSet()) {
525 const auto& itsTrf = recoData.getITSABRefs()[gidRefs[GTrackID::ITSAB]];
526 trcExt.nClITS = itsTrf.getNClusters();
527 for (int il = 0; il < 7; il++) {
528 if (itsTrf.hasHitOnLayer(il)) {
529 trcExt.pattITS |= 0x1 << il;
530 }
531 }
532 }
533 if (gidRefs[GTrackID::TPC].isIndexSet()) {
534 trcExt.q2ptTPC = recoData.getTrackParam(gidRefs[GTrackID::TPC]).getQ2Pt();
535 trcExt.nClTPC = nclTPC;
536 }
537 if (gidRefs[GTrackID::ITSTPC].isIndexSet()) {
538 const auto& trTPCITS = recoData.getTPCITSTrack(gidRefs[GTrackID::ITSTPC]);
539 trcExt.q2ptITSTPC = trTPCITS.getQ2Pt();
540 trcExt.chi2ITSTPC = trTPCITS.getChi2Match();
541 }
542 if (gidRefs[GTrackID::TRD].isIndexSet()) {
543 trcExt.q2ptITSTPCTRD = recoData.getTrackParam(gidRefs[GTrackID::TRD]).getQ2Pt();
544 }
545 if (gidRefs[GTrackID::TOF].isIndexSet()) {
546 trcExt.infoTOF = recoData.getTOFMatch(vid);
547 }
548 }
549 }
550 }
551 float tpcOccBef = 0., tpcOccAft = 0.;
552 if (iv != nv - 1) {
553 int tb = pveVec[iv].getTimeStamp().getTimeStamp() * mTPCTBinMUSInv * mNTPCOccBinLengthInv;
554 tpcOccBef = tb < 0 ? mTBinClOccBef[0] : (tb >= mTBinClOccBef.size() ? mTBinClOccBef.back() : mTBinClOccBef[tb]);
555 tpcOccAft = tb < 0 ? mTBinClOccAft[0] : (tb >= mTBinClOccAft.size() ? mTBinClOccAft.back() : mTBinClOccAft[tb]);
556 int tbc = pveVec[iv].getTimeStamp().getTimeStamp() * mTPCTBinMUSInv * mNTPCOccBinLengthInv - groupOcc / 2.;
557 for (int iob = 0; iob < mNOccBinsDrift; iob++) {
558 float sm = 0;
559 for (int ig = 0; ig < groupOcc; ig++) {
560 int ocb = tbc + ig + groupOcc * iob;
561 if (ocb < 0 || ocb >= (int)mMltHistTB.size()) {
562 sm = -1;
563 break;
564 }
565 sm += mMltHistTB[ocb];
566 }
567 tpcOccAftV[iob] = sm;
568 //
569 sm = 0;
570 for (int ig = 0; ig < groupOcc; ig++) {
571 int ocb = tbc + ig - groupOcc * iob;
572 if (ocb < 0 || ocb >= (int)mMltHistTB.size()) {
573 sm = -1;
574 break;
575 }
576 sm += mMltHistTB[ocb];
577 }
578 tpcOccBefV[iob] = sm;
579 }
580 }
581 (*mDBGOut) << "trpv"
582 << "orbit=" << recoData.startIR.orbit << "tfID=" << TFCount
583 << "tpcOccBef=" << tpcOccBef << "tpcOccAft=" << tpcOccAft
584 << "tpcOccBefV=" << tpcOccBefV << "tpcOccAftV=" << tpcOccAftV
585 << "pve=" << pveVec[iv] << "trc=" << trcExtVec << "\n";
586
587 if (mDoPairsCorr) {
588 for (int it0 = 0; it0 < (int)trcExtVec.size(); it0++) {
589 const auto& tr0 = trcExtVec[it0];
590 if (tr0.nClTPC < 1) {
591 continue;
592 }
593 for (int it1 = it0 + 1; it1 < (int)trcExtVec.size(); it1++) {
594 const auto& tr1 = trcExtVec[it1];
595 if (tr1.nClTPC < 1) {
596 continue;
597 }
598
599 if (std::abs(tr0.track.getTgl() - tr1.track.getTgl()) > 0.25) {
600 continue;
601 }
602 auto dphi = tr0.track.getPhi() - tr1.track.getPhi();
603 if (dphi < -o2::constants::math::PI) {
605 } else if (dphi > o2::constants::math::PI) {
607 }
608 if (std::abs(dphi) > 0.25) {
609 continue;
610 }
611 auto& pr = trcPairsVec.emplace_back();
612 assignRecTrack(tr0, pr.tr0);
613 assignRecTrack(tr1, pr.tr1);
614 auto shinfo = getTPCPairSharing(recoData.getTPCTrack(recoData.getTPCContributorGID(tr0.gid)), recoData.getTPCTrack(recoData.getTPCContributorGID(tr1.gid)));
615 pr.nshTPC = shinfo.first;
616 pr.nshTPCRow = shinfo.second;
617 }
618 }
619 (*mDBGOut) << "pairs" << "pr=" << trcPairsVec << "\n";
620 }
621 }
622
623 int nvtot = mMaxNeighbours < 0 ? -1 : (int)pveVec.size();
624
625 auto insSlot = [maxSlots = mMaxNeighbours](std::vector<float>& vc, float v, int slot, std::vector<int>& vid, int id) {
626 for (int i = maxSlots - 1; i > slot; i--) {
627 std::swap(vc[i], vc[i - 1]);
628 std::swap(vid[i], vid[i - 1]);
629 }
630 vc[slot] = v;
631 vid[slot] = id;
632 };
633
634 for (int cnt = 0; cnt < nvtot; cnt++) {
635 const auto& pve = pveVec[cnt];
636 float tv = pve.getTimeStamp().getTimeStamp();
637 std::vector<o2::dataformats::PrimaryVertexExt> pveT(mMaxNeighbours); // neighbours in time
638 std::vector<o2::dataformats::PrimaryVertexExt> pveZ(mMaxNeighbours); // neighbours in Z
639 std::vector<int> idT(mMaxNeighbours), idZ(mMaxNeighbours);
640 std::vector<float> dT(mMaxNeighbours), dZ(mMaxNeighbours);
641 for (int i = 0; i < mMaxNeighbours; i++) {
642 idT[i] = idZ[i] = -1;
643 dT[i] = mMaxVTTimeDiff;
644 dZ[i] = 1e9;
645 }
646 int cntM = cnt - 1, cntP = cnt + 1;
647 for (; cntM >= 0; cntM--) { // backward
648 const auto& vt = pveVec[cntM];
649 auto dtime = std::abs(tv - vt.getTimeStamp().getTimeStamp());
650 if (dtime > mMaxVTTimeDiff) {
651 continue;
652 }
653 for (int i = 0; i < mMaxNeighbours; i++) {
654 if (dT[i] > dtime) {
655 insSlot(dT, dtime, i, idT, cntM);
656 break;
657 }
658 }
659 auto dz = std::abs(pve.getZ() - vt.getZ());
660 for (int i = 0; i < mMaxNeighbours; i++) {
661 if (dZ[i] > dz) {
662 insSlot(dZ, dz, i, idZ, cntM);
663 break;
664 }
665 }
666 }
667 for (; cntP < nvtot; cntP++) { // forward
668 const auto& vt = pveVec[cntP];
669 auto dtime = std::abs(tv - vt.getTimeStamp().getTimeStamp());
670 if (dtime > mMaxVTTimeDiff) {
671 continue;
672 }
673 for (int i = 0; i < mMaxNeighbours; i++) {
674 if (dT[i] > dtime) {
675 insSlot(dT, dtime, i, idT, cntP);
676 break;
677 }
678 }
679 auto dz = std::abs(pve.getZ() - vt.getZ());
680 for (int i = 0; i < mMaxNeighbours; i++) {
681 if (dZ[i] > dz) {
682 insSlot(dZ, dz, i, idZ, cntP);
683 break;
684 }
685 }
686 }
687 for (int i = 0; i < mMaxNeighbours; i++) {
688 if (idT[i] != -1) {
689 pveT[i] = pveVec[idT[i]];
690 } else {
691 break;
692 }
693 }
694 for (int i = 0; i < mMaxNeighbours; i++) {
695 if (idZ[i] != -1) {
696 pveZ[i] = pveVec[idZ[i]];
697 } else {
698 break;
699 }
700 }
701 (*mDBGOutVtx) << "pvExt"
702 << "pve=" << pve
703 << "pveT=" << pveT
704 << "pveZ=" << pveZ
705 << "tfID=" << TFCount
706 << "\n";
707 }
708
709 TFCount++;
710}
711
713{
714 mDBGOut.reset();
715 mDBGOutVtx.reset();
716}
717
719{
721 return;
722 }
723 if (mTPCVDriftHelper.accountCCDBInputs(matcher, obj)) {
724 return;
725 }
726 if (mTPCCorrMapsLoader.accountCCDBInputs(matcher, obj)) {
727 return;
728 }
729 if (matcher == ConcreteDataMatcher("GLO", "MEANVERTEX", 0)) {
730 LOG(info) << "Imposing new MeanVertex: " << ((const o2::dataformats::MeanVertexObject*)obj)->asString();
731 mMeanVtx = *(const o2::dataformats::MeanVertexObject*)obj;
732 return;
733 }
734}
735
736float TrackingStudySpec::getDCAYCut(float pt) const
737{
738 static TF1 fun("dcayvspt", mDCAYFormula.c_str(), 0, 20);
739 return fun.Eval(pt);
740}
741
742float TrackingStudySpec::getDCAZCut(float pt) const
743{
744 static TF1 fun("dcazvspt", mDCAZFormula.c_str(), 0, 20);
745 return fun.Eval(pt);
746}
747
749{
750 std::vector<OutputSpec> outputs;
751 auto dataRequest = std::make_shared<DataRequest>();
752
753 dataRequest->requestTracks(srcTracks, useMC);
754 dataRequest->requestClusters(srcClusters, useMC);
755 dataRequest->requestPrimaryVertices(useMC);
756 dataRequest->inputs.emplace_back("meanvtx", "GLO", "MEANVERTEX", 0, Lifetime::Condition, ccdbParamSpec("GLO/Calib/MeanVertex", {}, 1));
757 auto ggRequest = std::make_shared<o2::base::GRPGeomRequest>(false, // orbitResetTime
758 true, // GRPECS=true
759 true, // GRPLHCIF
760 true, // GRPMagField
761 true, // askMatLUT
763 dataRequest->inputs,
764 true);
765
766 Options opts{
767 {"max-vtx-neighbours", VariantType::Int, 3, {"Max PV neighbours fill, no PV study if < 0"}},
768 {"max-vtx-timediff", VariantType::Float, 90.f, {"Max PV time difference to consider"}},
769 {"dcay-vs-pt", VariantType::String, "0.0105 + 0.0350 / pow(x, 1.1)", {"Formula for global tracks DCAy vs pT cut"}},
770 {"dcaz-vs-pt", VariantType::String, "0.0105 + 0.0350 / pow(x, 1.1)", {"Formula for global tracks DCAy vs pT cut"}},
771 {"min-tpc-clusters", VariantType::Int, 60, {"Cut on TPC clusters"}},
772 {"max-tpc-dcay", VariantType::Float, 5.f, {"Cut on TPC dcaY"}},
773 {"max-tpc-dcaz", VariantType::Float, 5.f, {"Cut on TPC dcaZ"}},
774 {"max-eta", VariantType::Float, 1.0f, {"Cut on track eta"}},
775 {"min-pt", VariantType::Float, 0.1f, {"Cut on track pT"}},
776 {"with-its-only", VariantType::Bool, false, {"Store tracks with ITS only"}},
777 {"pair-correlations", VariantType::Bool, false, {"Do pairs correlation"}},
778 {"occ-weight-fun", VariantType::String, "(x>=-40&&x<-5) ? (1./1225*pow(x+40,2)) : ((x>-5&&x<15) ? 1. : ((x>=15&&x<40) ? (-0.4/25*x+1.24 ) : ( (x>40&&x<100) ? -0.4/60*x+0.6+0.8/3 : 0)))", {"Occupancy weighting f-n vs time in musec"}},
779 {"noccbins", VariantType::Int, 10, {"Number of occupancy bins per full drift time"}},
780 {"min-x-prop", VariantType::Float, 100.f, {"track should be propagated to this X at least"}},
781 };
782 o2::tpc::VDriftHelper::requestCCDBInputs(dataRequest->inputs);
783 o2::tpc::CorrectionMapsLoader::requestCCDBInputs(dataRequest->inputs, opts, sclOpts);
784
785 return DataProcessorSpec{
786 "track-study",
787 dataRequest->inputs,
788 outputs,
789 AlgorithmSpec{adaptFromTask<TrackingStudySpec>(dataRequest, ggRequest, srcTracks, useMC, sclOpts)},
790 opts};
791}
792
793} // namespace o2::trackstudy
Helper class to access load maps from CCDB.
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.
Definition of the Names Generator class.
Definition of the parameter class for the detector electronics.
Wrapper container for different reconstructed object types.
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 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
math_utils::Point3D< float > & getPos()
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
ServiceRegistryRef services()
Definition InitContext.h:34
ConfigParamRegistry const & options()
Definition InitContext.h:33
decltype(auto) get(R binding, int part=0) const
InputRecord & inputs()
The inputs associated with this processing context.
void extractCCDBInputs(o2::framework::ProcessingContext &pc)
void updateVDrift(float vdriftCorr, float vdrifRef, float driftTimeOffset=0)
bool accountCCDBInputs(const o2::framework::ConcreteDataMatcher &matcher, void *obj)
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
bool accountCCDBInputs(const o2::framework::ConcreteDataMatcher &matcher, void *obj)
static std::string_view getSourceName(Source s)
bool isUpdated() const
void process(o2::globaltracking::RecoContainer &recoData)
void run(ProcessingContext &pc) final
TrackingStudySpec(std::shared_ptr< DataRequest > dr, std::shared_ptr< o2::base::GRPGeomRequest > gr, GTrackID::mask_t src, bool useMC, const o2::tpc::CorrectionMapsLoaderGloOpts &sclOpts)
void finaliseCCDB(ConcreteDataMatcher &matcher, void *obj) final
void endOfStream(EndOfStreamContext &ec) final
This is invoked whenever we have an EndOfStream event.
void init(InitContext &ic) final
GLenum src
Definition glcorearb.h:1767
const GLdouble * v
Definition glcorearb.h:832
GLenum GLenum dst
Definition glcorearb.h:1767
GLuint id
Definition glcorearb.h:650
constexpr double LHCBunchSpacingMUS
constexpr int LHCMaxBunches
constexpr double LHCBunchSpacingNS
constexpr float TwoPI
constexpr float PI
Defining PrimaryVertex explicitly as messageable.
Definition TFIDInfo.h:20
std::vector< ConfigParamSpec > ccdbParamSpec(std::string const &path, int runDependent, std::vector< CCDBMetadata > metadata={}, int qrate=0)
std::vector< ConfigParamSpec > Options
detail::Bracket< float > Bracketf_t
Definition Primitive2D.h:40
TrackParCovF TrackParCov
Definition Track.h:33
o2::math_utils::Bracketf_t TBracket
o2::framework::DataProcessorSpec getTrackingStudySpec(o2::dataformats::GlobalTrackID::mask_t srcTracks, o2::dataformats::GlobalTrackID::mask_t srcClus, bool useMC, const o2::tpc::CorrectionMapsLoaderGloOpts &sclOpts)
create a processor spec
o2::dataformats::VtxTrackRef V2TRef
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.
uint32_t orbit
LHC orbit.
bool isSelected(const RecPoints &rp) 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
const o2::dataformats::TrackTPCITS & getTPCITSTrack(GTrackID gid) const
void collectData(o2::framework::ProcessingContext &pc, const DataRequest &request)
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
const o2::dataformats::MatchInfoTOF & getTOFMatch(GTrackID id) const
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"
const std::string str