Skip to content

Commit 73c6a33

Browse files
romainschotterromainschotter
committed
Consider propagation of Tracks + add method to calculate DCA without propagation
1 parent 7660cad commit 73c6a33

File tree

1 file changed

+128
-44
lines changed

1 file changed

+128
-44
lines changed

‎Common/Tools/TrackPropagationModule.h‎

Lines changed: 128 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -272,58 +272,81 @@ class TrackPropagationModule
272272
// std::array<float, 3> trackPxPyPzTuned = {0.0, 0.0, 0.0};
273273
double q2OverPtNew = -9999.;
274274
// Only propagate tracks which have passed the innermost wall of the TPC (e.g. skipping loopers etc). Others fill unpropagated.
275-
if (track.trackType() == o2::aod::track::TrackIU && track.x() < cGroup.minPropagationRadius.value) {
276-
if (fillTracksCov) {
277-
if constexpr (isMc) { // checking MC and fillCovMat block begins
278-
// bool hasMcParticle = track.has_mcParticle();
279-
if (cGroup.useTrackTuner.value) {
280-
trackTunedTracks->Fill(1); // all tracks
281-
bool hasMcParticle = track.has_mcParticle();
282-
if (hasMcParticle) {
283-
auto mcParticle = track.mcParticle();
284-
trackTunerObj.tuneTrackParams(mcParticle, mTrackParCov, matCorr, &mDcaInfoCov, trackTunedTracks);
285-
q2OverPtNew = mTrackParCov.getQ2Pt();
275+
if (track.x() < cGroup.minPropagationRadius.value) {
276+
if (track.trackType() == o2::aod::track::TrackIU) {
277+
if (fillTracksCov) {
278+
if constexpr (isMc) { // checking MC and fillCovMat block begins
279+
// bool hasMcParticle = track.has_mcParticle();
280+
if (cGroup.useTrackTuner.value) {
281+
trackTunedTracks->Fill(1); // all tracks
282+
bool hasMcParticle = track.has_mcParticle();
283+
if (hasMcParticle) {
284+
auto mcParticle = track.mcParticle();
285+
trackTunerObj.tuneTrackParams(mcParticle, mTrackParCov, matCorr, &mDcaInfoCov, trackTunedTracks);
286+
q2OverPtNew = mTrackParCov.getQ2Pt();
287+
}
286288
}
289+
} // MC and fillCovMat block ends
290+
}
291+
bool isPropagationOK = true;
292+
293+
if (track.has_collision()) {
294+
auto const& collision = collisions.rawIteratorAt(track.collisionId());
295+
if (fillTracksCov) {
296+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
297+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
298+
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov);
299+
} else {
300+
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo);
287301
}
288-
} // MC and fillCovMat block ends
289-
}
290-
bool isPropagationOK = true;
291-
292-
if (track.has_collision()) {
293-
auto const& collision = collisions.rawIteratorAt(track.collisionId());
294-
if (fillTracksCov) {
295-
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
296-
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
297-
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov);
298302
} else {
299-
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo);
303+
if (fillTracksCov) {
304+
mVtx.setPos({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()});
305+
mVtx.setCov(ccdbLoader.mMeanVtx->getSigmaX() * ccdbLoader.mMeanVtx->getSigmaX(), 0.0f, ccdbLoader.mMeanVtx->getSigmaY() * ccdbLoader.mMeanVtx->getSigmaY(), 0.0f, 0.0f, ccdbLoader.mMeanVtx->getSigmaZ() * ccdbLoader.mMeanVtx->getSigmaZ());
306+
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov);
307+
} else {
308+
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo);
309+
}
300310
}
301-
} else {
311+
if (isPropagationOK) {
312+
trackType = o2::aod::track::Track;
313+
}
314+
// filling some QA histograms for track tuner test purpose
302315
if (fillTracksCov) {
303-
mVtx.setPos({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()});
304-
mVtx.setCov(ccdbLoader.mMeanVtx->getSigmaX() * ccdbLoader.mMeanVtx->getSigmaX(), 0.0f, ccdbLoader.mMeanVtx->getSigmaY() * ccdbLoader.mMeanVtx->getSigmaY(), 0.0f, 0.0f, ccdbLoader.mMeanVtx->getSigmaZ() * ccdbLoader.mMeanVtx->getSigmaZ());
305-
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov);
306-
} else {
307-
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo);
316+
if constexpr (isMc) { // checking MC and fillCovMat block begins
317+
if (track.has_mcParticle() && isPropagationOK) {
318+
auto mcParticle1 = track.mcParticle();
319+
// && abs(mcParticle1.pdgCode())==211
320+
if (mcParticle1.isPhysicalPrimary()) {
321+
registry.fill(HIST("hDCAxyVsPtRec"), mDcaInfoCov.getY(), mTrackParCov.getPt());
322+
registry.fill(HIST("hDCAxyVsPtMC"), mDcaInfoCov.getY(), mcParticle1.pt());
323+
registry.fill(HIST("hDCAzVsPtRec"), mDcaInfoCov.getZ(), mTrackParCov.getPt());
324+
registry.fill(HIST("hDCAzVsPtMC"), mDcaInfoCov.getZ(), mcParticle1.pt());
325+
}
326+
}
327+
} // MC and fillCovMat block ends
308328
}
309-
}
310-
if (isPropagationOK) {
311-
trackType = o2::aod::track::Track;
312-
}
313-
// filling some QA histograms for track tuner test purpose
314-
if (fillTracksCov) {
315-
if constexpr (isMc) { // checking MC and fillCovMat block begins
316-
if (track.has_mcParticle() && isPropagationOK) {
317-
auto mcParticle1 = track.mcParticle();
318-
// && abs(mcParticle1.pdgCode())==211
319-
if (mcParticle1.isPhysicalPrimary()) {
320-
registry.fill(HIST("hDCAxyVsPtRec"), mDcaInfoCov.getY(), mTrackParCov.getPt());
321-
registry.fill(HIST("hDCAxyVsPtMC"), mDcaInfoCov.getY(), mcParticle1.pt());
322-
registry.fill(HIST("hDCAzVsPtRec"), mDcaInfoCov.getZ(), mTrackParCov.getPt());
323-
registry.fill(HIST("hDCAzVsPtMC"), mDcaInfoCov.getZ(), mcParticle1.pt());
329+
} else {
330+
if (fillTracksDCA || fillTracksDCACov) {
331+
if (track.has_collision()) {
332+
auto const& collision = collisions.rawIteratorAt(track.collisionId());
333+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
334+
if (fillTracksDCACov) {
335+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
336+
}
337+
} else {
338+
mVtx.setPos({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()});
339+
if (fillTracksDCACov) {
340+
mVtx.setCov(ccdbLoader.mMeanVtx->getSigmaX() * ccdbLoader.mMeanVtx->getSigmaX(), 0.0f, ccdbLoader.mMeanVtx->getSigmaY() * ccdbLoader.mMeanVtx->getSigmaY(), 0.0f, 0.0f, ccdbLoader.mMeanVtx->getSigmaZ() * ccdbLoader.mMeanVtx->getSigmaZ());
324341
}
325342
}
326-
} // MC and fillCovMat block ends
343+
344+
if (fillTracksDCACov) {
345+
calculateDCA(mTrackParCov, mVtx, o2::base::Propagator::Instance()->getNominalBz(), &mDcaInfoCov, 999.f);
346+
} else {
347+
calculateDCA(mTrackPar, mVtx, o2::base::Propagator::Instance()->getNominalBz(), &mDcaInfo, 999.f);
348+
}
349+
}
327350
}
328351
}
329352
// Filling modified Q/Pt values at IU/production point by track tuner in track tuner table
@@ -356,6 +379,67 @@ class TrackPropagationModule
356379
}
357380
}
358381
}
382+
383+
//_______________________________________________________________________
384+
// TTrackPar type: either aod::TrackPar or aod::TrackParCov
385+
// TVertex type: either math_utils::Point3D<value_t> or o2::dataformats::VertexBase
386+
// TDCA type: either dim2_t or o2::dataformats::DCA
387+
template<typename TTrackPar, typename TVertex, typename TDCA>
388+
bool calculateDCA(TTrackPar trackPar, const TVertex& vtx, double b, TDCA* dca, double maxD)
389+
{
390+
// propagate track to DCA to the vertex
391+
double sn, cs, alp = trackPar.getAlpha();
392+
o2::math_utils::detail::sincos(alp, sn, cs);
393+
double x = trackPar.getX(), y = trackPar.getY(), snp = trackPar.getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp));
394+
double xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
395+
x -= xv;
396+
y -= yv;
397+
// Estimate the impact parameter neglecting the track curvature
398+
double d = gpu::CAMath::Abs(x * snp - y * csp);
399+
if (d > maxD) {
400+
// provide default DCA for failed propag
401+
if constexpr (requires { trackPar.getSigmaY2(); vtx.getSigmaX2(); }) {
402+
dca->set(o2::track::DefaultDCA, o2::track::DefaultDCA,
403+
o2::track::DefaultDCACov, o2::track::DefaultDCACov, o2::track::DefaultDCACov);
404+
} else {
405+
(*dca)[0] = o2::track::DefaultDCA;
406+
(*dca)[1] = o2::track::DefaultDCA;
407+
}
408+
return false;
409+
}
410+
double crv = trackPar.getCurvature(b);
411+
double tgfv = -(crv * x - snp) / (crv * y + csp);
412+
sn = tgfv / gpu::CAMath::Sqrt(1.f + tgfv * tgfv);
413+
cs = gpu::CAMath::Sqrt((1.f - sn) * (1.f + sn));
414+
cs = (gpu::CAMath::Abs(tgfv) > constants::math::Almost0) ? sn / tgfv : constants::math::Almost1;
415+
416+
x = xv * cs + yv * sn;
417+
yv = -xv * sn + yv * cs;
418+
xv = x;
419+
420+
alp += gpu::CAMath::ASin(sn);
421+
if (!trackPar.rotate(alp)) {
422+
// provide default DCA for failed propag
423+
if constexpr (requires { trackPar.getSigmaY2(); vtx.getSigmaX2(); }) {
424+
dca->set(o2::track::DefaultDCA, o2::track::DefaultDCA,
425+
o2::track::DefaultDCACov, o2::track::DefaultDCACov, o2::track::DefaultDCACov);
426+
} else {
427+
(*dca)[0] = o2::track::DefaultDCA;
428+
(*dca)[1] = o2::track::DefaultDCA;
429+
}
430+
431+
return false;
432+
}
433+
if constexpr (requires { trackPar.getSigmaY2(); vtx.getSigmaX2(); }) {
434+
o2::math_utils::detail::sincos(alp, sn, cs);
435+
auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
436+
dca->set(trackPar.getY() - yv, trackPar.getZ() - zv, trackPar.getSigmaY2() + s2ylocvtx, trackPar.getSigmaZY(), trackPar.getSigmaZ2() + vtx.getSigmaZ2());
437+
} else {
438+
(*dca)[0] = trackPar.getY() - yv;
439+
(*dca)[1] = trackPar.getZ() - zv;
440+
}
441+
return true;
442+
}
359443
};
360444

361445
} // namespace common

0 commit comments

Comments
 (0)