Skip to content
Draft
122 changes: 103 additions & 19 deletions Common/Tools/TrackPropagationModule.h
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@

struct TrackPropagationConfigurables : o2::framework::ConfigurableGroup {
std::string prefix = "trackPropagation";
o2::framework::Configurable<float> minPropagationRadius{"minPropagationDistance", o2::constants::geom::XTPCInnerRef + 0.1, "Only tracks which are at a smaller radius will be propagated, defaults to TPC inner wall"};

Check failure on line 74 in Common/Tools/TrackPropagationModule.h

View workflow job for this annotation

GitHub Actions / O2 linter

[name/configurable]

Use lowerCamelCase for names of configurables and use the same name for the struct member as for the JSON string. (Declare the type and names on the same line.)
// for TrackTuner only (MC smearing)
o2::framework::Configurable<bool> useTrackTuner{"useTrackTuner", false, "Apply track tuner corrections to MC"};
o2::framework::Configurable<bool> useTrkPid{"useTrkPid", false, "use pid in tracking"};
Expand Down Expand Up @@ -176,7 +176,7 @@
/// to understand whether the TrackTuner::getDcaGraphs function can be called here (input path from string/configurables)
/// or inside the process function, to "auto-detect" the input file based on the run number
const auto& workflows = initContext.services().template get<o2::framework::RunningWorkflowInfo const>();
for (const o2::framework::DeviceSpec& device : workflows.devices) { /// loop over devices

Check failure on line 179 in Common/Tools/TrackPropagationModule.h

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop]

Use constant references for non-modified iterators in range-based for loops.
if (device.name == "propagation-service") {
// loop over the options
// to find the value of TrackTuner::autoDetectDcaCalib
Expand Down Expand Up @@ -271,8 +271,9 @@
// std::array<float, 3> trackPxPyPz;
// std::array<float, 3> trackPxPyPzTuned = {0.0, 0.0, 0.0};
double q2OverPtNew = -9999.;
bool isPropagationOK = true;
// Only propagate tracks which have passed the innermost wall of the TPC (e.g. skipping loopers etc). Others fill unpropagated.
if (track.trackType() == o2::aod::track::TrackIU && track.x() < cGroup.minPropagationRadius.value) {
if (track.x() < cGroup.minPropagationRadius.value) {
if (fillTracksCov) {
if constexpr (isMc) { // checking MC and fillCovMat block begins
// bool hasMcParticle = track.has_mcParticle();
Expand All @@ -287,29 +288,51 @@
}
} // MC and fillCovMat block ends
}
bool isPropagationOK = true;

if (track.has_collision()) {
auto const& collision = collisions.rawIteratorAt(track.collisionId());
if (fillTracksCov) {
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov);

if (track.trackType() == o2::aod::track::TrackIU) {
if (track.has_collision()) {
auto const& collision = collisions.rawIteratorAt(track.collisionId());
if (fillTracksCov) {
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov);
} else {
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo);
}
} else {
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo);
if (fillTracksCov) {
mVtx.setPos({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()});
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());
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov);
} else {
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo);
}
}
if (isPropagationOK) {
trackType = o2::aod::track::Track;
}
} else {
if (fillTracksCov) {
mVtx.setPos({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()});
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());
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov);
} else {
isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo);
if (fillTracksDCA || fillTracksDCACov) {
if (track.has_collision()) {
auto const& collision = collisions.rawIteratorAt(track.collisionId());
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
if (fillTracksDCACov) {
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
}
} else {
mVtx.setPos({ccdbLoader.mMeanVtx->getX(), ccdbLoader.mMeanVtx->getY(), ccdbLoader.mMeanVtx->getZ()});
if (fillTracksDCACov) {
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());
}
}

if (fillTracksDCACov) {
calculateDCA(mTrackParCov, mVtx, o2::base::Propagator::Instance()->getNominalBz(), &mDcaInfoCov, 999.f);
} else {
calculateDCA(mTrackPar, mVtx, o2::base::Propagator::Instance()->getNominalBz(), &mDcaInfo, 999.f);
}
}
}
if (isPropagationOK) {
trackType = o2::aod::track::Track;
}
// filling some QA histograms for track tuner test purpose
if (fillTracksCov) {
if constexpr (isMc) { // checking MC and fillCovMat block begins
Expand Down Expand Up @@ -356,6 +379,67 @@
}
}
}

