diff --git a/Detectors/Align/include/Align/AlignConfig.h b/Detectors/Align/include/Align/AlignConfig.h index 91b503c2c923e..e72d436a14e3b 100644 --- a/Detectors/Align/include/Align/AlignConfig.h +++ b/Detectors/Align/include/Align/AlignConfig.h @@ -85,6 +85,7 @@ struct AlignConfig : public o2::conf::ConfigurableParamHelper { float controlFraction = -1.; // fraction for which control output is requested, if negative - only 1st instance of device will write them float MPRecOutFraction = -1.; // compact Millepede2Record fraction, if negative - only 1st instance of device will write them + bool useLinRef = true; // use initial track for lienarization reference point bool MilleOut = true; // Mille output bool KalmanResid = true; // Kalman residuals bool MilleOutBin = true; // text vs binary output for mille data diff --git a/Detectors/Align/include/Align/AlignmentTrack.h b/Detectors/Align/include/Align/AlignmentTrack.h index ef4552cb9a37d..cb69f11cbf85c 100644 --- a/Detectors/Align/include/Align/AlignmentTrack.h +++ b/Detectors/Align/include/Align/AlignmentTrack.h @@ -39,6 +39,7 @@ class AlignmentTrack : public trackParam_t, public TObject { public: using trackParam_t = o2::track::TrackParametrizationWithError; + using trackPar_t = o2::track::TrackParametrization; using PropagatorD = o2::base::PropagatorD; using MatCorrType = PropagatorD::MatCorrType; using GTrackID = o2::dataformats::GlobalTrackID; @@ -83,9 +84,9 @@ class AlignmentTrack : public trackParam_t, public TObject // template void copyFrom(const o2::track::TrackParametrizationWithError

& trc); - bool propagateToPoint(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tLT = nullptr, int signCorr = 0); - bool propagateParamToPoint(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep = 3, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, int signCorr = 0); // param only - bool propagateParamToPoint(trackParam_t* trSet, int nTr, const AlignmentPoint* pnt, double maxStep = 3, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, int signCorr = 0); // params only + bool propagateToPoint(trackParam_t& tr, trackPar_t* linRef, const AlignmentPoint* pnt, double maxStep, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tLT = nullptr, int signCorr = 0); + bool propagateParamToPoint(trackPar_t& tr, const AlignmentPoint* pnt, double maxStep = 3, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, int signCorr = 0); // param only + bool propagateParamToPoint(trackPar_t* trSet, int nTr, const AlignmentPoint* pnt, double maxStep = 3, double maxSnp = 0.95, MatCorrType mt = MatCorrType::USEMatCorrLUT, int signCorr = 0); // params only // bool calcResiduals(const double* params = nullptr); bool calcResidDeriv(double* params = nullptr); @@ -119,23 +120,23 @@ class AlignmentTrack : public trackParam_t, public TObject void imposePtBOff(double pt) { setQ2Pt(1. / pt); } // propagation methods void copyFrom(const trackParam_t* etp); - bool applyMatCorr(trackParam_t& trPar, const double* corrDiag, const AlignmentPoint* pnt); - bool applyMatCorr(trackParam_t* trSet, int ntr, const double* corrDiaf, const AlignmentPoint* pnt); - bool applyMatCorr(trackParam_t& trPar, const double* corrpar); + bool applyMatCorr(trackPar_t& trPar, const double* corrDiag, const AlignmentPoint* pnt); + bool applyMatCorr(trackPar_t* trSet, int ntr, const double* corrDiaf, const AlignmentPoint* pnt); + bool applyMatCorr(trackPar_t& trPar, const double* corrpar); // double getResidual(int dim, int pntID) const { return mResid[dim][pntID]; } const double* getDResDLoc(int dim, int pntID) const { return mDResDLoc[dim].data() + (pntID * mNLocPar); } const double* getDResDGlo(int dim, int id) const { return mDResDGlo[dim].data() + id; } const int* getGloParID() const { return mGloParID.data(); } // - void setParams(trackParam_t& tr, double x, double alp, const double* par, bool add); - void setParams(trackParam_t* trSet, int ntr, double x, double alp, const double* par, bool add); - void setParam(trackParam_t& tr, int par, double val); - void setParam(trackParam_t* trSet, int ntr, int par, double val); - void modParam(trackParam_t& tr, int par, double delta); - void modParam(trackParam_t* trSet, int ntr, int par, double delta); + void setParams(trackPar_t& tr, double x, double alp, const double* par, bool add); + void setParams(trackPar_t* trSet, int ntr, double x, double alp, const double* par, bool add); + void setParam(trackPar_t& tr, int par, double val); + void setParam(trackPar_t* trSet, int ntr, int par, double val); + void modParam(trackPar_t& tr, int par, double delta); + void modParam(trackPar_t* trSet, int ntr, int par, double delta); // - void richardsonDeriv(const trackParam_t* trSet, const double* delta, + void richardsonDeriv(const trackPar_t* trSet, const double* delta, const AlignmentPoint* pnt, double& derY, double& derZ); // const double* getLocPars() const { return mLocPar.data(); } @@ -179,13 +180,14 @@ class AlignmentTrack : public trackParam_t, public TObject std::vector mLocPar; // local parameters array std::vector mGloParID; // IDs of relevant global params private: - bool propagate(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr = 0); + bool propagate(trackParam_t& tr, trackPar_t* linRef, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr = 0); + bool propagate(trackPar_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr = 0); // ClassDefOverride(AlignmentTrack, 2) }; //____________________________________________________________________________________________ -inline void AlignmentTrack::setParams(trackParam_t& tr, double x, double alp, const double* par, bool add) +inline void AlignmentTrack::setParams(trackPar_t& tr, double x, double alp, const double* par, bool add) { // set track params const double kDefQ2PtCosm = 1; @@ -205,7 +207,7 @@ inline void AlignmentTrack::setParams(trackParam_t& tr, double x, double alp, co } //____________________________________________________________________________________________ -inline void AlignmentTrack::setParams(trackParam_t* trSet, int ntr, double x, double alp, const double* par, bool add) +inline void AlignmentTrack::setParams(trackPar_t* trSet, int ntr, double x, double alp, const double* par, bool add) { // set parames for multiple tracks (VECTORIZE THIS) if (!add) { // full parameter supplied @@ -224,14 +226,14 @@ inline void AlignmentTrack::setParams(trackParam_t* trSet, int ntr, double x, do } //____________________________________________________________________________________________ -inline void AlignmentTrack::setParam(trackParam_t& tr, int par, double val) +inline void AlignmentTrack::setParam(trackPar_t& tr, int par, double val) { // set track parameter tr.setParam(val, par); } //____________________________________________________________________________________________ -inline void AlignmentTrack::setParam(trackParam_t* trSet, int ntr, int par, double val) +inline void AlignmentTrack::setParam(trackPar_t* trSet, int ntr, int par, double val) { // set parames for multiple tracks (VECTORIZE THIS) for (int i = 0; i < ntr; ++i) { @@ -240,7 +242,7 @@ inline void AlignmentTrack::setParam(trackParam_t* trSet, int ntr, int par, doub } //____________________________________________________________________________________________ -inline void AlignmentTrack::modParam(trackParam_t& tr, int par, double delta) +inline void AlignmentTrack::modParam(trackPar_t& tr, int par, double delta) { // modify track parameter const auto val = tr.getParam(par) + delta; @@ -248,7 +250,7 @@ inline void AlignmentTrack::modParam(trackParam_t& tr, int par, double delta) } //____________________________________________________________________________________________ -inline void AlignmentTrack::modParam(trackParam_t* trSet, int ntr, int par, double delta) +inline void AlignmentTrack::modParam(trackPar_t* trSet, int ntr, int par, double delta) { // modify track parameter (VECTORIZE THOS) for (int i = 0; i < ntr; ++i) { diff --git a/Detectors/Align/src/AlignmentTrack.cxx b/Detectors/Align/src/AlignmentTrack.cxx index 554d30e246e29..644ee07c64984 100644 --- a/Detectors/Align/src/AlignmentTrack.cxx +++ b/Detectors/Align/src/AlignmentTrack.cxx @@ -168,7 +168,7 @@ bool AlignmentTrack::calcResidDeriv(double* extendedParams, bool invert, int pFr // (like http://root.cern.ch/root/html/ROOT__Math__RichardsonDerivator.html) // const auto& algConf = AlignConfig::Instance(); - trackParam_t probD[kNRDClones]; // use this to vary supplied param for derivative calculation + trackPar_t probD[kNRDClones]; // use this to vary supplied param for derivative calculation double varDelta[kRichardsonN]; const int kInvElem[kNKinParBON] = {-1, 1, 1, -1, -1}; // @@ -511,7 +511,7 @@ bool AlignmentTrack::calcResiduals(const double* extendedParams, bool invert, in } //______________________________________________________ -bool AlignmentTrack::propagateParamToPoint(trackParam_t* tr, int nTr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, int signCorr) +bool AlignmentTrack::propagateParamToPoint(trackPar_t* tr, int nTr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, int signCorr) { // Propagate set of tracks to the point (only parameters, no error matrix) // VECTORIZE this @@ -521,7 +521,7 @@ bool AlignmentTrack::propagateParamToPoint(trackParam_t* tr, int nTr, const Alig if (!propagateParamToPoint(tr[itr], pnt, maxStep, maxSnp, mt, signCorr)) { if (algConf.verbose > 2) { LOG(error) << "Failed on clone " << itr << " propagation "; - tr[itr].print(); + tr[itr].printParam(); pnt->print(AlignmentPoint::kMeasurementBit | AlignmentPoint::kMaterialBit); } return false; @@ -531,21 +531,33 @@ bool AlignmentTrack::propagateParamToPoint(trackParam_t* tr, int nTr, const Alig } //______________________________________________________ -bool AlignmentTrack::propagateParamToPoint(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, int signCorr) +bool AlignmentTrack::propagateParamToPoint(trackPar_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, int signCorr) { // propagate tracks to the point (only parameters, no error matrix) return propagate(tr, pnt, maxStep, maxSnp, mt, nullptr, signCorr); } //______________________________________________________ -bool AlignmentTrack::propagateToPoint(trackParam_t& tr, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr) +bool AlignmentTrack::propagateToPoint(trackParam_t& tr, trackPar_t* linRef, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr) { // propagate tracks to the point. If matCor is true, then material corrections will be applied. // if matPar pointer is provided, it will be filled by total x2x0 and signed xrho - return propagate(tr, pnt, maxStep, maxSnp, mt, tLT, signCorr); + return propagate(tr, linRef, pnt, maxStep, maxSnp, mt, tLT, signCorr); } -bool AlignmentTrack::propagate(trackParam_t& track, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr) +bool AlignmentTrack::propagate(trackParam_t& track, trackPar_t* linRef, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr) +{ + if (signCorr == 0) { // auto + // calculate the sign of the energy loss correction and ensure the upper leg of cosmics is calculated correctly. + double dx = pnt->getXTracking() - track.getX(); + int dir = dx > 0.f ? 1 : -1; + signCorr = pnt->isInvDir() ? dir : -dir; // propagation along the track direction should have signCorr=-1 + } + // do propagation in at least 2 step to reveal eventual effect of MS on the position + return PropagatorD::Instance()->propagateToAlphaX(track, linRef, pnt->getAlphaSens(), pnt->getXTracking(), pnt->getUseBzOnly(), maxSnp, maxStep, 2, mt, tLT, signCorr); +} + +bool AlignmentTrack::propagate(trackPar_t& track, const AlignmentPoint* pnt, double maxStep, double maxSnp, MatCorrType mt, track::TrackLTIntegral* tLT, int signCorr) { if (signCorr == 0) { // auto // calculate the sign of the energy loss correction and ensure the upper leg of cosmics is calculated correctly. @@ -603,7 +615,7 @@ bool AlignmentTrack::ApplyMS(trackParam_t& trPar, double tms,double pms) */ //______________________________________________________ -bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corrPar, const AlignmentPoint* pnt) +bool AlignmentTrack::applyMatCorr(trackPar_t& trPar, const double* corrPar, const AlignmentPoint* pnt) { // Modify track param (e.g. trackParam_t) in the tracking frame // by delta accounting for material effects @@ -630,7 +642,7 @@ bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corrPar, co } //______________________________________________________ -bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corr) +bool AlignmentTrack::applyMatCorr(trackPar_t& trPar, const double* corr) { // Modify track param (e.g. trackParam_t) in the tracking frame // by delta accounting for material effects @@ -645,7 +657,7 @@ bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corr) printf("%+.3e ", corr[i]); } printf("\n"); - trPar.print(); + trPar.printParam(); } return false; } @@ -656,7 +668,7 @@ bool AlignmentTrack::applyMatCorr(trackParam_t& trPar, const double* corr) } //______________________________________________________ -bool AlignmentTrack::applyMatCorr(trackParam_t* trSet, int ntr, const double* corrDiag, const AlignmentPoint* pnt) +bool AlignmentTrack::applyMatCorr(trackPar_t* trSet, int ntr, const double* corrDiag, const AlignmentPoint* pnt) { // Modify set of track params (e.g. trackParam_t) in the tracking frame // by delta accounting for material effects @@ -683,7 +695,7 @@ bool AlignmentTrack::applyMatCorr(trackParam_t* trSet, int ntr, const double* co if (!applyMatCorr(trSet[itr], corr)) { if (algConf.verbose > 2) { LOGP(error, "Failed on clone {} materials", itr); - trSet[itr].print(); + trSet[itr].printParam(); } return false; } @@ -732,7 +744,7 @@ double AlignmentTrack::richardsonExtrap(const double* val, int ord) } //______________________________________________ -void AlignmentTrack::richardsonDeriv(const trackParam_t* trSet, const double* delta, const AlignmentPoint* pnt, double& derY, double& derZ) +void AlignmentTrack::richardsonDeriv(const trackPar_t* trSet, const double* delta, const AlignmentPoint* pnt, double& derY, double& derZ) { // Calculate Richardson derivatives for diagonalized Y and Z from a set of kRichardsonN pairs // of tracks with same parameter of i-th pair varied by +-delta[i] @@ -882,7 +894,7 @@ bool AlignmentTrack::iniFit() // // propagate to reference point, which is the inner point of lower leg const AlignmentPoint* refP = getPoint(getInnerPointID()); - if (!propagateToPoint(trcU, refP, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, -1)) { // moving along the track: energy is lost + if (!propagateToPoint(trcU, nullptr, refP, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, -1)) { // moving along the track: energy is lost return false; } // @@ -1024,6 +1036,7 @@ bool AlignmentTrack::fitLeg(trackParam_t& trc, int pFrom, int pTo, bool& inv) } return false; } + trackPar_t linRef(trc), *linRefP = algConf.useLinRef ? &linRef : nullptr; trc.setCov(kIniErr); trc.setCov(16 * trc.getQ2Pt() * trc.getQ2Pt(), 4, 4); // lowest diagonal element (Q2Pt2) // @@ -1042,7 +1055,7 @@ bool AlignmentTrack::fitLeg(trackParam_t& trc, int pFrom, int pTo, bool& inv) int pntCnt = 0; for (int ip = pFrom; ip != pTo; ip += pinc) { // inward fit from outer point AlignmentPoint* pnt = getPoint(ip); - if (!propagateToPoint(trc, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { // against track direction : e.loss is compensated + if (!propagateToPoint(trc, linRefP, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { // against track direction : e.loss is compensated if (algConf.verbose > 2) { LOGF(warn, "Failed on propagateToPoint %d (%d : %d) %f", ip, pFrom, pTo, pnt->getXTracking()); trc.print(); @@ -1139,7 +1152,7 @@ bool AlignmentTrack::residKalman() trc.invert(); inv = !inv; } - if (!propagateToPoint(trc, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { + if (!propagateToPoint(trc, nullptr, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { return false; } if (!pnt->containsMeasurement()) { @@ -1178,7 +1191,7 @@ bool AlignmentTrack::residKalman() trc.invert(); inv = !inv; } - if (!propagateToPoint(trc, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { // we are going along track direction, e.loss is applied + if (!propagateToPoint(trc, nullptr, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), nullptr, signELoss)) { // we are going along track direction, e.loss is applied return false; } if (!pnt->containsMeasurement()) { @@ -1335,7 +1348,7 @@ bool AlignmentTrack::processMaterials(trackParam_t& trc, int pFrom, int pTo) // matTL.clearFast(); // printf("-> ProcMat %d (%d->%d)\n",ip,pFrom,pTo); - if (!propagateToPoint(trc, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), &matTL, signELoss)) { // with material corrections + if (!propagateToPoint(trc, nullptr, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType(algConf.matCorType), &matTL, signELoss)) { // with material corrections if (algConf.verbose > 2) { LOG(error) << "Failed to take track to point" << ip << " (dir: " << pFrom << "->" << pTo << ") with mat.corr."; trc.print(); @@ -1346,7 +1359,7 @@ bool AlignmentTrack::processMaterials(trackParam_t& trc, int pFrom, int pTo) // // is there enough material to consider the point as a scatterer? bool hasMaterial = matTL.getX2X0() > minX2X0; - if (!propagateToPoint(tr0, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType::USEMatCorrNONE, nullptr, signELoss)) { // no material corrections + if (!propagateToPoint(tr0, nullptr, pnt, algConf.maxStep, algConf.maxSnp, MatCorrType::USEMatCorrNONE, nullptr, signELoss)) { // no material corrections if (algConf.verbose > 2) { LOG(error) << "Failed to take track to point" << ip << " (dir: " << pFrom << "->" << pTo << ") with mat.corr."; tr0.print(); diff --git a/Detectors/Base/include/DetectorsBase/Propagator.h b/Detectors/Base/include/DetectorsBase/Propagator.h index 6fa750577255d..75b9446aebade 100644 --- a/Detectors/Base/include/DetectorsBase/Propagator.h +++ b/Detectors/Base/include/DetectorsBase/Propagator.h @@ -127,6 +127,9 @@ class PropagatorImpl GPUd() bool propagateToAlphaX(track_T& track, value_type alpha, value_type x, bool bzOnly = false, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP, int minSteps = 1, MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const; + GPUd() bool propagateToAlphaX(TrackParCov_t& track, TrackPar_t* linRef, value_type alpha, value_type x, bool bzOnly = false, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP, int minSteps = 1, + MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const; + template GPUd() bool propagateToR(track_T& track, value_type r, bool bzOnly = false, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP, MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const; diff --git a/Detectors/Base/src/Propagator.cxx b/Detectors/Base/src/Propagator.cxx index 02e7a05080ac5..a5983cab8e257 100644 --- a/Detectors/Base/src/Propagator.cxx +++ b/Detectors/Base/src/Propagator.cxx @@ -648,6 +648,23 @@ GPUd() bool PropagatorImpl::propagateToR(track_T& track, value_type r, return propagateToX(track, xfin, bzOnly, maxSnp, maxStep, matCorr, tofInfo, signCorr); } +template +GPUd() bool PropagatorImpl::propagateToAlphaX(TrackParCov_t& track, TrackPar_t* linRef, value_type alpha, value_type x, bool bzOnly, value_type maxSnp, value_type maxStep, int minSteps, + MatCorrType matCorr, track::TrackLTIntegral* tofInfo, int signCorr) const +{ + // propagate to alpha,X, if needed in a few steps + auto snp = track.getSnpAt(alpha, x, getNominalBz()); + // apply safety factor 0.9 for crude rotation estimate + if (math_utils::detail::abs(snp) < maxSnp * 0.9 && (linRef ? track.rotate(alpha, *linRef, getNominalBz()) : track.rotate(alpha))) { + auto dx = math_utils::detail::abs(x - track.getX()); + if (dx < Epsilon) { + return true; + } + return propagateTo(track, linRef, x, bzOnly, maxSnp, math_utils::detail::min(dx / minSteps, maxStep), matCorr, tofInfo, signCorr); + } + return false; +} + //_______________________________________________________________________ template template