//_______________________________________________________________________
// TTrackPar type: either aod::TrackPar or aod::TrackParCov
// TVertex type: either math_utils::Point3D<value_t> or o2::dataformats::VertexBase
// TDCA type: either dim2_t or o2::dataformats::DCA
template <typename TTrackPar, typename TVertex, typename TDCA>
bool calculateDCA(TTrackPar trackPar, const TVertex& vtx, double b, TDCA* dca, double maxD)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the reason to pass the trackPar by copy rather than by reference?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @shahor02 ! That's a mistake from my side, thank you for pointing that out! I will change it tomorrow morning

{
// propagate track to DCA to the vertex
double sn, cs, alp = trackPar.getAlpha();
o2::math_utils::detail::sincos(alp, sn, cs);
double x = trackPar.getX(), y = trackPar.getY(), snp = trackPar.getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp));
double xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
x -= xv;
y -= yv;
// Estimate the impact parameter neglecting the track curvature
double d = gpu::CAMath::Abs(x * snp - y * csp);
if (d > maxD) {
// provide default DCA for failed propag
if constexpr (requires { trackPar.getSigmaY2(); vtx.getSigmaX2(); }) {
dca->set(o2::track::DefaultDCA, o2::track::DefaultDCA,
o2::track::DefaultDCACov, o2::track::DefaultDCACov, o2::track::DefaultDCACov);
} else {
(*dca)[0] = o2::track::DefaultDCA;
(*dca)[1] = o2::track::DefaultDCA;
}
return false;
}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I understand correctly, the assumption is that the track is already propagated to the PCA to the vertex and just needs a rotation to the frame where the track is normal to the DCA vector. Then (1) it is not clear why the code above is needed (2) how the track can be in another alpha frame if it was already propagated to the DCA.
If it is already in that frame, then one just rotated the PV to the same frame, i.e. L394 is enough.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function is essentially a generalization of the propagateToDCABxByBz in the Detectors/Base/src/Propagator.cxx for TrackPar and TrackParCov, where the propagation call has been removed. All the rest remained the same.
I agree that l.398-409 and l.421-432 can be removed. However, shouldn't we keep the remaining parts for the DCA covariance calculation?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the track was already propagated to the PCA and rotated to the frame normal to the DCA vector, then the only code needed is for the rotation of the PV to the same frame, i.e.

float sn, cs, alp = trackPar.getAlpha();
o2::math_utils::detail::sincos(alp, sn, cs);
float xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
if constexpr (requires { trackPar.getSigmaY2(); vtx.getSigmaX2(); }) {
  o2::math_utils::detail::sincos(alp, sn, cs);
  auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
  dca->set(trackPar.getY() - yv, trackPar.getZ() - zv, trackPar.getSigmaY2() + s2ylocvtx, trackPar.getSigmaZY(), trackPar.getSigmaZ2() + vtx.getSigmaZ2());
} else {
  (*dca)[0] = trackPar.getY() - yv;
  (*dca)[1] = trackPar.getZ() - zv;
}

If the track was propagated to the PCA but not rotated (but I don't see under which circumstances this may be the case, unless there was some rotation applied after the initial propagateToDCA), then the full code is indeed needed.
And in all other case (including track modification by the MC tuner mentioned by @mfaggin ) the full propagation and rotation will be needed (though there is no reason to do this in double).

double crv = trackPar.getCurvature(b);
double tgfv = -(crv * x - snp) / (crv * y + csp);
sn = tgfv / gpu::CAMath::Sqrt(1.f + tgfv * tgfv);
cs = gpu::CAMath::Sqrt((1.f - sn) * (1.f + sn));
cs = (gpu::CAMath::Abs(tgfv) > constants::math::Almost0) ? sn / tgfv : constants::math::Almost1;

x = xv * cs + yv * sn;
yv = -xv * sn + yv * cs;
xv = x;

alp += gpu::CAMath::ASin(sn);
if (!trackPar.rotate(alp)) {
Comment on lines +418 to +419
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Did you simply suppress these 2 lines? If the track would not be rotated to the frame of the PCA after the propagation, the difference would be much larger than what your pdf shows.
Is it possible that you are related the same preprogated track to two different compatible PVs?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I actually tried different approaches:

  • I considered only the rotation of the PV to the same frame as the track at the PCA (using the exactly the code you suggested 3days ago, i.e.
float sn, cs, alp = trackPar.getAlpha();
o2::math_utils::detail::sincos(alp, sn, cs);
float xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
if constexpr (requires { trackPar.getSigmaY2(); vtx.getSigmaX2(); }) {
  o2::math_utils::detail::sincos(alp, sn, cs);
  auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
  dca->set(trackPar.getY() - yv, trackPar.getZ() - zv, trackPar.getSigmaY2() + s2ylocvtx, trackPar.getSigmaZY(), trackPar.getSigmaZ2() + vtx.getSigmaZ2());
} else {
  (*dca)[0] = trackPar.getY() - yv;
  (*dca)[1] = trackPar.getZ() - zv;
}

)

  • I also tried removing the following line entirely:
if (!trackPar.rotate(alp)) {

In both cases, I got the same DCA distribution which corresponds to the one shown in open star markers in the PDF attached in my previous message.

I suppose it is possible that I am calculating the DCA of a prepropagated track with respect to another compatible PV. Is there a way to check that?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see that some analyses are using track.compatibleCollIds().size() > 0 to check if the track is ambiguous, but I don't know if getting this track attribute requires some special. Btw, in my code snippet, which you quoted, I've wrongly left unnecessary sincos call, should be

float sn, cs, alp = trackPar.getAlpha();
o2::math_utils::detail::sincos(alp, sn, cs);
float xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
if constexpr (requires { trackPar.getSigmaY2(); vtx.getSigmaX2(); }) {
  auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
  dca->set(trackPar.getY() - yv, trackPar.getZ() - zv, trackPar.getSigmaY2() + s2ylocvtx, trackPar.getSigmaZY(), trackPar.getSigmaZ2() + vtx.getSigmaZ2());
} else {
  (*dca)[0] = trackPar.getY() - yv;
  (*dca)[1] = trackPar.getZ() - zv;
}

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @shahor02 ! Apologies for the late reply!
I tried the correct code snippet this time (without the remaining unnecessary sincos call) but the DCA distribution looked the same as before, meaning that I still obtain the open star markers in the PDF attached in my previous message.

Also, I repeated the comparison focusing only on collisions without any other compatible PV (by requiring track.compatibleCollIds().size() > 0) and now the DCA distribution with your code snipped (without rotation) looks the same as with rotation and the same as without prepropagated tracks.
DCAxy_StdVsPrepropagate_WithoutRotation_withoutCompatibleColl_Strictly0.pdf

It seems that, indeed, I was calculating the DCA of some prepropagated tracks with respect to another compatible PV.
However, I am not entirely sure to understand how that's possible. I am wondering if it is not possible that the tracks are always prepropagated to the mean vertex (https://github.com/AliceO2Group/AliceO2/blob/4090041b401c7aa6c919ca923126fff950cbccd1/Detectors/AOD/src/AODProducerWorkflowSpec.cxx#L2858) instead of the PV of its assigned collision. If I am saying something wrong, please don't hesitate to correct me

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @romainschotter ,
The extra sincos could not change the results, it was simply not needed.

Concerning this:

Also, I repeated the comparison focusing only on collisions without any other compatible PV (by requiring track.compatibleCollIds().size() > 0) and now the DCA distribution with your code snipped (without rotation) looks the same as with rotation and the same as without prepropagated tracks.

Do you mean track.compatibleCollIds().size() == 1 (for the tracks compatible with 1 PV only)?

The track is propagated to the meanvertex only if it is not associated to any real PV. But in case it is associated to >1 PV, the propagation will be done only to 1st one (since the track is saved only once).

But if this is the case, then in your method just rotation would not be enough;
in case of the ambiguous vertex assignment also the propagation would be needed, as it is done in case of trackTuning.
I am wondering if your reference distribution is not suffering from a similar problem.
I think it is worth comparing for ambiguous tracks only (i.e. track.compatibleCollIds().size() > 0) the DCA distribution with old propagate to dca, and your method with rotation only and with rotation + propagation.

// provide default DCA for failed propag
if constexpr (requires { trackPar.getSigmaY2(); vtx.getSigmaX2(); }) {
dca->set(o2::track::DefaultDCA, o2::track::DefaultDCA,
o2::track::DefaultDCACov, o2::track::DefaultDCACov, o2::track::DefaultDCACov);
} else {
(*dca)[0] = o2::track::DefaultDCA;
(*dca)[1] = o2::track::DefaultDCA;
}

return false;
}
if constexpr (requires { trackPar.getSigmaY2(); vtx.getSigmaX2(); }) {
o2::math_utils::detail::sincos(alp, sn, cs);
auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
dca->set(trackPar.getY() - yv, trackPar.getZ() - zv, trackPar.getSigmaY2() + s2ylocvtx, trackPar.getSigmaZY(), trackPar.getSigmaZ2() + vtx.getSigmaZ2());
} else {
(*dca)[0] = trackPar.getY() - yv;
(*dca)[1] = trackPar.getZ() - zv;
}
return true;
}
};

} // namespace common
Expand Down
Loading