diff --git a/Common/DCAFitter/CMakeLists.txt b/Common/DCAFitter/CMakeLists.txt deleted file mode 100644 index e15c2885f2259..0000000000000 --- a/Common/DCAFitter/CMakeLists.txt +++ /dev/null @@ -1,40 +0,0 @@ -# Copyright 2019-2020 CERN and copyright holders of ALICE O2. -# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -# All rights not expressly granted are reserved. -# -# This software is distributed under the terms of the GNU General Public -# License v3 (GPL Version 3), copied verbatim in the file "COPYING". -# -# In applying this license CERN does not waive the privileges and immunities -# granted to it by virtue of its status as an Intergovernmental Organization -# or submit itself to any jurisdiction. - -o2_add_library(DCAFitter - TARGETVARNAME targetName - SOURCES src/DCAFitterN.cxx - src/FwdDCAFitterN.cxx - PUBLIC_LINK_LIBRARIES ROOT::Core - O2::CommonUtils - O2::ReconstructionDataFormats - O2::DataFormatsParameters - O2::DetectorsBase) - -o2_target_root_dictionary(DCAFitter - HEADERS include/DCAFitter/HelixHelper.h - include/DCAFitter/DCAFitterN.h - include/DCAFitter/FwdDCAFitterN.h) - -if (OpenMP_CXX_FOUND) - target_compile_definitions(${targetName} PRIVATE WITH_OPENMP) - target_link_libraries(${targetName} PRIVATE OpenMP::OpenMP_CXX) -endif() - - -o2_add_test( - DCAFitterN - SOURCES test/testDCAFitterN.cxx - COMPONENT_NAME DCAFitter - PUBLIC_LINK_LIBRARIES O2::DCAFitter ROOT::Core ROOT::Physics - LABELS vertexing - ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage - VMCWORKDIR=${CMAKE_BINARY_DIR}/stage/${CMAKE_INSTALL_DATADIR}) diff --git a/Common/DCAFitter/README.md b/Common/DCAFitter/README.md deleted file mode 100644 index c952a2e33124f..0000000000000 --- a/Common/DCAFitter/README.md +++ /dev/null @@ -1,80 +0,0 @@ - - -## DCAFitterN - -Templated class to fit the Point of Closest Approach (PCA) of secondary vertex with N prongs. Allows minimization of either absolute or weighted Distances of Closest Approach (DCA) of N tracks to their common PCA. - -For every N (prongs) a separate specialization must be instantiated, e.g. -```cpp -using Track = o2::track::TrackParCov; -o2::vertexing::DCAFitterN<2,Track> ft2; // 2-prongs fitter -// or, to set at once some parameters -float bz = 5.; // field in kGauss -bool useAbsDCA = true; // use abs. DCA minimizaition instead of default weighted -bool propToDCA = true; // after fit, create a copy of tracks at the found PCA -o2::vertexing::DCAFitterN<3,Track> ft3(bz, useAbsDCA, propToDCA); // 3-prongs fitter -``` -One can also use predefined aliases ``o2::vertexing::DCAFitter2`` and ``o2::vertexing::DCAFitter3``; -The main processing method is -```cpp -o2::vertexing::DCAFitterN::process(const Track& trc1,..., cons Track& trcN); -``` - -The typical use case is (for e.g. 3-prong fitter): -```cpp -using Vec3D = ROOT::Math::SVector; // this is a type of the fitted vertex -o2::vertexing::DCAFitter3 ft; -ft.setBz(5.0); -ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway -ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway -ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway -ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway -ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway -ft.setMinRelChi2Change(0.9);// stop iterations if chi2 improves by less that this factor -ft.setMaxChi2(10); // discard vertices with chi2/Nprongs (or sum{DCAi^2}/Nprongs for abs. distance minimization) - -Track tr0,tr1,tr2; // decide candidate tracks -int nc = ft.process(tr0,tr1,tr2); // one can have up to 2 candidates, though the 2nd (if any) will have worse quality -if (nc) { - Vec3D vtx = ft.getPCACandidate(); // same as ft.getPCACandidate(0); - LOG(info) << "found vertex " << vtx[0] << ' ' << vtx[1] << ' ' << vtx[2]; - // access the track's X parameters at PCA - for (int i=0;i<3;i++) { - LOG(info) << "Track " << i << " at PCA for X = " << ft.getTrackX(i); - } - // access directly the tracks propagated to the DCA - for (int i=0;i<3;i++) { - const auto& track = ft.getTrack(i); - track.print(); - } -} -``` - -By default the propagation is done with bZ provided by the user and w/o material corrections applied. -One can request the propagation with full local field and/or material corrections by setting -``` -ft.setUsePropagator(true); // use must take care of initialization of the propagator (loading geometry and magnetic field) -ft.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrLUT); // of USEMatCorrTGeo -``` - -Note that if material correction is not default USEMatCorrNone, then the propagator will be used even if not requested (hence must be initialized by the user). - -To get the most precise results one can request `ft.setRefitWithMatCorr(true)`: in this case when `propagateTracksToVertex()` is called, the tracks will be propagated -to the V0 with requested material corrections, one new V0 minimization will be done and only after that the final propagation to final V0 position will be done. -Since this is CPU consiming, it is reccomended to disable propagation to V0 by default (`ft.setPropagateToPCA(false)`) and call separately `ft.propagateTracksToVertex()` -after preliminary checks on the V0 candidate. - -By default the final V0 position is defined as -1) With `useAbsDCA = true`: simple average of tracks position propagated to respective `X_dca` parameters and rotated to the lab. frame. -2) With `useAbsDCA = false`: weighted (by tracks covariances) average of tracks position propagated to respective `X_dca` parameters and rotated to the lab. frame. - -Extra method `setWeightedFinalPCA(bool)` is provided for the "mixed" mode: if `setWeightedFinalPCA(true)` is set with `useAbsDCA = true` before the `process` call, the minimization will be done neglecting the track covariances, -but the final V0 position will be calculated using weighted average. One can also recalculate the V0 position by the weighted average method by calling explicitly -`ft.recalculatePCAWithErrors(int icand=0)`, w/o prior call of `setWeightedFinalPCA(true)`: this will update the position returned by the `getPCACandidate(int cand = 0)`. - -The covariance matrix of the V0 position is calculated as an inversed sum of tracks inversed covariances at respective `X_dca` points. - -See ``O2/Detectors/Base/test/testDCAFitterN.cxx`` for more extended example. -Currently only 2 and 3 prongs permitted, thought this can be changed by modifying ``DCAFitterN::NMax`` constant. diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h deleted file mode 100644 index 70b5d820423d3..0000000000000 --- a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h +++ /dev/null @@ -1,1078 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file DCAFitterN.h -/// \brief Defintions for N-prongs secondary vertex fit -/// \author ruben.shahoyan@cern.ch -/// For the formulae derivation see /afs/cern.ch/user/s/shahoian/public/O2/DCAFitter/DCAFitterN.pdf - -#ifndef _ALICEO2_DCA_FITTERN_ -#define _ALICEO2_DCA_FITTERN_ -#include -#include "MathUtils/Cartesian.h" -#include "ReconstructionDataFormats/Track.h" -#include "DCAFitter/HelixHelper.h" -#include "DetectorsBase/Propagator.h" - -namespace o2 -{ -namespace vertexing -{ -///__________________________________________________________________________________ -///< Inverse cov matrix (augmented by a dummy X error) of the point defined by the track -struct TrackCovI { - float sxx, syy, syz, szz; - - TrackCovI(const o2::track::TrackParCov& trc, float xerrFactor = 1.) { set(trc, xerrFactor); } - - TrackCovI() = default; - - void set(const o2::track::TrackParCov& trc, float xerrFactor = 1) - { - // we assign Y error to X for DCA calculation - // (otherwise for quazi-collinear tracks the X will not be constrained) - float cyy = trc.getSigmaY2(), czz = trc.getSigmaZ2(), cyz = trc.getSigmaZY(), cxx = cyy * xerrFactor; - float detYZ = cyy * czz - cyz * cyz; - if (detYZ > 0.) { - auto detYZI = 1. / detYZ; - sxx = 1. / cxx; - syy = czz * detYZI; - syz = -cyz * detYZI; - szz = cyy * detYZI; - } else { - throw std::runtime_error("invalid track covariance"); - } - } -}; - -///__________________________________________________________________________ -///< Derivative (up to 2) of the TrackParam position over its running param X -struct TrackDeriv { - float dydx, dzdx, d2ydx2, d2zdx2; - TrackDeriv() = default; - TrackDeriv(const o2::track::TrackPar& trc, float bz) { set(trc, bz); } - void set(const o2::track::TrackPar& trc, float bz) - { - float snp = trc.getSnp(), csp = std::sqrt((1. - snp) * (1. + snp)), cspI = 1. / csp, crv2c = trc.getCurvature(bz) * cspI; - dydx = snp * cspI; // = snp/csp - dzdx = trc.getTgl() * cspI; // = tgl/csp - d2ydx2 = crv2c * cspI * cspI; // = crv/csp^3 - d2zdx2 = crv2c * dzdx * dydx; // = crv*tgl*snp/csp^3 - } -}; - -template -class DCAFitterN -{ - static constexpr double NMin = 2; - static constexpr double NMax = 4; - static constexpr double NInv = 1. / N; - static constexpr int MAXHYP = 2; - static constexpr float XerrFactor = 5.; // factor for conversion of track covYY to dummy covXX - using Track = o2::track::TrackParCov; - using TrackAuxPar = o2::track::TrackAuxPar; - using CrossInfo = o2::track::CrossInfo; - - using Vec3D = ROOT::Math::SVector; - using VecND = ROOT::Math::SVector; - using MatSym3D = ROOT::Math::SMatrix>; - using MatStd3D = ROOT::Math::SMatrix>; - using MatSymND = ROOT::Math::SMatrix>; - using MatStdND = ROOT::Math::SMatrix>; - using TrackCoefVtx = MatStd3D; - using ArrTrack = std::array; // container for prongs (tracks) at single vertex cand. - using ArrTrackCovI = std::array; // container for inv.cov.matrices at single vertex cand. - using ArrTrCoef = std::array; // container of TrackCoefVtx coefficients at single vertex cand. - using ArrTrDer = std::array; // container of Track 1st and 2nd derivative over their X param - using ArrTrPos = std::array; // container of Track positions - - public: - static constexpr int getNProngs() { return N; } - - DCAFitterN() = default; - DCAFitterN(float bz, bool useAbsDCA, bool prop2DCA) : mBz(bz), mUseAbsDCA(useAbsDCA), mPropagateToPCA(prop2DCA) - { - static_assert(N >= NMin && N <= NMax, "N prongs outside of allowed range"); - } - - //========================================================================= - ///< return PCA candidate, by default best on is provided (no check for the index validity) - const Vec3D& getPCACandidate(int cand = 0) const { return mPCA[mOrder[cand]]; } - const auto getPCACandidatePos(int cand = 0) const - { - const auto& vd = mPCA[mOrder[cand]]; - return std::array{float(vd[0]), float(vd[1]), float(vd[2])}; - } - - ///< return Chi2 at PCA candidate (no check for its validity) - float getChi2AtPCACandidate(int cand = 0) const { return mChi2[mOrder[cand]]; } - - ///< prepare copies of tracks at the V0 candidate (no check for the candidate validity) - /// must be called before getTrack(i,cand) query - bool propagateTracksToVertex(int cand = 0); - - ///< check if propagation of tracks to candidate vertex was done - bool isPropagateTracksToVertexDone(int cand = 0) const { return mTrPropDone[mOrder[cand]]; } - - ///< track param propagated to V0 candidate (no check for the candidate validity) - /// propagateTracksToVertex must be called in advance - Track& getTrack(int i, int cand = 0) - { - if (!mTrPropDone[mOrder[cand]]) { - throw std::runtime_error("propagateTracksToVertex was not called yet"); - } - return mCandTr[mOrder[cand]][i]; - } - - const Track& getTrack(int i, int cand = 0) const - { - if (!mTrPropDone[mOrder[cand]]) { - throw std::runtime_error("propagateTracksToVertex was not called yet"); - } - return mCandTr[mOrder[cand]][i]; - } - - ///< create parent track param with errors for decay vertex - o2::track::TrackParCov createParentTrackParCov(int cand = 0, bool sectorAlpha = true) const; - - ///< create parent track param w/o errors for decay vertex - o2::track::TrackPar createParentTrackPar(int cand = 0, bool sectorAlpha = true) const; - - ///< calculate on the fly track param (no cov mat) at candidate, check isValid to make sure propagation was successful - o2::track::TrackPar getTrackParamAtPCA(int i, int cand = 0) const; - - ///< recalculate PCA as a cov-matrix weighted mean, even if absDCA method was used - bool recalculatePCAWithErrors(int cand = 0); - - MatSym3D calcPCACovMatrix(int cand = 0) const; - - std::array calcPCACovMatrixFlat(int cand = 0) const - { - auto m = calcPCACovMatrix(cand); - return {float(m(0, 0)), float(m(1, 0)), float(m(1, 1)), float(m(2, 0)), float(m(2, 1)), float(m(2, 2))}; - } - - const Track* getOrigTrackPtr(int i) const { return mOrigTrPtr[i]; } - - ///< return number of iterations during minimization (no check for its validity) - int getNIterations(int cand = 0) const { return mNIters[mOrder[cand]]; } - void setPropagateToPCA(bool v = true) { mPropagateToPCA = v; } - void setMaxIter(int n = 20) { mMaxIter = n > 2 ? n : 2; } - void setMaxR(float r = 200.) { mMaxR2 = r * r; } - void setMaxDZIni(float d = 4.) { mMaxDZIni = d; } - void setMaxDXYIni(float d = 4.) { mMaxDXYIni = d > 0 ? d : 1e9; } - void setMaxChi2(float chi2 = 999.) { mMaxChi2 = chi2; } - void setBz(float bz) { mBz = std::abs(bz) > o2::constants::math::Almost0 ? bz : 0.f; } - void setMinParamChange(float x = 1e-3) { mMinParamChange = x > 1e-4 ? x : 1.e-4; } - void setMinRelChi2Change(float r = 0.9) { mMinRelChi2Change = r > 0.1 ? r : 999.; } - void setUseAbsDCA(bool v) { mUseAbsDCA = v; } - void setWeightedFinalPCA(bool v) { mWeightedFinalPCA = v; } - void setMaxDistance2ToMerge(float v) { mMaxDist2ToMergeSeeds = v; } - void setMatCorrType(o2::base::Propagator::MatCorrType m = o2::base::Propagator::MatCorrType::USEMatCorrLUT) { mMatCorr = m; } - void setUsePropagator(bool v) { mUsePropagator = v; } - void setRefitWithMatCorr(bool v) { mRefitWithMatCorr = v; } - void setMaxSnp(float s) { mMaxSnp = s; } - void setMaxStep(float s) { mMaxStep = s; } - void setMinXSeed(float x) { mMinXSeed = x; } - - int getNCandidates() const { return mCurHyp; } - int getMaxIter() const { return mMaxIter; } - float getMaxR() const { return std::sqrt(mMaxR2); } - float getMaxDZIni() const { return mMaxDZIni; } - float getMaxDXYIni() const { return mMaxDXYIni; } - float getMaxChi2() const { return mMaxChi2; } - float getMinParamChange() const { return mMinParamChange; } - float getBz() const { return mBz; } - float getMaxDistance2ToMerge() const { return mMaxDist2ToMergeSeeds; } - bool getUseAbsDCA() const { return mUseAbsDCA; } - bool getWeightedFinalPCA() const { return mWeightedFinalPCA; } - bool getPropagateToPCA() const { return mPropagateToPCA; } - o2::base::Propagator::MatCorrType getMatCorrType() const { return mMatCorr; } - bool getUsePropagator() const { return mUsePropagator; } - bool getRefitWithMatCorr() const { return mRefitWithMatCorr; } - float getMaxSnp() const { return mMaxSnp; } - float getMasStep() const { return mMaxStep; } - float getMinXSeed() const { return mMinXSeed; } - - template - int process(const Tr&... args); - void print() const; - - protected: - bool calcPCACoefs(); - bool calcInverseWeight(); - void calcResidDerivatives(); - void calcResidDerivativesNoErr(); - void calcRMatrices(); - void calcChi2Derivatives(); - void calcChi2DerivativesNoErr(); - void calcPCA(); - void calcPCANoErr(); - void calcTrackResiduals(); - void calcTrackDerivatives(); - double calcChi2() const; - double calcChi2NoErr() const; - bool correctTracks(const VecND& corrX); - bool minimizeChi2(); - bool minimizeChi2NoErr(); - bool roughDZCut() const; - bool closerToAlternative() const; - bool propagateToX(o2::track::TrackParCov& t, float x) const; - bool propagateParamToX(o2::track::TrackPar& t, float x) const; - - static double getAbsMax(const VecND& v); - ///< track param positions at V0 candidate (no check for the candidate validity) - const Vec3D& getTrackPos(int i, int cand = 0) const { return mTrPos[mOrder[cand]][i]; } - - ///< track X-param at V0 candidate (no check for the candidate validity) - float getTrackX(int i, int cand = 0) const { return getTrackPos(i, cand)[0]; } - - MatStd3D getTrackRotMatrix(int i) const // generate 3D matrix for track rotation to global frame - { - MatStd3D mat; - mat(2, 2) = 1; - mat(0, 0) = mat(1, 1) = mTrAux[i].c; - mat(0, 1) = -mTrAux[i].s; - mat(1, 0) = mTrAux[i].s; - return mat; - } - - MatSym3D getTrackCovMatrix(int i, int cand = 0) const // generate covariance matrix of track position, adding fake X error - { - const auto& trc = mCandTr[mOrder[cand]][i]; - MatSym3D mat; - mat(0, 0) = trc.getSigmaY2() * XerrFactor; - mat(1, 1) = trc.getSigmaY2(); - mat(2, 2) = trc.getSigmaZ2(); - mat(2, 1) = trc.getSigmaZY(); - return mat; - } - - void assign(int) {} - template - void assign(int i, const T& t, const Tr&... args) - { - static_assert(std::is_convertible(), "Wrong track type"); - mOrigTrPtr[i] = &t; - assign(i + 1, args...); - } - - void clear() - { - mCurHyp = 0; - mAllowAltPreference = true; - } - - static void setTrackPos(Vec3D& pnt, const Track& tr) - { - pnt[0] = tr.getX(); - pnt[1] = tr.getY(); - pnt[2] = tr.getZ(); - } - - private: - // vectors of 1st derivatives of track local residuals over X parameters - std::array, N> mDResidDx; - // vectors of 1nd derivatives of track local residuals over X parameters - // (cross-derivatives DR/(dx_j*dx_k) = 0 for j!=k, therefore the hessian is diagonal) - std::array, N> mD2ResidDx2; - VecND mDChi2Dx; // 1st derivatives of chi2 over tracks X params - MatSymND mD2Chi2Dx2; // 2nd derivatives of chi2 over tracks X params (symmetric matrix) - MatSymND mCosDif; // matrix with cos(alp_j-alp_i) for j mOrigTrPtr; - std::array mTrAux; // Aux track info for each track at each cand. vertex - CrossInfo mCrossings; // info on track crossing - - std::array mTrcEInv; // errors for each track at each cand. vertex - std::array mCandTr; // tracks at each cond. vertex (Note: Errors are at seed XY point) - std::array mTrCFVT; // TrackCoefVtx for each track at each cand. vertex - std::array mTrDer; // Track derivativse - std::array mTrPos; // Track positions - std::array mTrRes; // Track residuals - std::array mPCA; // PCA for each vertex candidate - std::array mChi2 = {0}; // Chi2 at PCA candidate - std::array mNIters; // number of iterations for each seed - std::array mTrPropDone; // Flag that the tracks are fully propagated to PCA - MatSym3D mWeightInv; // inverse weight of single track, [sum{M^T E M}]^-1 in EQ.T - std::array mOrder{0}; - int mCurHyp = 0; - int mCrossIDCur = 0; - int mCrossIDAlt = -1; - bool mAllowAltPreference = true; // if the fit converges to alternative PCA seed, abandon the current one - bool mUseAbsDCA = false; // use abs. distance minimization rather than chi2 - bool mWeightedFinalPCA = false; // recalculate PCA as a cov-matrix weighted mean, even if absDCA method was used - bool mPropagateToPCA = true; // create tracks version propagated to PCA - bool mUsePropagator = false; // use propagator with 3D B-field, set automatically if material correction is requested - bool mRefitWithMatCorr = false; // when doing propagateTracksToVertex, propagate tracks to V0 with material corrections and rerun minimization again - o2::base::Propagator::MatCorrType mMatCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; // material corrections type - int mMaxIter = 20; // max number of iterations - float mBz = 0; // bz field, to be set by user - float mMaxR2 = 200. * 200.; // reject PCA's above this radius - float mMinXSeed = -50.; // reject seed if it corresponds to X-param < mMinXSeed for one of candidates (e.g. X becomes strongly negative) - float mMaxDZIni = 4.; // reject (if>0) PCA candidate if tracks DZ exceeds threshold - float mMaxDXYIni = 4.; // reject (if>0) PCA candidate if tracks dXY exceeds threshold - float mMinParamChange = 1e-3; // stop iterations if largest change of any X is smaller than this - float mMinRelChi2Change = 0.9; // stop iterations is chi2/chi2old > this - float mMaxChi2 = 100; // abs cut on chi2 or abs distance - float mMaxDist2ToMergeSeeds = 1.; // merge 2 seeds to their average if their distance^2 is below the threshold - float mMaxSnp = 0.95; // Max snp for propagation with Propagator - float mMaxStep = 2.0; // Max step for propagation with Propagator - - ClassDefNV(DCAFitterN, 1); -}; - -///_________________________________________________________________________ -template -template -int DCAFitterN::process(const Tr&... args) -{ - // This is a main entry point: fit PCA of N tracks - static_assert(sizeof...(args) == N, "incorrect number of input tracks"); - assign(0, args...); - clear(); - for (int i = 0; i < N; i++) { - mTrAux[i].set(*mOrigTrPtr[i], mBz); - } - if (!mCrossings.set(mTrAux[0], *mOrigTrPtr[0], mTrAux[1], *mOrigTrPtr[1], mMaxDXYIni)) { // even for N>2 it should be enough to test just 1 loop - return 0; // no crossing - } - if (mUseAbsDCA) { - calcRMatrices(); // needed for fast residuals derivatives calculation in case of abs. distance minimization - } - if (mCrossings.nDCA == MAXHYP) { // if there are 2 candidates and they are too close, chose their mean as a starting point - auto dst2 = (mCrossings.xDCA[0] - mCrossings.xDCA[1]) * (mCrossings.xDCA[0] - mCrossings.xDCA[1]) + - (mCrossings.yDCA[0] - mCrossings.yDCA[1]) * (mCrossings.yDCA[0] - mCrossings.yDCA[1]); - if (dst2 < mMaxDist2ToMergeSeeds) { - mCrossings.nDCA = 1; - mCrossings.xDCA[0] = 0.5 * (mCrossings.xDCA[0] + mCrossings.xDCA[1]); - mCrossings.yDCA[0] = 0.5 * (mCrossings.yDCA[0] + mCrossings.yDCA[1]); - } - } - // check all crossings - for (int ic = 0; ic < mCrossings.nDCA; ic++) { - // check if radius is acceptable - if (mCrossings.xDCA[ic] * mCrossings.xDCA[ic] + mCrossings.yDCA[ic] * mCrossings.yDCA[ic] > mMaxR2) { - continue; - } - mCrossIDCur = ic; - mCrossIDAlt = (mCrossings.nDCA == 2 && mAllowAltPreference) ? 1 - ic : -1; // works for max 2 crossings - mNIters[mCurHyp] = 0; - mTrPropDone[mCurHyp] = false; - mChi2[mCurHyp] = -1.; - mPCA[mCurHyp][0] = mCrossings.xDCA[ic]; - mPCA[mCurHyp][1] = mCrossings.yDCA[ic]; - - if (mUseAbsDCA ? minimizeChi2NoErr() : minimizeChi2()) { - mOrder[mCurHyp] = mCurHyp; - if (mPropagateToPCA && !propagateTracksToVertex(mCurHyp)) { - continue; // discard candidate if failed to propagate to it - } - mCurHyp++; - } - } - - for (int i = mCurHyp; i--;) { // order in quality - for (int j = i; j--;) { - if (mChi2[mOrder[i]] < mChi2[mOrder[j]]) { - std::swap(mOrder[i], mOrder[j]); - } - } - } - if (mUseAbsDCA && mWeightedFinalPCA) { - for (int i = mCurHyp; i--;) { - recalculatePCAWithErrors(i); - } - } - - return mCurHyp; -} - -//__________________________________________________________________________ -template -bool DCAFitterN::calcPCACoefs() -{ - //< calculate Ti matrices for global vertex decomposition to V = sum_{0 -bool DCAFitterN::calcInverseWeight() -{ - //< calculate [sum_{0 -void DCAFitterN::calcResidDerivatives() -{ - //< calculate matrix of derivatives for weighted chi2: residual i vs parameter X of track j - MatStd3D matMT; - for (int i = N; i--;) { // residual being differentiated - const auto& taux = mTrAux[i]; - for (int j = N; j--;) { // track over which we differentiate - const auto& matT = mTrCFVT[mCurHyp][j]; // coefficient matrix for track J - const auto& trDx = mTrDer[mCurHyp][j]; // track point derivs over track X param - auto& dr1 = mDResidDx[i][j]; - auto& dr2 = mD2ResidDx2[i][j]; - // calculate M_i^tr * T_j - matMT[0][0] = taux.c * matT[0][0] + taux.s * matT[1][0]; - matMT[0][1] = taux.c * matT[0][1] + taux.s * matT[1][1]; - matMT[0][2] = taux.c * matT[0][2] + taux.s * matT[1][2]; - matMT[1][0] = -taux.s * matT[0][0] + taux.c * matT[1][0]; - matMT[1][1] = -taux.s * matT[0][1] + taux.c * matT[1][1]; - matMT[1][2] = -taux.s * matT[0][2] + taux.c * matT[1][2]; - matMT[2][0] = matT[2][0]; - matMT[2][1] = matT[2][1]; - matMT[2][2] = matT[2][2]; - - // calculate DResid_i/Dx_j = (delta_ij - M_i^tr * T_j) * DTrack_k/Dx_k - dr1[0] = -(matMT[0][0] + matMT[0][1] * trDx.dydx + matMT[0][2] * trDx.dzdx); - dr1[1] = -(matMT[1][0] + matMT[1][1] * trDx.dydx + matMT[1][2] * trDx.dzdx); - dr1[2] = -(matMT[2][0] + matMT[2][1] * trDx.dydx + matMT[2][2] * trDx.dzdx); - - // calculate D2Resid_I/(Dx_J Dx_K) = (delta_ijk - M_i^tr * T_j * delta_jk) * D2Track_k/dx_k^2 - dr2[0] = -(matMT[0][1] * trDx.d2ydx2 + matMT[0][2] * trDx.d2zdx2); - dr2[1] = -(matMT[1][1] * trDx.d2ydx2 + matMT[1][2] * trDx.d2zdx2); - dr2[2] = -(matMT[2][1] * trDx.d2ydx2 + matMT[2][2] * trDx.d2zdx2); - - if (i == j) { - dr1[0] += 1.; - dr1[1] += trDx.dydx; - dr1[2] += trDx.dzdx; - - dr2[1] += trDx.d2ydx2; - dr2[2] += trDx.d2zdx2; - } - } // track over which we differentiate - } // residual being differentiated -} - -//__________________________________________________________________________ -template -void DCAFitterN::calcResidDerivativesNoErr() -{ - //< calculate matrix of derivatives for absolute distance chi2: residual i vs parameter X of track j - constexpr double NInv1 = 1. - NInv; // profit from Rii = I/Ninv - for (int i = N; i--;) { // residual being differentiated - const auto& trDxi = mTrDer[mCurHyp][i]; // track point derivs over track X param - auto& dr1ii = mDResidDx[i][i]; - auto& dr2ii = mD2ResidDx2[i][i]; - dr1ii[0] = NInv1; - dr1ii[1] = NInv1 * trDxi.dydx; - dr1ii[2] = NInv1 * trDxi.dzdx; - - dr2ii[0] = 0; - dr2ii[1] = NInv1 * trDxi.d2ydx2; - dr2ii[2] = NInv1 * trDxi.d2zdx2; - - for (int j = i; j--;) { // track over which we differentiate - auto& dr1ij = mDResidDx[i][j]; - auto& dr1ji = mDResidDx[j][i]; - const auto& trDxj = mTrDer[mCurHyp][j]; // track point derivs over track X param - auto cij = mCosDif[i][j], sij = mSinDif[i][j]; // M_i^T*M_j / N matrices non-trivial elements = {ci*cj+si*sj , si*cj-ci*sj }, see 5 in ref. - - // calculate DResid_i/Dx_j = (delta_ij - R_ij) * DTrack_j/Dx_j for j -void DCAFitterN::calcRMatrices() -{ - //< calculate Rij = 1/N M_i^T * M_j matrices (rotation from j-th track to i-th track frame) - for (int i = N; i--;) { - const auto& mi = mTrAux[i]; - for (int j = i; j--;) { - const auto& mj = mTrAux[j]; - mCosDif[i][j] = (mi.c * mj.c + mi.s * mj.s) * NInv; // cos(alp_i-alp_j) / N - mSinDif[i][j] = (mi.s * mj.c - mi.c * mj.s) * NInv; // sin(alp_i-alp_j) / N - } - } -} - -//__________________________________________________________________________ -template -void DCAFitterN::calcChi2Derivatives() -{ - //< calculate 1st and 2nd derivatives of wighted DCA (chi2) over track parameters X, see EQ.Chi2 in the ref - std::array, N> covIDrDx; // tempory vectors of covI_j * dres_j/dx_i - - // chi2 1st derivative - for (int i = N; i--;) { - auto& dchi1 = mDChi2Dx[i]; // DChi2/Dx_i = sum_j { res_j * covI_j * Dres_j/Dx_i } - dchi1 = 0; - for (int j = N; j--;) { - const auto& res = mTrRes[mCurHyp][j]; // vector of residuals of track j - const auto& covI = mTrcEInv[mCurHyp][j]; // inverse cov matrix of track j - const auto& dr1 = mDResidDx[j][i]; // vector of j-th residuals 1st derivative over X param of track i - auto& cidr = covIDrDx[i][j]; // vector covI_j * dres_j/dx_i, save for 2nd derivative calculation - cidr[0] = covI.sxx * dr1[0]; - cidr[1] = covI.syy * dr1[1] + covI.syz * dr1[2]; - cidr[2] = covI.syz * dr1[1] + covI.szz * dr1[2]; - // calculate res_i * covI_j * dres_j/dx_i - dchi1 += ROOT::Math::Dot(res, cidr); - } - } - // chi2 2nd derivative - for (int i = N; i--;) { - for (int j = i + 1; j--;) { // symmetric matrix - auto& dchi2 = mD2Chi2Dx2[i][j]; // D2Chi2/Dx_i/Dx_j = sum_k { Dres_k/Dx_j * covI_k * Dres_k/Dx_i + res_k * covI_k * D2res_k/Dx_i/Dx_j } - dchi2 = 0; - for (int k = N; k--;) { - const auto& dr1j = mDResidDx[k][j]; // vector of k-th residuals 1st derivative over X param of track j - const auto& cidrkj = covIDrDx[i][k]; // vector covI_k * dres_k/dx_i - dchi2 += ROOT::Math::Dot(dr1j, cidrkj); - if (k == j) { - const auto& res = mTrRes[mCurHyp][k]; // vector of residuals of track k - const auto& covI = mTrcEInv[mCurHyp][k]; // inverse cov matrix of track k - const auto& dr2ij = mD2ResidDx2[k][j]; // vector of k-th residuals 2nd derivative over X params of track j - dchi2 += res[0] * covI.sxx * dr2ij[0] + res[1] * (covI.syy * dr2ij[1] + covI.syz * dr2ij[2]) + res[2] * (covI.syz * dr2ij[1] + covI.szz * dr2ij[2]); - } - } - } - } -} - -//__________________________________________________________________________ -template -void DCAFitterN::calcChi2DerivativesNoErr() -{ - //< calculate 1st and 2nd derivatives of abs DCA (chi2) over track parameters X, see (6) in the ref - for (int i = N; i--;) { - auto& dchi1 = mDChi2Dx[i]; // DChi2/Dx_i = sum_j { res_j * Dres_j/Dx_i } - dchi1 = 0; // chi2 1st derivative - for (int j = N; j--;) { - const auto& res = mTrRes[mCurHyp][j]; // vector of residuals of track j - const auto& dr1 = mDResidDx[j][i]; // vector of j-th residuals 1st derivative over X param of track i - dchi1 += ROOT::Math::Dot(res, dr1); - if (i >= j) { // symmetrix matrix - // chi2 2nd derivative - auto& dchi2 = mD2Chi2Dx2[i][j]; // D2Chi2/Dx_i/Dx_j = sum_k { Dres_k/Dx_j * covI_k * Dres_k/Dx_i + res_k * covI_k * D2res_k/Dx_i/Dx_j } - dchi2 = ROOT::Math::Dot(mTrRes[mCurHyp][i], mD2ResidDx2[i][j]); - for (int k = N; k--;) { - dchi2 += ROOT::Math::Dot(mDResidDx[k][i], mDResidDx[k][j]); - } - } - } - } -} - -//___________________________________________________________________ -template -void DCAFitterN::calcPCA() -{ - // calculate point of closest approach for N prongs - mPCA[mCurHyp] = mTrCFVT[mCurHyp][N - 1] * mTrPos[mCurHyp][N - 1]; - for (int i = N - 1; i--;) { - mPCA[mCurHyp] += mTrCFVT[mCurHyp][i] * mTrPos[mCurHyp][i]; - } -} - -//___________________________________________________________________ -template -bool DCAFitterN::recalculatePCAWithErrors(int cand) -{ - // recalculate PCA as a cov-matrix weighted mean, even if absDCA method was used - if (isPropagateTracksToVertexDone(cand) && !propagateTracksToVertex(cand)) { - return false; - } - int saveCurHyp = mCurHyp; - mCurHyp = mOrder[cand]; - if (mUseAbsDCA) { - for (int i = N; i--;) { - mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor); // prepare inverse cov.matrices at starting point - } - if (!calcPCACoefs()) { - mCurHyp = saveCurHyp; - return false; - } - } - auto oldPCA = mPCA[mOrder[cand]]; - calcPCA(); - mCurHyp = saveCurHyp; - return true; -} - -//___________________________________________________________________ -template -void DCAFitterN::calcPCANoErr() -{ - // calculate point of closest approach for N prongs w/o errors - auto& pca = mPCA[mCurHyp]; - o2::math_utils::rotateZd(mTrPos[mCurHyp][N - 1][0], mTrPos[mCurHyp][N - 1][1], pca[0], pca[1], mTrAux[N - 1].s, mTrAux[N - 1].c); - // RRRR mTrAux[N-1].loc2glo(mTrPos[mCurHyp][N-1][0], mTrPos[mCurHyp][N-1][1], pca[0], pca[1] ); - pca[2] = mTrPos[mCurHyp][N - 1][2]; - for (int i = N - 1; i--;) { - double x, y; - o2::math_utils::rotateZd(mTrPos[mCurHyp][i][0], mTrPos[mCurHyp][i][1], x, y, mTrAux[i].s, mTrAux[i].c); - // RRRR mTrAux[i].loc2glo(mTrPos[mCurHyp][i][0], mTrPos[mCurHyp][i][1], x, y ); - pca[0] += x; - pca[1] += y; - pca[2] += mTrPos[mCurHyp][i][2]; - } - pca[0] *= NInv; - pca[1] *= NInv; - pca[2] *= NInv; -} - -//___________________________________________________________________ -template -ROOT::Math::SMatrix> DCAFitterN::calcPCACovMatrix(int cand) const -{ - // calculate covariance matrix for the point of closest approach - MatSym3D covm; - int nAdded = 0; - for (int i = N; i--;) { // calculate sum of inverses - // MatSym3D covTr = ROOT::Math::Similarity(mUseAbsDCA ? getTrackRotMatrix(i) : mTrCFVT[mOrder[cand]][i], getTrackCovMatrix(i, cand)); - // RS by using Similarity(mTrCFVT[mOrder[cand]][i], getTrackCovMatrix(i, cand)) we underestimate the error, use simple rotation - MatSym3D covTr = ROOT::Math::Similarity(getTrackRotMatrix(i), getTrackCovMatrix(i, cand)); - if (covTr.Invert()) { - covm += covTr; - nAdded++; - } - } - if (nAdded && covm.Invert()) { - return covm; - } - // correct way has failed, use simple sum - MatSym3D covmSum; - for (int i = N; i--;) { - MatSym3D covTr = ROOT::Math::Similarity(getTrackRotMatrix(i), getTrackCovMatrix(i, cand)); - } - return covmSum; -} - -//___________________________________________________________________ -template -void DCAFitterN::calcTrackResiduals() -{ - // calculate residuals - Vec3D vtxLoc; - for (int i = N; i--;) { - mTrRes[mCurHyp][i] = mTrPos[mCurHyp][i]; - vtxLoc = mPCA[mCurHyp]; - o2::math_utils::rotateZInvd(vtxLoc[0], vtxLoc[1], vtxLoc[0], vtxLoc[1], mTrAux[i].s, mTrAux[i].c); // glo->loc - mTrRes[mCurHyp][i] -= vtxLoc; - } -} - -//___________________________________________________________________ -template -inline void DCAFitterN::calcTrackDerivatives() -{ - // calculate track derivatives over X param - for (int i = N; i--;) { - mTrDer[mCurHyp][i].set(mCandTr[mCurHyp][i], mBz); - } -} - -//___________________________________________________________________ -template -inline double DCAFitterN::calcChi2() const -{ - // calculate current chi2 - double chi2 = 0; - for (int i = N; i--;) { - const auto& res = mTrRes[mCurHyp][i]; - const auto& covI = mTrcEInv[mCurHyp][i]; - chi2 += res[0] * res[0] * covI.sxx + res[1] * res[1] * covI.syy + res[2] * res[2] * covI.szz + 2. * res[1] * res[2] * covI.syz; - } - return chi2; -} - -//___________________________________________________________________ -template -inline double DCAFitterN::calcChi2NoErr() const -{ - // calculate current chi2 of abs. distance minimization - double chi2 = 0; - for (int i = N; i--;) { - const auto& res = mTrRes[mCurHyp][i]; - chi2 += res[0] * res[0] + res[1] * res[1] + res[2] * res[2]; - } - return chi2; -} - -//___________________________________________________________________ -template -bool DCAFitterN::correctTracks(const VecND& corrX) -{ - // propagate tracks to updated X - for (int i = N; i--;) { - const auto& trDer = mTrDer[mCurHyp][i]; - auto dx2h = 0.5 * corrX[i] * corrX[i]; - mTrPos[mCurHyp][i][0] -= corrX[i]; - mTrPos[mCurHyp][i][1] -= trDer.dydx * corrX[i] - dx2h * trDer.d2ydx2; - mTrPos[mCurHyp][i][2] -= trDer.dzdx * corrX[i] - dx2h * trDer.d2zdx2; - } - return true; -} - -//___________________________________________________________________ -template -bool DCAFitterN::propagateTracksToVertex(int icand) -{ - // propagate tracks to current vertex - int ord = mOrder[icand]; - if (mTrPropDone[ord]) { - return true; - } - - // need to refit taking as a seed already found vertex - if (mRefitWithMatCorr) { - int curHypSav = mCurHyp, curCrosIDAlt = mCrossIDAlt; // save - mCurHyp = ord; - mCrossIDAlt = -1; // disable alternative check - auto restore = [this, curHypSav, curCrosIDAlt]() { this->mCurHyp = curHypSav; this->mCrossIDAlt = curCrosIDAlt; }; - if (!(mUseAbsDCA ? minimizeChi2NoErr() : minimizeChi2())) { // do final propagation - restore(); - return false; - } - restore(); - } - - for (int i = N; i--;) { - if (mUseAbsDCA || mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) { - mCandTr[ord][i] = *mOrigTrPtr[i]; // fetch the track again, as mCandTr might have been propagated w/o errors or material corrections might be wrong - } - auto x = mTrAux[i].c * mPCA[ord][0] + mTrAux[i].s * mPCA[ord][1]; // X of PCA in the track frame - if (!propagateToX(mCandTr[ord][i], x)) { - return false; - } - } - - mTrPropDone[ord] = true; - return true; -} - -//___________________________________________________________________ -template -inline o2::track::TrackPar DCAFitterN::getTrackParamAtPCA(int i, int icand) const -{ - // propagate tracks param only to current vertex (if not already done) - int ord = mOrder[icand]; - o2::track::TrackPar trc(mCandTr[ord][i]); - if (!mTrPropDone[ord]) { - auto x = mTrAux[i].c * mPCA[ord][0] + mTrAux[i].s * mPCA[ord][1]; // X of PCA in the track frame - if (!propagateParamToX(trc, x)) { - trc.invalidate(); - } - } - return std::move(trc); -} - -//___________________________________________________________________ -template -inline double DCAFitterN::getAbsMax(const VecND& v) -{ - double mx = -1; - for (int i = N; i--;) { - auto vai = std::abs(v[i]); - if (mx < vai) { - mx = vai; - } - } - return mx; -} - -//___________________________________________________________________ -template -bool DCAFitterN::minimizeChi2() -{ - // find best chi2 (weighted DCA) of N tracks in the vicinity of the seed PCA - for (int i = N; i--;) { - mCandTr[mCurHyp][i] = *mOrigTrPtr[i]; - auto x = mTrAux[i].c * mPCA[mCurHyp][0] + mTrAux[i].s * mPCA[mCurHyp][1]; // X of PCA in the track frame - if (x < mMinXSeed || !propagateToX(mCandTr[mCurHyp][i], x)) { - return false; - } - setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions - mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor); // prepare inverse cov.matrices at starting point - } - - if (mMaxDZIni > 0 && !roughDZCut()) { // apply rough cut on tracks Z difference - return false; - } - - if (!calcPCACoefs()) { // prepare tracks contribution matrices to the global PCA - return false; - } - calcPCA(); // current PCA - calcTrackResiduals(); // current track residuals - float chi2Upd, chi2 = calcChi2(); - do { - calcTrackDerivatives(); // current track derivatives (1st and 2nd) - calcResidDerivatives(); // current residals derivatives (1st and 2nd) - calcChi2Derivatives(); // current chi2 derivatives (1st and 2nd) - - // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1 - if (!mD2Chi2Dx2.Invert()) { - LOG(error) << "InversionFailed"; - return false; - } - VecND dx = mD2Chi2Dx2 * mDChi2Dx; - if (!correctTracks(dx)) { - return false; - } - calcPCA(); // updated PCA - if (mCrossIDAlt >= 0 && closerToAlternative()) { - mAllowAltPreference = false; - return false; - } - calcTrackResiduals(); // updated residuals - chi2Upd = calcChi2(); // updated chi2 - if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) { - chi2 = chi2Upd; - break; // converged - } - chi2 = chi2Upd; - } while (++mNIters[mCurHyp] < mMaxIter); - // - mChi2[mCurHyp] = chi2 * NInv; - return mChi2[mCurHyp] < mMaxChi2; -} - -//___________________________________________________________________ -template -bool DCAFitterN::minimizeChi2NoErr() -{ - // find best chi2 (absolute DCA) of N tracks in the vicinity of the PCA seed - - for (int i = N; i--;) { - mCandTr[mCurHyp][i] = *mOrigTrPtr[i]; - auto x = mTrAux[i].c * mPCA[mCurHyp][0] + mTrAux[i].s * mPCA[mCurHyp][1]; // X of PCA in the track frame - if (x < mMinXSeed || !propagateParamToX(mCandTr[mCurHyp][i], x)) { - return false; - } - setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions - } - if (mMaxDZIni > 0 && !roughDZCut()) { // apply rough cut on tracks Z difference - return false; - } - - calcPCANoErr(); // current PCA - calcTrackResiduals(); // current track residuals - float chi2Upd, chi2 = calcChi2NoErr(); - do { - calcTrackDerivatives(); // current track derivatives (1st and 2nd) - calcResidDerivativesNoErr(); // current residals derivatives (1st and 2nd) - calcChi2DerivativesNoErr(); // current chi2 derivatives (1st and 2nd) - - // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1 - if (!mD2Chi2Dx2.Invert()) { - LOG(error) << "InversionFailed"; - return false; - } - VecND dx = mD2Chi2Dx2 * mDChi2Dx; - if (!correctTracks(dx)) { - return false; - } - calcPCANoErr(); // updated PCA - if (mCrossIDAlt >= 0 && closerToAlternative()) { - mAllowAltPreference = false; - return false; - } - calcTrackResiduals(); // updated residuals - chi2Upd = calcChi2NoErr(); // updated chi2 - if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) { - chi2 = chi2Upd; - break; // converged - } - chi2 = chi2Upd; - } while (++mNIters[mCurHyp] < mMaxIter); - // - mChi2[mCurHyp] = chi2 * NInv; - return mChi2[mCurHyp] < mMaxChi2; -} - -//___________________________________________________________________ -template -bool DCAFitterN::roughDZCut() const -{ - // apply rough cut on DZ between the tracks in the seed point - bool accept = true; - for (int i = N; accept && i--;) { - for (int j = i; j--;) { - if (std::abs(mCandTr[mCurHyp][i].getZ() - mCandTr[mCurHyp][j].getZ()) > mMaxDZIni) { - accept = false; - break; - } - } - } - return accept; -} - -//___________________________________________________________________ -template -bool DCAFitterN::closerToAlternative() const -{ - // check if the point current PCA point is closer to the seeding XY point being tested or to alternative see (if any) - auto dxCur = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDCur], dyCur = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDCur]; - auto dxAlt = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDAlt], dyAlt = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDAlt]; - return dxCur * dxCur + dyCur * dyCur > dxAlt * dxAlt + dyAlt * dyAlt; -} - -//___________________________________________________________________ -template -void DCAFitterN::print() const -{ - LOG(info) << N << "-prong vertex fitter in " << (mUseAbsDCA ? "abs." : "weighted") << " distance minimization mode"; - LOG(info) << "Bz: " << mBz << " MaxIter: " << mMaxIter << " MaxChi2: " << mMaxChi2; - LOG(info) << "Stopping condition: Max.param change < " << mMinParamChange << " Rel.Chi2 change > " << mMinRelChi2Change; - LOG(info) << "Discard candidates for : Rvtx > " << getMaxR() << " DZ between tracks > " << mMaxDZIni; -} - -//___________________________________________________________________ -template -o2::track::TrackParCov DCAFitterN::createParentTrackParCov(int cand, bool sectorAlpha) const -{ - const auto& trP = getTrack(0, cand); - const auto& trN = getTrack(1, cand); - std::array covV = {0.}; - std::array pvecV = {0.}; - int q = 0; - for (int it = 0; it < N; it++) { - const auto& trc = getTrack(it, cand); - std::array pvecT = {0.}; - std::array covT = {0.}; - trc.getPxPyPzGlo(pvecT); - trc.getCovXYZPxPyPzGlo(covT); - constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component - for (int i = 0; i < 6; i++) { - covV[MomInd[i]] += covT[MomInd[i]]; - } - for (int i = 0; i < 3; i++) { - pvecV[i] += pvecT[i]; - } - q += trc.getCharge(); - } - auto covVtxV = calcPCACovMatrix(cand); - covV[0] = covVtxV(0, 0); - covV[1] = covVtxV(1, 0); - covV[2] = covVtxV(1, 1); - covV[3] = covVtxV(2, 0); - covV[4] = covVtxV(2, 1); - covV[5] = covVtxV(2, 2); - return std::move(o2::track::TrackParCov(getPCACandidatePos(cand), pvecV, covV, q, sectorAlpha)); -} - -//___________________________________________________________________ -template -o2::track::TrackPar DCAFitterN::createParentTrackPar(int cand, bool sectorAlpha) const -{ - const auto& trP = getTrack(0, cand); - const auto& trN = getTrack(1, cand); - const auto& wvtx = getPCACandidate(cand); - std::array pvecV = {0.}; - int q = 0; - for (int it = 0; it < N; it++) { - const auto& trc = getTrack(it, cand); - std::array pvecT = {0.}; - trc.getPxPyPzGlo(pvecT); - for (int i = 0; i < 3; i++) { - pvecV[i] += pvecT[i]; - } - q += trc.getCharge(); - } - const std::array vertex = {(float)wvtx[0], (float)wvtx[1], (float)wvtx[2]}; - return std::move(o2::track::TrackPar(vertex, pvecV, q, sectorAlpha)); -} - -//___________________________________________________________________ -template -inline bool DCAFitterN::propagateParamToX(o2::track::TrackPar& t, float x) const -{ - if (mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) { - return o2::base::Propagator::Instance()->PropagateToXBxByBz(t, x, mMaxSnp, mMaxStep, mMatCorr); - } else { - return t.propagateParamTo(x, mBz); - } -} - -//___________________________________________________________________ -template -inline bool DCAFitterN::propagateToX(o2::track::TrackParCov& t, float x) const -{ - if (mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) { - return o2::base::Propagator::Instance()->PropagateToXBxByBz(t, x, mMaxSnp, mMaxStep, mMatCorr); - } else { - return t.propagateTo(x, mBz); - } -} - -using DCAFitter2 = DCAFitterN<2, o2::track::TrackParCov>; -using DCAFitter3 = DCAFitterN<3, o2::track::TrackParCov>; - -} // namespace vertexing -} // namespace o2 -#endif // _ALICEO2_DCA_FITTERN_ diff --git a/Common/DCAFitter/include/DCAFitter/FwdDCAFitterN.h b/Common/DCAFitter/include/DCAFitter/FwdDCAFitterN.h deleted file mode 100644 index 7bcb7966d3fe1..0000000000000 --- a/Common/DCAFitter/include/DCAFitter/FwdDCAFitterN.h +++ /dev/null @@ -1,1265 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file FwdDCAFitterN.h -/// \brief Defintions for N-prongs secondary vertex fit -/// \author ruben.shahoyan@cern.ch, adapted from central barrel to fwd rapidities by Rita Sadek, rita.sadek@cern.ch -/// For the formulae derivation see /afs/cern.ch/user/s/shahoian/public/O2/DCAFitter/DCAFitterN.pdf - -#ifndef _ALICEO2_DCA_FWDFITTERN_ -#define _ALICEO2_DCA_FWDFITTERN_ -#include -#include "MathUtils/Cartesian.h" -#include "ReconstructionDataFormats/TrackFwd.h" -#include "ReconstructionDataFormats/Track.h" -#include "DCAFitter/HelixHelper.h" -#include - -namespace o2 -{ -namespace vertexing -{ - -///__________________________________________________________________________________ -///< Fwd Inverse cov matrix (augmented by a dummy Z error) of the point defined by the track -struct FwdTrackCovI { - float sxx, syy, sxy, szz; - - FwdTrackCovI(const o2::track::TrackParCovFwd& trc, float zerrFactor = 1.) { set(trc, zerrFactor); } - FwdTrackCovI() = default; - void set(const o2::track::TrackParCovFwd& trc, float zerrFactor = 1) - { - float cxx = trc.getSigma2X(), cyy = trc.getSigma2Y(), cxy = trc.getSigmaXY(), czz = cyy * zerrFactor; - float detXY = cxx * cyy - cxy * cxy; - if (detXY > 0.) { - auto detXYI = 1. / detXY; - sxx = cyy * detXYI; - syy = cxx * detXYI; - sxy = -cxy * detXYI; - szz = 1. / czz; - } else { - throw std::runtime_error("invalid track covariance"); - } - } -}; - -///__________________________________________________________________________ -///< Fwd derivative (up to 2) of the TrackParam position over its running param Z -struct FwdTrackDeriv { - float dxdz, dydz, d2xdz2, d2ydz2; - FwdTrackDeriv() = default; - FwdTrackDeriv(const o2::track::TrackParFwd& trc, float bz) { set(trc, bz); } - void set(const o2::track::TrackParFwd& trc, float bz) - { - float snp = trc.getSnp(), csp = std::sqrt((1. - snp) * (1. + snp)), cspI = 1. / csp, crv2c = trc.getCurvature(bz), tgl = trc.getTanl(), tglI = 1. / tgl; - if (crv2c == 0.) { - crv2c = (trc.getCharge()) * 0.3 * bz * (-1e-3); - } - - dxdz = csp * tglI; - dydz = snp * tglI; - d2xdz2 = crv2c * snp * tglI * tglI; - d2ydz2 = -crv2c * csp * tglI * tglI; - } -}; - -template -class FwdDCAFitterN -{ - static constexpr double NMin = 2; - static constexpr double NMax = 4; - static constexpr double NInv = 1. / N; - static constexpr int MAXHYP = 2; - static constexpr float ZerrFactor = 5.; // factor for conversion of track covXX to dummy covZZ - using Track = o2::track::TrackParCovFwd; - using TrackAuxPar = o2::track::TrackAuxPar; - using CrossInfo = o2::track::CrossInfo; - using Vec3D = ROOT::Math::SVector; - using VecND = ROOT::Math::SVector; - using MatSym3D = ROOT::Math::SMatrix>; - using MatStd3D = ROOT::Math::SMatrix>; - using MatSymND = ROOT::Math::SMatrix>; - using MatStdND = ROOT::Math::SMatrix>; - using SMatrix55 = ROOT::Math::SMatrix>; - using TrackCoefVtx = MatStd3D; - using ArrTrack = std::array; // container for prongs (tracks) at single vertex cand. - using ArrTrackCovI = std::array; // container for inv.cov.matrices at single vertex cand. - using ArrTrCoef = std::array; // container of TrackCoefVtx coefficients at single vertex cand. - using ArrTrDer = std::array; // container of Track 1st and 2nd derivative over their Z param - using ArrTrPos = std::array; // container of Track positions - - public: - static constexpr int getNProngs() { return N; } - - FwdDCAFitterN() = default; - FwdDCAFitterN(float bz, bool useAbsDCA, bool prop2DCA) : mBz(bz), mUseAbsDCA(useAbsDCA), mPropagateToPCA(prop2DCA) - { - static_assert(N >= NMin && N <= NMax, "N prongs outside of allowed range"); - } - - //========================================================================= - ///< return PCA candidate, by default best on is provided (no check for the index validity) - const Vec3D& getPCACandidate(int cand = 0) const { return mPCA[mOrder[cand]]; } - const auto getPCACandidatePos(int cand = 0) const - { - const auto& vd = mPCA[mOrder[cand]]; - return std::array{float(vd[0]), float(vd[1]), float(vd[2])}; - } - - ///< return Chi2 at PCA candidate (no check for its validity) - float getChi2AtPCACandidate(int cand = 0) const { return mChi2[mOrder[cand]]; } - - ///< prepare copies of tracks at the V0 candidate (no check for the candidate validity) - /// must be called before getTrack(i,cand) query - bool FwdpropagateTracksToVertex(int cand = 0); - - ///< check if propagation of tracks to candidate vertex was done - bool isPropagateTracksToVertexDone(int cand = 0) const { return mTrPropDone[mOrder[cand]]; } - - ///< track param propagated to V0 candidate (no check for the candidate validity) - /// propagateTracksToVertex must be called in advance - Track& getTrack(int i, int cand = 0) - { - if (!mTrPropDone[mOrder[cand]]) { - throw std::runtime_error("propagateTracksToVertex was not called yet"); - } - return mCandTr[mOrder[cand]][i]; - } - - ///< calculate on the fly track param (no cov mat) at candidate - o2::track::TrackParFwd FwdgetTrackParamAtPCA(int i, int cand = 0) const; - - MatSym3D calcPCACovMatrix(int cand = 0) const; - - std::array calcPCACovMatrixFlat(int cand = 0) const - { - auto m = calcPCACovMatrix(cand); - return {float(m(0, 0)), float(m(1, 0)), float(m(1, 1)), float(m(2, 0)), float(m(2, 1)), float(m(2, 2))}; - } - - const Track* getOrigTrackPtr(int i) const { return mOrigTrPtr[i]; } - - ///< return number of iterations during minimization (no check for its validity) - int getNIterations(int cand = 0) const { return mNIters[mOrder[cand]]; } - void setPropagateToPCA(bool v = true) { mPropagateToPCA = v; } - void setMaxIter(int n = 60) { mMaxIter = n > 2 ? n : 2; } - void setMaxR(float r = 200.) { mMaxR2 = r * r; } - void setMaxDXIni(float d = 4.) { mMaxDXIni = d; } - void setMaxChi2(float chi2 = 999.) { mMaxChi2 = chi2; } - void setBz(float bz) { mBz = std::abs(bz) > o2::constants::math::Almost0 ? bz : 0.f; } - void setMinParamChange(float x = 1e-3) { mMinParamChange = x > 1e-4 ? x : 1.e-4; } - void setMinRelChi2Change(float r = 0.9) { mMinRelChi2Change = r > 0.1 ? r : 999.; } - void setUseAbsDCA(bool v) { mUseAbsDCA = v; } - void setMaxDistance2ToMerge(float v) { mMaxDist2ToMergeSeeds = v; } - - int getNCandidates() const { return mCurHyp; } - int getMaxIter() const { return mMaxIter; } - float getMaxR() const { return std::sqrt(mMaxR2); } - float getMaxDXIni() const { return mMaxDXIni; } - float getMaxChi2() const { return mMaxChi2; } - float getMinParamChange() const { return mMinParamChange; } - float getBz() const { return mBz; } - double getK(double b) const { return std::abs(o2::constants::math::B2C * b); } - double getHz(double b) const { return std::copysign(1, b); } - - float getMaxDistance2ToMerge() const { return mMaxDist2ToMergeSeeds; } - bool getUseAbsDCA() const { return mUseAbsDCA; } - bool getPropagateToPCA() const { return mPropagateToPCA; } - - template - int process(const Tr&... args); - void print() const; - - protected: - bool FwdcalcPCACoefs(); - bool FwdcalcInverseWeight(); - void FwdcalcResidDerivatives(); - void FwdcalcResidDerivativesNoErr(); - void FwdcalcChi2Derivatives(); - void FwdcalcChi2DerivativesNoErr(); - void FwdcalcPCA(); - void FwdcalcPCANoErr(); - void FwdcalcTrackResiduals(); - void calcTrackDerivatives(); - float findZatXY(int cand = 0); - void findZatXY_mid(int cand = 0); - void findZatXY_lineApprox(int cand = 0); - void findZatXY_quad(int cand = 0); - void findZatXY_linear(int cand = 0); - double FwdcalcChi2() const; - double FwdcalcChi2NoErr() const; - bool FwdcorrectTracks(const VecND& corrZ); - bool minimizeChi2(); - bool minimizeChi2NoErr(); - bool roughDXCut() const; - bool closerToAlternative() const; - static double getAbsMax(const VecND& v); - - ///< track param positions at V0 candidate (no check for the candidate validity) - const Vec3D& getTrackPos(int i, int cand = 0) const { return mTrPos[mOrder[cand]][i]; } - - ///< track Z-param at V0 candidate (no check for the candidate validity) - float getTrackZ(int i, int cand = 0) const { return getTrackPos(i, cand)[2]; } - - MatStd3D getTrackRotMatrix(int i) const // generate 3D matrix for track rotation to global frame - // no rotation for fwd: mat=I - { - MatStd3D mat; - mat(0, 0) = 1; - mat(1, 1) = 1; - mat(2, 2) = 1; - return mat; - } - - MatSym3D getTrackCovMatrix(int i, int cand = 0) const // generate covariance matrix of track position, adding fake Z error - { - const auto& trc = mCandTr[mOrder[cand]][i]; - MatSym3D mat; - mat(0, 0) = trc.getSigma2X(); - mat(1, 1) = trc.getSigma2Y(); - mat(1, 0) = trc.getSigmaXY(); - mat(2, 2) = trc.getSigma2Y() * ZerrFactor; - return mat; - } - - void assign(int) {} - template - void assign(int i, const T& t, const Tr&... args) - { - static_assert(std::is_convertible(), "Wrong track type"); - mOrigTrPtr[i] = &t; - assign(i + 1, args...); - } - - void clear() - { - mCurHyp = 0; - mAllowAltPreference = true; - } - - static void setTrackPos(Vec3D& pnt, const Track& tr) - { - pnt[0] = tr.getX(); - pnt[1] = tr.getY(); - pnt[2] = tr.getZ(); - } - - private: - // vectors of 1st derivatives of track local residuals over Z parameters - std::array, N> mDResidDz; - // vectors of 1nd derivatives of track local residuals over Z parameters - std::array, N> mD2ResidDz2; - VecND mDChi2Dz; // 1st derivatives of chi2 over tracks Z params - MatSymND mD2Chi2Dz2; // 2nd derivatives of chi2 over tracks Z params (symmetric matrix) - - std::array mOrigTrPtr; - std::array mTrAux; // Aux track info for each track at each cand. vertex - CrossInfo mCrossings; // info on track crossing - - std::array mTrcEInv; // errors for each track at each cand. vertex - std::array mCandTr; // tracks at each cond. vertex (Note: Errors are at seed XY point) - std::array mTrCFVT; // TrackCoefVtx for each track at each cand. vertex - std::array mTrDer; // Track derivativse - std::array mTrPos; // Track positions - std::array mTrRes; // Track residuals - std::array mPCA; // PCA for each vertex candidate - std::array mChi2 = {0}; // Chi2 at PCA candidate - std::array mNIters; // number of iterations for each seed - std::array mTrPropDone; // Flag that the tracks are fully propagated to PCA - MatSym3D mWeightInv; // inverse weight of single track, [sum{M^T E M}]^-1 in EQ.T - std::array mOrder{0}; - int mCurHyp = 0; - int mCrossIDCur = 0; - int mCrossIDAlt = -1; - bool mAllowAltPreference = true; // if the fit converges to alternative PCA seed, abandon the current one - bool mUseAbsDCA = false; // use abs. distance minimization rather than chi2 - bool mPropagateToPCA = true; // create tracks version propagated to PCA - int mMaxIter = 60; // max number of iterations - float mBz = 0; // bz field, to be set by user - float mMaxR2 = 200. * 200.; // reject PCA's above this radius - float mMaxDXIni = 4.; // reject (if>0) PCA candidate if tracks DZ exceeds threshold - float mMinParamChange = 1e-5; // stop iterations if largest change of any X is smaller than this - float mMinRelChi2Change = 0.98; // stop iterations is chi2/chi2old > this - float mMaxChi2 = 100; // abs cut on chi2 or abs distance - float mMaxDist2ToMergeSeeds = 1.; // merge 2 seeds to their average if their distance^2 is below the threshold - - ClassDefNV(FwdDCAFitterN, 1); -}; - -///_________________________________________________________________________ -template -template -int FwdDCAFitterN::process(const Tr&... args) -{ - - static_assert(sizeof...(args) == N, "incorrect number of input tracks"); - assign(0, args...); - clear(); - - for (int i = 0; i < N; i++) { - mTrAux[i].set(*mOrigTrPtr[i], mBz); - } - - if (!mCrossings.set(mTrAux[0], *mOrigTrPtr[0], mTrAux[1], *mOrigTrPtr[1])) { // even for N>2 it should be enough to test just 1 loop - return 0; // no crossing - } - - if (mCrossings.nDCA == MAXHYP) { // if there are 2 candidates - auto dst2 = (mCrossings.xDCA[0] - mCrossings.xDCA[1]) * (mCrossings.xDCA[0] - mCrossings.xDCA[1]) + - (mCrossings.yDCA[0] - mCrossings.yDCA[1]) * (mCrossings.yDCA[0] - mCrossings.yDCA[1]); - - if (dst2 < mMaxDist2ToMergeSeeds) { - mCrossings.nDCA = 1; - mCrossings.xDCA[0] = 0.5 * (mCrossings.xDCA[0] + mCrossings.xDCA[1]); - mCrossings.yDCA[0] = 0.5 * (mCrossings.yDCA[0] + mCrossings.yDCA[1]); - } - } - - // check all crossings - for (int ic = 0; ic < mCrossings.nDCA; ic++) { - // check if radius is acceptable - if (mCrossings.xDCA[ic] * mCrossings.xDCA[ic] + mCrossings.yDCA[ic] * mCrossings.yDCA[ic] > mMaxR2) { - continue; - } - - mCrossIDCur = ic; - mCrossIDAlt = (mCrossings.nDCA == 2 && mAllowAltPreference) ? 1 - ic : -1; // works for max 2 crossings - mNIters[mCurHyp] = 0; - mTrPropDone[mCurHyp] = false; - mChi2[mCurHyp] = -1.; - - findZatXY_mid(mCurHyp); - - if (mUseAbsDCA ? minimizeChi2NoErr() : minimizeChi2()) { - mOrder[mCurHyp] = mCurHyp; - if (mPropagateToPCA && !FwdpropagateTracksToVertex(mCurHyp)) { - continue; - } - mCurHyp++; - } - } - - for (int i = mCurHyp; i--;) { // order in quality - for (int j = i; j--;) { - if (mChi2[mOrder[i]] < mChi2[mOrder[j]]) { - std::swap(mOrder[i], mOrder[j]); - } - } - } - - return mCurHyp; -} - -//__________________________________________________________________________ -template -bool FwdDCAFitterN::FwdcalcPCACoefs() -{ - //< calculate Ti matrices for global vertex decomposition to V = sum_{0 -bool FwdDCAFitterN::FwdcalcInverseWeight() -{ - //< calculate [sum_{0 -void FwdDCAFitterN::FwdcalcResidDerivatives() -{ - //< calculate matrix of derivatives for weighted chi2: residual i vs parameter Z of track j - MatStd3D matMT; - for (int i = N; i--;) { // residual being differentiated - // const auto& taux = mTrAux[i]; - for (int j = N; j--;) { // track over which we differentiate - const auto& matT = mTrCFVT[mCurHyp][j]; // coefficient matrix for track J - const auto& trDz = mTrDer[mCurHyp][j]; // track point derivs over track Z param - auto& dr1 = mDResidDz[i][j]; - auto& dr2 = mD2ResidDz2[i][j]; - // calculate M_i^transverse * T_j , M_i^transverse=I -> MT=T - matMT[0][0] = matT[0][0]; - matMT[0][1] = matT[0][1]; - matMT[0][2] = matT[0][2]; - matMT[1][0] = matT[1][0]; - matMT[1][1] = matT[1][1]; - matMT[1][2] = matT[1][2]; - matMT[2][0] = matT[2][0]; - matMT[2][1] = matT[2][1]; - matMT[2][2] = matT[2][2]; - - // calculate DResid_i/Dz_j = (delta_ij - M_i^tr * T_j) * DTrack_k/Dz_k - dr1[0] = -(matMT[0][0] * trDz.dxdz + matMT[0][1] * trDz.dydz + matMT[0][2]); - dr1[1] = -(matMT[1][0] * trDz.dxdz + matMT[1][1] * trDz.dydz + matMT[1][2]); - dr1[2] = -(matMT[2][0] * trDz.dxdz + matMT[2][1] * trDz.dydz + matMT[2][2]); - - // calculate D2Resid_I/(Dz_J Dz_K) = (delta_ijk - M_i^tr * T_j * delta_jk) * D2Track_k/dz_k^2 - dr2[0] = -(matMT[0][1] * trDz.d2ydz2 + matMT[0][0] * trDz.d2xdz2); - dr2[1] = -(matMT[1][1] * trDz.d2ydz2 + matMT[1][0] * trDz.d2xdz2); - dr2[2] = -(matMT[2][1] * trDz.d2ydz2 + matMT[2][0] * trDz.d2xdz2); - - if (i == j) { - dr1[0] += trDz.dxdz; - dr1[1] += trDz.dydz; - dr1[2] += 1.; - - dr2[0] += trDz.d2xdz2; - dr2[1] += trDz.d2ydz2; - } - } // track over which we differentiate - } // residual being differentiated -} - -//__________________________________________________________________________ -template -void FwdDCAFitterN::FwdcalcResidDerivativesNoErr() -{ - //< calculate matrix of derivatives for absolute distance chi2: residual i vs parameter Z of track j - constexpr double NInv1 = 1. - NInv; // profit from Rii = I/Ninv - for (int i = N; i--;) { // residual being differentiated - const auto& trDzi = mTrDer[mCurHyp][i]; // track point derivs over track Z param - auto& dr1ii = mDResidDz[i][i]; - auto& dr2ii = mD2ResidDz2[i][i]; - - dr1ii[0] = NInv1 * trDzi.dxdz; - dr1ii[1] = NInv1 * trDzi.dydz; - dr1ii[2] = NInv1; - - dr2ii[0] = NInv1 * trDzi.d2xdz2; - dr2ii[1] = NInv1 * trDzi.d2ydz2; - dr2ii[2] = 0; - - for (int j = i; j--;) { // track over which we differentiate - auto& dr1ij = mDResidDz[i][j]; - auto& dr1ji = mDResidDz[j][i]; - const auto& trDzj = mTrDer[mCurHyp][j]; // track point derivs over track Z param - - // calculate DResid_i/Dz_j = (delta_ij - R_ij) * DTrack_j/Dz_j for j -void FwdDCAFitterN::FwdcalcChi2Derivatives() -{ - //< calculate 1st and 2nd derivatives of wighted DCA (chi2) over track parameters Z - std::array, N> covIDrDz; // tempory vectors of covI_j * dres_j/dz_i - - // chi2 1st derivative - for (int i = N; i--;) { - auto& dchi1 = mDChi2Dz[i]; // DChi2/Dz_i = sum_j { res_j * covI_j * Dres_j/Dz_i } - dchi1 = 0; - for (int j = N; j--;) { - const auto& res = mTrRes[mCurHyp][j]; // vector of residuals of track j - const auto& covI = mTrcEInv[mCurHyp][j]; // inverse cov matrix of track j - const auto& dr1 = mDResidDz[j][i]; // vector of j-th residuals 1st derivative over Z param of track i - auto& cidr = covIDrDz[i][j]; // vector covI_j * dres_j/dz_i, save for 2nd derivative calculation - cidr[0] = covI.sxx * dr1[0] + covI.sxy * dr1[1]; - cidr[1] = covI.sxy * dr1[0] + covI.syy * dr1[1]; - cidr[2] = covI.szz * dr1[2]; - - dchi1 += ROOT::Math::Dot(res, cidr); - } - } - - // chi2 2nd derivative - for (int i = N; i--;) { - for (int j = i + 1; j--;) { // symmetric matrix - auto& dchi2 = mD2Chi2Dz2[i][j]; // D2Chi2/Dz_i/Dz_j = sum_k { Dres_k/Dz_j * covI_k * Dres_k/Dz_i + res_k * covI_k * D2res_k/Dz_i/Dz_j } - dchi2 = 0; - for (int k = N; k--;) { - const auto& dr1j = mDResidDz[k][j]; // vector of k-th residuals 1st derivative over Z param of track j - const auto& cidrkj = covIDrDz[i][k]; // vector covI_k * dres_k/dz_i - dchi2 += ROOT::Math::Dot(dr1j, cidrkj); - if (k == j) { - const auto& res = mTrRes[mCurHyp][k]; // vector of residuals of track k - const auto& covI = mTrcEInv[mCurHyp][k]; // inverse cov matrix of track k - const auto& dr2ij = mD2ResidDz2[k][j]; // vector of k-th residuals 2nd derivative over Z params of track j - dchi2 += res[0] * (covI.sxx * dr2ij[0] + covI.sxy * dr2ij[1]) + res[1] * (covI.sxy * dr2ij[0] + covI.syy * dr2ij[1]) + res[2] * covI.szz * dr2ij[2]; - } - } - } - } -} - -//__________________________________________________________________________ -template -void FwdDCAFitterN::FwdcalcChi2DerivativesNoErr() -{ - //< calculate 1st and 2nd derivatives of abs DCA (chi2) over track parameters Z - for (int i = N; i--;) { - auto& dchi1 = mDChi2Dz[i]; // DChi2/Dz_i = sum_j { res_j * Dres_j/Dz_i } - dchi1 = 0; // chi2 1st derivative - for (int j = N; j--;) { - const auto& res = mTrRes[mCurHyp][j]; // vector of residuals of track j - const auto& dr1 = mDResidDz[j][i]; // vector of j-th residuals 1st derivative over Z param of track i - dchi1 += ROOT::Math::Dot(res, dr1); - if (i >= j) { // symmetrix matrix - // chi2 2nd derivative - auto& dchi2 = mD2Chi2Dz2[i][j]; // D2Chi2/Dz_i/Dz_j = sum_k { Dres_k/Dz_j * covI_k * Dres_k/Dz_i + res_k * covI_k * D2res_k/Dz_i/Dz_j } - dchi2 = ROOT::Math::Dot(mTrRes[mCurHyp][i], mD2ResidDz2[i][j]); - for (int k = N; k--;) { - dchi2 += ROOT::Math::Dot(mDResidDz[k][i], mDResidDz[k][j]); - } - } - } - } -} - -//___________________________________________________________________ -template -void FwdDCAFitterN::FwdcalcPCA() -{ - // calculate point of closest approach for N prongs - // calculating V = sum (Ti*Pi) - mPCA[mCurHyp] = mTrCFVT[mCurHyp][N - 1] * mTrPos[mCurHyp][N - 1]; - for (int i = N - 1; i--;) { - mPCA[mCurHyp] += mTrCFVT[mCurHyp][i] * mTrPos[mCurHyp][i]; - } -} - -//___________________________________________________________________ -template -void FwdDCAFitterN::FwdcalcPCANoErr() -{ - // calculate point of closest approach for N prongs w/o errors - auto& pca = mPCA[mCurHyp]; - - pca[0] = mTrPos[mCurHyp][N - 1][0]; - pca[1] = mTrPos[mCurHyp][N - 1][1]; - pca[2] = mTrPos[mCurHyp][N - 1][2]; - - for (int i = N - 1; i--;) { - pca[0] += mTrPos[mCurHyp][i][0]; - pca[1] += mTrPos[mCurHyp][i][1]; - pca[2] += mTrPos[mCurHyp][i][2]; - } - pca[0] *= NInv; - pca[1] *= NInv; - pca[2] *= NInv; -} - -//___________________________________________________________________ -template -ROOT::Math::SMatrix> FwdDCAFitterN::calcPCACovMatrix(int cand) const -{ - // calculate covariance matrix for the point of closest approach - MatSym3D covm; - for (int i = N; i--;) { - covm += ROOT::Math::Similarity(mUseAbsDCA ? getTrackRotMatrix(i) : mTrCFVT[mOrder[cand]][i], getTrackCovMatrix(i, cand)); - } - return covm; -} - -//___________________________________________________________________ -template -void FwdDCAFitterN::FwdcalcTrackResiduals() -{ - // calculate residuals, res = Pi - V - Vec3D vtxLoc; - for (int i = N; i--;) { - mTrRes[mCurHyp][i] = mTrPos[mCurHyp][i]; - vtxLoc = mPCA[mCurHyp]; - mTrRes[mCurHyp][i] -= vtxLoc; - } -} - -//___________________________________________________________________ -template -inline void FwdDCAFitterN::calcTrackDerivatives() -{ - // calculate track derivatives over Z param - for (int i = N; i--;) { - mTrDer[mCurHyp][i].set(mCandTr[mCurHyp][i], mBz); - } -} - -//___________________________________________________________________ -template -inline double FwdDCAFitterN::FwdcalcChi2() const -{ - // calculate current chi2 - double chi2 = 0; - for (int i = N; i--;) { - const auto& res = mTrRes[mCurHyp][i]; - const auto& covI = mTrcEInv[mCurHyp][i]; - chi2 += res[0] * res[0] * covI.sxx + res[1] * res[1] * covI.syy + res[2] * res[2] * covI.szz + 2. * res[0] * res[1] * covI.sxy; - } - return chi2; -} - -//___________________________________________________________________ -template -inline double FwdDCAFitterN::FwdcalcChi2NoErr() const -{ - // calculate current chi2 of abs. distance minimization - double chi2 = 0; - for (int i = N; i--;) { - const auto& res = mTrRes[mCurHyp][i]; - chi2 += res[0] * res[0] + res[1] * res[1] + res[2] * res[2]; - } - return chi2; -} - -//___________________________________________________________________ -template -bool FwdDCAFitterN::FwdcorrectTracks(const VecND& corrZ) -{ - // propagate tracks to updated Z - for (int i = N; i--;) { - const auto& trDer = mTrDer[mCurHyp][i]; - auto dz2h = 0.5 * corrZ[i] * corrZ[i]; - mTrPos[mCurHyp][i][0] -= trDer.dxdz * corrZ[i] - dz2h * trDer.d2xdz2; - mTrPos[mCurHyp][i][1] -= trDer.dydz * corrZ[i] - dz2h * trDer.d2ydz2; - mTrPos[mCurHyp][i][2] -= corrZ[i]; - } - - return true; -} - -//___________________________________________________________________ -template -bool FwdDCAFitterN::FwdpropagateTracksToVertex(int icand) -{ - // propagate on z axis to vertex - int ord = mOrder[icand]; - if (mTrPropDone[ord]) { - return true; - } - const Vec3D& pca = mPCA[ord]; - for (int i = N; i--;) { - if (mUseAbsDCA) { - mCandTr[ord][i] = *mOrigTrPtr[i]; // fetch the track again, as mCandTr might have been propagated w/o errors - } - auto& trc = mCandTr[ord][i]; - auto z = pca[2]; - - trc.propagateToZquadratic(z, mBz); - } - - mTrPropDone[ord] = true; - return true; -} - -//___________________________________________________________________ -template -float FwdDCAFitterN::findZatXY(int mCurHyp) // Between 2 tracks -{ - - double step = 0.001; // initial step - double startPoint = 20.; // first MFT disk - - double z[2] = {startPoint, startPoint}; - double newX[2], newY[2]; - - double X = mPCA[mCurHyp][0]; // X seed - double Y = mPCA[mCurHyp][1]; // Y seed - - mCandTr[mCurHyp][0] = *mOrigTrPtr[0]; - mCandTr[mCurHyp][1] = *mOrigTrPtr[1]; - - double dstXY[2][3] = {{999., 999., 999.}, {999., 999., 999.}}; - - double Z[2]; - double finalZ[2]; - - double newDstXY; - - for (int i = 0; i < 2; i++) { - - while (z[i] > -10) { - - mCandTr[mCurHyp][i].propagateParamToZquadratic(z[i], mBz); - newX[i] = mCandTr[mCurHyp][i].getX(); - newY[i] = mCandTr[mCurHyp][i].getY(); - - newDstXY = std::sqrt((newX[i] - X) * (newX[i] - X) + - (newY[i] - Y) * (newY[i] - Y)); - - // Update points - dstXY[i][0] = dstXY[i][1]; - dstXY[i][1] = dstXY[i][2]; - dstXY[i][2] = newDstXY; - - if (dstXY[i][2] > dstXY[i][1] && dstXY[i][1] < dstXY[i][0]) { - finalZ[i] = z[i] + step; - break; - } - - z[i] -= step; - } - } - - float rez = 0.5 * (finalZ[0] + finalZ[1]); - return rez; -} - -//___________________________________________________________________ -template -void FwdDCAFitterN::findZatXY_mid(int mCurHyp) -{ - // look into dXY of T0 - T1 between 2 points(0,40cm); the one with the highest dXY is moved to mid - - double startPoint = -40.; - double endPoint = 50.; - double midPoint = 0.5 * (startPoint + endPoint); - - double z[2][2] = {{startPoint, endPoint}, {startPoint, endPoint}}; // z for tracks 0/1 on starting poing and endpoint - - double DeltaZ = std::abs(endPoint - startPoint); - - double newX[2][2]; - double newY[2][2]; - - double epsilon = 0.0001; - - double X = mPCA[mCurHyp][0]; // X seed - double Y = mPCA[mCurHyp][1]; // Y seed - - mCandTr[mCurHyp][0] = *mOrigTrPtr[0]; - mCandTr[mCurHyp][1] = *mOrigTrPtr[1]; - - double finalZ; - - double dstXY[2]; // 0 -> distance btwn both tracks at startPoint - - while (DeltaZ > epsilon) { - - midPoint = 0.5 * (startPoint + endPoint); - - for (int i = 0; i < 2; i++) { - mCandTr[mCurHyp][i].propagateParamToZquadratic(startPoint, mBz); - newX[i][0] = mCandTr[mCurHyp][i].getX(); - newY[i][0] = mCandTr[mCurHyp][i].getY(); - - mCandTr[mCurHyp][i].propagateParamToZquadratic(endPoint, mBz); - newX[i][1] = mCandTr[mCurHyp][i].getX(); - newY[i][1] = mCandTr[mCurHyp][i].getY(); - } - - dstXY[0] = (newX[0][0] - newX[1][0]) * (newX[0][0] - newX[1][0]) + - (newY[0][0] - newY[1][0]) * (newY[0][0] - newY[1][0]); - - dstXY[1] = (newX[0][1] - newX[1][1]) * (newX[0][1] - newX[1][1]) + - (newY[0][1] - newY[1][1]) * (newY[0][1] - newY[1][1]); - - DeltaZ = std::abs(endPoint - startPoint); - - if (DeltaZ < epsilon) { - finalZ = 0.5 * (startPoint + endPoint); - break; - } - - // chose new start and end Point according to the smallest D_XY - if (dstXY[1] > dstXY[0]) { - endPoint = midPoint; - } else { - startPoint = midPoint; - } - } - - mPCA[mCurHyp][2] = finalZ; -} - -//___________________________________________________________________ -template -void FwdDCAFitterN::findZatXY_lineApprox(int mCurHyp) -{ - // approx method: z=(b-b')/(a'-a) -> tracks to lines with y0,1=az0,1+b for each track (in YZ and XZ plane) - - double startPoint = 1.; - double endPoint = 50.; // first disk - - double X = mPCA[mCurHyp][0]; // X seed - double Y = mPCA[mCurHyp][1]; // Y seed - - mCandTr[mCurHyp][0] = *mOrigTrPtr[0]; - mCandTr[mCurHyp][1] = *mOrigTrPtr[1]; - - double y[2][2]; // Y00: y track 0 at point 0; Y01: y track 0 at point 1 - double z[2][2]; - double x[2][2]; - - double aYZ[2]; - double bYZ[2]; - - double aXZ[2]; - double bXZ[2]; - - double finalZ; - - // find points of the tracks = 2 straight lines - for (int i = 0; i < 2; i++) { - - mCandTr[mCurHyp][i].propagateToZquadratic(startPoint, mBz); - // mCandTr[mCurHyp][i].propagateToZlinear(startPoint); - z[i][0] = startPoint; - y[i][0] = mCandTr[mCurHyp][i].getY(); - x[i][0] = mCandTr[mCurHyp][i].getX(); - - mCandTr[mCurHyp][i].propagateToZquadratic(endPoint, mBz); - // mCandTr[mCurHyp][i].propagateToZlinear(endPoint); - z[i][1] = endPoint; - y[i][1] = mCandTr[mCurHyp][i].getY(); - x[i][1] = mCandTr[mCurHyp][i].getX(); - - bYZ[i] = (y[i][1] - y[i][0] * z[i][1] / z[i][0]) / (1 - z[i][1] / z[i][0]); - aYZ[i] = (y[i][0] - bYZ[i]) / z[i][0]; - - bXZ[i] = (x[i][1] - x[i][0] * z[i][1] / z[i][0]) / (1 - z[i][1] / z[i][0]); - aXZ[i] = (x[i][0] - bXZ[i]) / z[i][0]; - } - - // z seed: equ. for intersection of these lines - finalZ = 0.5 * ((bYZ[0] - bYZ[1]) / (aYZ[1] - aYZ[0]) + (bXZ[0] - bXZ[1]) / (aXZ[1] - aXZ[0])); - - mPCA[mCurHyp][2] = finalZ; -} - -//___________________________________________________________________ -template -void FwdDCAFitterN::findZatXY_quad(int mCurHyp) -{ - double startPoint = 0.; - double endPoint = 40.; // first disk - - double X = mPCA[mCurHyp][0]; // X seed - double Y = mPCA[mCurHyp][1]; // Y seed - - mCandTr[mCurHyp][0] = *mOrigTrPtr[0]; - mCandTr[mCurHyp][1] = *mOrigTrPtr[1]; - - double x[2]; - double y[2]; - double sinPhi0[2]; - double cosPhi0[2]; - double tanL0[2]; - double qpt0[2]; - - double k[2]; // B2C *abs(mBz) - double Hz[2]; // mBz/abs(mBz) - - double Ax[2], Bx[2], Cx[2]; - double Ay[2], By[2], Cy[2]; - - double deltaX[2], deltaY[2]; - - bool posX[2], nulX[2], negX[2]; - double z1X[2], z2X[2], z12X[2]; - - bool posY[2], nulY[2], negY[2]; - double z1Y[2], z2Y[2], z12Y[2]; - - double finalZ[2]; - - // find all variables for 2 tracks at z0 = startPoint - // set A, B, C variables for x/y equation for 2 tracks - // calculate Deltax/y for both and roots - - for (int i = 0; i < 2; i++) { - mCandTr[mCurHyp][i].propagateToZquadratic(startPoint, mBz); - x[i] = mCandTr[mCurHyp][i].getX(); - y[i] = mCandTr[mCurHyp][i].getY(); - sinPhi0[i] = mCandTr[mCurHyp][i].getSnp(); - cosPhi0[i] = std::sqrt((1. - sinPhi0[i]) * (1. + sinPhi0[i])); - tanL0[i] = mCandTr[mCurHyp][i].getTanl(); - qpt0[i] = mCandTr[mCurHyp][i].getInvQPt(); - k[i] = getK(mBz); - Hz[i] = getHz(mBz); - - Ax[i] = qpt0[i] * Hz[i] * k[i] * sinPhi0[i] / (2 * tanL0[i] * tanL0[i]); - Bx[i] = cosPhi0[i] / tanL0[i]; - Cx[i] = x[i] - X; - - Ay[i] = -qpt0[i] * Hz[i] * k[i] * cosPhi0[i] / (2 * tanL0[i] * tanL0[i]); - By[i] = sinPhi0[i] / tanL0[i]; - Cy[i] = y[i] - Y; // - - deltaX[i] = Bx[i] * Bx[i] - 4 * Ax[i] * Cx[i]; - deltaY[i] = By[i] * By[i] - 4 * Ay[i] * Cy[i]; - - if (deltaX[i] > 0) { - posX[i] = true; - z1X[i] = (-Bx[i] - std::sqrt(deltaX[i])) / (2 * Ax[i]); - z2X[i] = (-Bx[i] + std::sqrt(deltaX[i])) / (2 * Ax[i]); - } else if (deltaX[i] == 0) { - nulX[i] = true; - z12X[i] = -Bx[i] / (2 * Ax[i]); - } else { - negX[i] = true; - z12X[i] = 0; - } // discard - - if (deltaY[i] > 0) { - posY[i] = true; - z1Y[i] = (-By[i] - std::sqrt(deltaY[i])) / (2 * Ay[i]); - z2Y[i] = (-By[i] + std::sqrt(deltaY[i])) / (2 * Ay[i]); - } else if (deltaX[i] == 0) { - nulY[i] = true; - z12Y[i] = -By[i] / (2 * Ay[i]); - } else { - negY[i] = true; - z12Y[i] = 0; - } - - // find the z located in an acceptable interval - if (posX[i]) { - if (z1X[i] < endPoint && z1X[i] > startPoint) { - z12X[i] = z1X[i]; - } else { - z12X[i] = z2X[i]; - } - } - - if (posY[i]) { - if (z1Y[i] < endPoint && z1Y[i] > startPoint) { - z12Y[i] = z1Y[i]; - } else { - z12Y[i] = z2Y[i]; - } - } - - finalZ[i] = 0.5 * (z12X[i] + z12Y[i]); - } - - mPCA[mCurHyp][2] = 0.5 * (finalZ[0] + finalZ[1]); -} - -//___________________________________________________________________ -template -void FwdDCAFitterN::findZatXY_linear(int mCurHyp) -{ - - double startPoint = 0.; - - double X = mPCA[mCurHyp][0]; // X seed - double Y = mPCA[mCurHyp][1]; // Y seed - - mCandTr[mCurHyp][0] = *mOrigTrPtr[0]; - mCandTr[mCurHyp][1] = *mOrigTrPtr[1]; - - double x[2]; - double y[2]; - double sinPhi0[2]; - double cosPhi0[2]; - double tanL0[2]; - - double Ax[2], Bx[2]; - double Ay[2], By[2]; - - double z12X[2]; - double z12Y[2]; - - double finalZ[2]; - - // find all variables for 2 tracks at z0 = startPoint - // set A, B variables for x/y equation for 2 tracks - // calculate root - - for (int i = 0; i < 2; i++) { - mCandTr[mCurHyp][i].propagateToZlinear(startPoint); - x[i] = mCandTr[mCurHyp][i].getX(); - y[i] = mCandTr[mCurHyp][i].getY(); - sinPhi0[i] = mCandTr[mCurHyp][i].getSnp(); - cosPhi0[i] = std::sqrt((1. - sinPhi0[i]) * (1. + sinPhi0[i])); - tanL0[i] = mCandTr[mCurHyp][i].getTanl(); - - Ax[i] = cosPhi0[i] / tanL0[i]; - Bx[i] = x[i] - X; - - Ay[i] = sinPhi0[i] / tanL0[i]; - By[i] = y[i] - Y; - - z12X[i] = -Bx[i] / Ax[i]; - z12Y[i] = -By[i] / Ay[i]; - - finalZ[i] = 0.5 * (z12X[i] + z12Y[i]); - } - - mPCA[mCurHyp][2] = 0.5 * (finalZ[0] + finalZ[1]); -} - -//___________________________________________________________________ -template -inline o2::track::TrackParFwd FwdDCAFitterN::FwdgetTrackParamAtPCA(int i, int icand) const -{ - // propagate tracks param only to current vertex (if not already done) - int ord = mOrder[icand]; - o2::track::TrackParFwd trc(mCandTr[ord][i]); - if (!mTrPropDone[ord]) { - auto z = mPCA[ord][2]; - trc.propagateParamToZquadratic(z, mBz); - } - - return {trc}; -} - -//___________________________________________________________________ -template -inline double FwdDCAFitterN::getAbsMax(const VecND& v) -{ - double mx = -1; - for (int i = N; i--;) { - auto vai = std::abs(v[i]); - if (mx < vai) { - mx = vai; - } - } - return mx; -} - -//___________________________________________________________________ -template -bool FwdDCAFitterN::minimizeChi2() -{ - // find best chi2 (weighted DCA) of N tracks in the vicinity of the seed PCA - double x[2], y[2]; - double sumX = 0.; - double sumY = 0.; - - for (int i = N; i--;) { - mCandTr[mCurHyp][i] = *mOrigTrPtr[i]; - auto z = mPCA[mCurHyp][2]; - - mCandTr[mCurHyp][i].propagateToZquadratic(z, mBz); - - x[i] = mCandTr[mCurHyp][i].getX(); - y[i] = mCandTr[mCurHyp][i].getY(); - - setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions - mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], ZerrFactor); // prepare inverse cov.matrices at starting point - - sumX = sumX + x[i]; - sumY = sumY + y[i]; - } - - mPCA[mCurHyp][0] = sumX / N; - mPCA[mCurHyp][1] = sumY / N; - - if (mMaxDXIni > 0 && !roughDXCut()) { // apply rough cut on tracks X difference - return false; - } - - if (!FwdcalcPCACoefs()) { // prepare tracks contribution matrices to the global PCA - return false; - } - FwdcalcPCA(); // current PCA - FwdcalcTrackResiduals(); // current track residuals - float chi2Upd, chi2 = FwdcalcChi2(); - do { - calcTrackDerivatives(); // current track derivatives (1st and 2nd) - FwdcalcResidDerivatives(); // current residals derivatives (1st and 2nd) - FwdcalcChi2Derivatives(); // current chi2 derivatives (1st and 2nd) to proceed for dz calculation - - // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1 - if (!mD2Chi2Dz2.Invert()) { - return false; - } - - VecND dz = mD2Chi2Dz2 * mDChi2Dz; - - if (!FwdcorrectTracks(dz)) { // calculate new Pi (mTrPos) following Newton-Rapson iteration - return false; - } - - FwdcalcPCA(); // updated mPCA (new V coordinates with new mTrPos (Pi)) - if (mCrossIDAlt >= 0 && closerToAlternative()) { - mAllowAltPreference = false; - return false; - } - - FwdcalcTrackResiduals(); // updated residuals - chi2Upd = FwdcalcChi2(); // updated chi2 - - if (getAbsMax(dz) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) { - chi2 = chi2Upd; - break; // converged - } - - chi2 = chi2Upd; - } while (++mNIters[mCurHyp] < mMaxIter); - - mChi2[mCurHyp] = chi2 * NInv; - return mChi2[mCurHyp] < mMaxChi2; -} - -//___________________________________________________________________ -template -bool FwdDCAFitterN::minimizeChi2NoErr() -{ - // find best chi2 (absolute DCA) of N tracks in the vicinity of the PCA seed - double x[2], y[2]; - double sumX = 0.; - double sumY = 0.; - - for (int i = N; i--;) { - - mCandTr[mCurHyp][i] = *mOrigTrPtr[i]; - - auto z = mPCA[mCurHyp][2]; - mCandTr[mCurHyp][i].propagateParamToZquadratic(z, mBz); - - x[i] = mCandTr[mCurHyp][i].getX(); - y[i] = mCandTr[mCurHyp][i].getY(); - - mPCA[mCurHyp][2] = z; - - setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions - - sumX = sumX + x[i]; - sumY = sumY + y[i]; - } - - mPCA[mCurHyp][0] = sumX / N; - mPCA[mCurHyp][1] = sumY / N; - - if (mMaxDXIni > 0 && !roughDXCut()) { // apply rough cut on tracks Z difference - return false; - } - - FwdcalcPCANoErr(); // current PCA - FwdcalcTrackResiduals(); // current track residuals - float chi2Upd, chi2 = FwdcalcChi2NoErr(); - do { - calcTrackDerivatives(); // current track derivatives (1st and 2nd) - FwdcalcResidDerivativesNoErr(); // current residals derivatives (1st and 2nd) - FwdcalcChi2DerivativesNoErr(); // current chi2 derivatives (1st and 2nd) - - // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1 - if (!mD2Chi2Dz2.Invert()) { - return false; - } - VecND dz = mD2Chi2Dz2 * mDChi2Dz; - - if (!FwdcorrectTracks(dz)) { - return false; - } - FwdcalcPCANoErr(); // updated PCA - if (mCrossIDAlt >= 0 && closerToAlternative()) { - mAllowAltPreference = false; - return false; - } - FwdcalcTrackResiduals(); // updated residuals - chi2Upd = FwdcalcChi2NoErr(); // updated chi2 - if (getAbsMax(dz) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) { - chi2 = chi2Upd; - break; // converged - } - chi2 = chi2Upd; - } while (++mNIters[mCurHyp] < mMaxIter); - // - mChi2[mCurHyp] = chi2 * NInv; - return mChi2[mCurHyp] < mMaxChi2; -} - -//___________________________________________________________________ -template -bool FwdDCAFitterN::roughDXCut() const -{ - // apply rough cut on DX between the tracks in the seed point - - bool accept = true; - for (int i = N; accept && i--;) { - for (int j = i; j--;) { - if (std::abs(mCandTr[mCurHyp][i].getX() - mCandTr[mCurHyp][j].getX()) > mMaxDXIni) { - accept = false; - break; - } - } - } - return accept; -} - -//___________________________________________________________________ -template -bool FwdDCAFitterN::closerToAlternative() const -{ - // check if the point current PCA point is closer to the seeding XY point being tested or to alternative see (if any) - auto dxCur = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDCur], dyCur = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDCur]; - auto dxAlt = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDAlt], dyAlt = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDAlt]; - return dxCur * dxCur + dyCur * dyCur > dxAlt * dxAlt + dyAlt * dyAlt; -} - -//___________________________________________________________________ -template -void FwdDCAFitterN::print() const -{ - LOG(info) << N << "-prong vertex fitter in " << (mUseAbsDCA ? "abs." : "weighted") << " distance minimization mode"; - LOG(info) << "Bz: " << mBz << " MaxIter: " << mMaxIter << " MaxChi2: " << mMaxChi2; - LOG(info) << "Stopping condition: Max.param change < " << mMinParamChange << " Rel.Chi2 change > " << mMinRelChi2Change; - LOG(info) << "Discard candidates for : Rvtx > " << getMaxR() << " DZ between tracks > " << mMaxDXIni; -} - -using FwdDCAFitter2 = FwdDCAFitterN<2, o2::track::TrackParCovFwd>; -using FwdDCAFitter3 = FwdDCAFitterN<3, o2::track::TrackParCovFwd>; - -} // namespace vertexing -} // namespace o2 -#endif // _ALICEO2_DCA_FWDFITTERN_ diff --git a/Common/DCAFitter/include/DCAFitter/HelixHelper.h b/Common/DCAFitter/include/DCAFitter/HelixHelper.h deleted file mode 100644 index 87575367e1af5..0000000000000 --- a/Common/DCAFitter/include/DCAFitter/HelixHelper.h +++ /dev/null @@ -1,282 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file HelixHelper.h -/// \brief Helper classes for helical tracks manipulations -/// \author ruben.shahoyan@cern.ch - -#ifndef _ALICEO2_HELIX_HELPER_ -#define _ALICEO2_HELIX_HELPER_ - -#include "CommonConstants/MathConstants.h" -#include "MathUtils/Utils.h" -#include "MathUtils/Primitive2D.h" - -namespace o2 -{ -namespace track -{ - -///__________________________________________________________________________ -//< precalculated track radius, center, alpha sin,cos and their combinations -struct TrackAuxPar : public o2::math_utils::CircleXYf_t { - float c, s, cc, ss, cs; // cos ans sin of track alpha and their products - - TrackAuxPar() = default; - - template - TrackAuxPar(const T& trc, float bz) - { - set(trc, bz); - } - float cosDif(const TrackAuxPar& t) const { return c * t.c + s * t.s; } // cos(alpha_this - alha_t) - float sinDif(const TrackAuxPar& t) const { return s * t.c - c * t.s; } // sin(alpha_this - alha_t) - - template - void set(const T& trc, float bz) - { - trc.getCircleParams(bz, *this, s, c); - cc = c * c; - ss = s * s; - cs = c * s; - } - ClassDefNV(TrackAuxPar, 1); -}; - -//__________________________________________________________ -//< crossing coordinates of 2 circles -struct CrossInfo { - static constexpr float MaxDistXYDef = 10.; - float xDCA[2] = {}; - float yDCA[2] = {}; - int nDCA = 0; - - int circlesCrossInfo(const TrackAuxPar& trax0, const TrackAuxPar& trax1, float maxDistXY = MaxDistXYDef) - { - const auto& trcA = trax0.rC > trax1.rC ? trax0 : trax1; // designate the largest circle as A - const auto& trcB = trax0.rC > trax1.rC ? trax1 : trax0; - float xDist = trcB.xC - trcA.xC, yDist = trcB.yC - trcA.yC; - float dist2 = xDist * xDist + yDist * yDist, dist = std::sqrt(dist2), rsum = trcA.rC + trcB.rC; - if (std::abs(dist) < 1e-12) { - return nDCA; // circles are concentric? - } - if (dist > rsum) { // circles don't touch, chose a point in between - // the parametric equation of lines connecting the centers is - // x = x0 + t/dist * (x1-x0), y = y0 + t/dist * (y1-y0) - if (dist - rsum > maxDistXY) { // too large distance - return nDCA; - } - notTouchingXY(dist, xDist, yDist, trcA, trcB.rC); - } else if (dist + trcB.rC < trcA.rC) { // the small circle is nestled into large one w/o touching - // select the point of closest approach of 2 circles - notTouchingXY(dist, xDist, yDist, trcA, -trcB.rC); - } else { // 2 intersection points - // to simplify calculations, we move to new frame x->x+Xc0, y->y+Yc0, so that - // the 1st one is centered in origin - if (std::abs(xDist) < std::abs(yDist)) { - float a = (trcA.rC * trcA.rC - trcB.rC * trcB.rC + dist2) / (2. * yDist), b = -xDist / yDist, ab = a * b, bb = b * b; - float det = ab * ab - (1. + bb) * (a * a - trcA.rC * trcA.rC); - if (det > 0.) { - det = std::sqrt(det); - xDCA[0] = (-ab + det) / (1. + b * b); - yDCA[0] = a + b * xDCA[0] + trcA.yC; - xDCA[0] += trcA.xC; - xDCA[1] = (-ab - det) / (1. + b * b); - yDCA[1] = a + b * xDCA[1] + trcA.yC; - xDCA[1] += trcA.xC; - nDCA = 2; - } else { // due to the finite precision the det<=0, i.e. the circles are barely touching, fall back to this special case - notTouchingXY(dist, xDist, yDist, trcA, trcB.rC); - } - } else { - float a = (trcA.rC * trcA.rC - trcB.rC * trcB.rC + dist2) / (2. * xDist), b = -yDist / xDist, ab = a * b, bb = b * b; - float det = ab * ab - (1. + bb) * (a * a - trcA.rC * trcA.rC); - if (det > 0.) { - det = std::sqrt(det); - yDCA[0] = (-ab + det) / (1. + bb); - xDCA[0] = a + b * yDCA[0] + trcA.xC; - yDCA[0] += trcA.yC; - yDCA[1] = (-ab - det) / (1. + bb); - xDCA[1] = a + b * yDCA[1] + trcA.xC; - yDCA[1] += trcA.yC; - nDCA = 2; - } else { // due to the finite precision the det<=0, i.e. the circles are barely touching, fall back to this special case - notTouchingXY(dist, xDist, yDist, trcA, trcB.rC); - } - } - } - return nDCA; - } - - void notTouchingXY(float dist, float xDist, float yDist, const TrackAuxPar& trcA, float rBSign) - { - // fast method to calculate DCA between 2 circles, assuming that they don't touch each outer: - // the parametric equation of lines connecting the centers is x = xA + t/dist * xDist, y = yA + t/dist * yDist - // with xA,yY being the center of the circle A ( = trcA.xC, trcA.yC ), xDist = trcB.xC = trcA.xC ... - // There are 2 special cases: - // (a) small circle is inside the large one: provide rBSign as -trcB.rC - // (b) circle are side by side: provide rBSign as trcB.rC - nDCA = 1; - auto t2d = (dist + trcA.rC - rBSign) / dist; - xDCA[0] = trcA.xC + 0.5 * (xDist * t2d); - yDCA[0] = trcA.yC + 0.5 * (yDist * t2d); - } - - template - int linesCrossInfo(const TrackAuxPar& trax0, const T& tr0, - const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef) - { - /// closest approach of 2 straight lines - /// TrackParam propagation can be parameterized in lab in a form - /// xLab(t) = (x*cosAlp - y*sinAlp) + t*(cosAlp - sinAlp* snp/csp) = xLab0 + t*(cosAlp - sinAlp* snp/csp) - /// yLab(t) = (x*sinAlp + y*cosAlp) + t*(sinAlp + cosAlp* snp/csp) = yLab0 + t*(sinAlp + cosAlp* snp/csp) - /// zLab(t) = z + t * tgl / csp = zLab0 + t * tgl / csp - /// where t is the x-step in the track alpha-frame, xLab,yLab,zLab are reference track coordinates in lab - /// frame (filled by TrackAuxPar for straight line tracks). - /// - /// Therefore, for the parametric track equation in lab 3D we have (wrt tracking-X increment t) - /// xL(t) = xL + t Kx; Kx = (cosAlp - sinAlp* snp/csp) - /// yL(t) = yL + t Ky; Ky = (sinAlp + cosAlp* snp/csp) - /// zL(t) = zL + t Kz; Kz = tgl / csp - /// Note that Kx^2 + Ky^2 + Kz^2 = (1+tgl^2) / csp^2 - - float dx = trax1.xC - trax0.xC; // for straight line TrackAuxPar stores lab coordinates at referene point!!! - float dy = trax1.yC - trax0.yC; // - float dz = tr1.getZ() - tr0.getZ(); - auto csp0i2 = 1. / tr0.getCsp2(); // 1 / csp^2 - auto csp0i = std::sqrt(csp0i2); - auto tgp0 = tr0.getSnp() * csp0i; - float kx0 = trax0.c - trax0.s * tgp0; - float ky0 = trax0.s + trax0.c * tgp0; - float kz0 = tr0.getTgl() * csp0i; - auto csp1i2 = 1. / tr1.getCsp2(); // 1 / csp^2 - auto csp1i = std::sqrt(csp1i2); - auto tgp1 = tr1.getSnp() * std::sqrt(csp1i2); - float kx1 = trax1.c - trax1.s * tgp1; - float ky1 = trax1.s + trax1.c * tgp1; - float kz1 = tr1.getTgl() * csp1i; - /// Minimize |vecL1 - vecL0|^2 wrt t0 and t1: point of closest approach - /// Leads to system - /// A Dx = B with Dx = {dx0, dx1} - /// with A = - /// | kx0^2+ky0^2+kz0^2 -(kx0*kx1+ky0*ky1+kz0*kz1) | = (1+tgl0^2) / csp0^2 .... - /// | -(kx0*kx1+ky0*ky1+kz0*kz1) kx0^2+ky0^2+kz0^2 | ..... (1+tgl1^2) / csp1^2 - /// and B = {(dx Kx0 + dy Ky0 + dz Kz0), -(dx Kx1 + dy Ky1 + dz Kz1) } - /// - float a00 = (1.f + tr0.getTgl() * tr0.getTgl()) * csp0i2, a11 = (1.f + tr1.getTgl() * tr1.getTgl()) * csp1i2, a01 = -(kx0 * kx1 + ky0 * ky1 + kz0 * kz1); - float b0 = dx * kx0 + dy * ky0 + dz * kz0, b1 = -(dx * kx1 + dy * ky1 + dz * kz1); - float det = a00 * a11 - a01 * a01, det0 = b0 * a11 - b1 * a01, det1 = a00 * b1 - a01 * b0; - if (std::abs(det) > o2::constants::math::Almost0) { - auto detI = 1. / det; - auto t0 = det0 * detI; - auto t1 = det1 * detI; - float addx0 = kx0 * t0, addy0 = ky0 * t0, addx1 = kx1 * t1, addy1 = ky1 * t1; - dx += addx1 - addx0; // recalculate XY distance at DCA - dy += addy1 - addy0; - if (dx * dx + dy * dy > maxDistXY * maxDistXY) { - return nDCA; - } - xDCA[0] = (trax0.xC + addx0 + trax1.xC + addx1) * 0.5; - yDCA[0] = (trax0.yC + addy0 + trax1.yC + addy1) * 0.5; - nDCA = 1; - } - return nDCA; - } - - template - int circleLineCrossInfo(const TrackAuxPar& trax0, const T& tr0, - const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef) - { - /// closest approach of line and circle - /// TrackParam propagation can be parameterized in lab in a form - /// xLab(t) = (x*cosAlp - y*sinAlp) + t*(cosAlp - sinAlp* snp/csp) = xLab0 + t*(cosAlp - sinAlp* snp/csp) - /// yLab(t) = (x*sinAlp + y*cosAlp) + t*(sinAlp + cosAlp* snp/csp) = yLab0 + t*(sinAlp + cosAlp* snp/csp) - /// zLab(t) = z + t * tgl / csp = zLab0 + t * tgl / csp - /// where t is the x-step in the track alpha-frame, xLab,yLab,zLab are reference track coordinates in lab - /// frame (filled by TrackAuxPar for straight line tracks). - /// - /// Therefore, for the parametric track equation in lab 3D we have (wrt tracking-X increment t) - /// xL(t) = xL + t Kx; Kx = (cosAlp - sinAlp* snp/csp) - /// yL(t) = yL + t Ky; Ky = (sinAlp + cosAlp* snp/csp) - /// zL(t) = zL + t Kz; Kz = tgl / csp - /// Note that Kx^2 + Ky^2 = 1 / csp^2 - - const auto& traxH = trax0.rC > trax1.rC ? trax0 : trax1; // circle (for the line rC is set to 0) - const auto& traxL = trax0.rC > trax1.rC ? trax1 : trax0; // line - const auto& trcL = trax0.rC > trax1.rC ? tr1 : tr0; // track of the line - - // solve quadratic equation of line crossing the circle - float dx = traxL.xC - traxH.xC; // X distance between the line lab reference and circle center - float dy = traxL.yC - traxH.yC; // Y... - // t^2(kx^2+ky^2) + 2t(dx*kx+dy*ky) + dx^2 + dy^2 - r^2 = 0 - auto cspi2 = 1. / trcL.getCsp2(); // 1 / csp^2 == kx^2 + ky^2 - auto cspi = std::sqrt(cspi2); - auto tgp = trcL.getSnp() * cspi; - float kx = traxL.c - traxL.s * tgp; - float ky = traxL.s + traxL.c * tgp; - double dk = dx * kx + dy * ky; - double det = dk * dk - cspi2 * (dx * dx + dy * dy - traxH.rC * traxH.rC); - if (det > 0) { // 2 crossings - det = std::sqrt(det); - float t0 = (-dk + det) * cspi2; - float t1 = (-dk - det) * cspi2; - xDCA[0] = traxL.xC + kx * t0; - yDCA[0] = traxL.yC + ky * t0; - xDCA[1] = traxL.xC + kx * t1; - yDCA[1] = traxL.yC + ky * t1; - nDCA = 2; - } else { - // there is no crossing, find the point of the closest approach on the line which is closest to the circle center - float t = -dk * cspi2; - float xL = traxL.xC + kx * t, yL = traxL.yC + ky * t; // point on the line, need to average with point on the circle - float dxc = xL - traxH.xC, dyc = yL - traxH.yC, dist = std::sqrt(dxc * dxc + dyc * dyc); - if (dist - traxH.rC > maxDistXY) { // too large distance - return nDCA; - } - float drcf = traxH.rC / dist; // radius / distance to circle center - float xH = traxH.xC + dxc * drcf, yH = traxH.yC + dyc * drcf; - xDCA[0] = (xL + xH) * 0.5; - yDCA[0] = (yL + yH) * 0.5; - nDCA = 1; - } - return nDCA; - } - - template - int set(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef) - { - // calculate up to 2 crossings between 2 circles - nDCA = 0; - if (trax0.rC > o2::constants::math::Almost0 && trax1.rC > o2::constants::math::Almost0) { // both are not straight lines - nDCA = circlesCrossInfo(trax0, trax1, maxDistXY); - } else if (trax0.rC < o2::constants::math::Almost0 && trax1.rC < o2::constants::math::Almost0) { // both are straigt lines - nDCA = linesCrossInfo(trax0, tr0, trax1, tr1, maxDistXY); - } else { - nDCA = circleLineCrossInfo(trax0, tr0, trax1, tr1, maxDistXY); - } - // - return nDCA; - } - - CrossInfo() = default; - - template - CrossInfo(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef) - { - set(trax0, tr0, trax1, tr1, maxDistXY); - } - ClassDefNV(CrossInfo, 1); -}; - -} // namespace track -} // namespace o2 - -#endif diff --git a/Common/DCAFitter/src/DCAFitterLinkDef.h b/Common/DCAFitter/src/DCAFitterLinkDef.h deleted file mode 100644 index 3589ffe559e96..0000000000000 --- a/Common/DCAFitter/src/DCAFitterLinkDef.h +++ /dev/null @@ -1,27 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -#ifdef __CLING__ - -#pragma link off all globals; -#pragma link off all classes; -#pragma link off all functions; - -#pragma link C++ class o2::vertexing::DCAFitterN < 2, o2::track::TrackParCov> + ; -#pragma link C++ class o2::vertexing::DCAFitterN < 3, o2::track::TrackParCov> + ; - -#pragma link C++ class o2::track::TrackAuxPar + ; -#pragma link C++ class o2::track::CrossInfo + ; - -#pragma link C++ function o2::vertexing::DCAFitter2::process(const o2::track::TrackParCov&, const o2::track::TrackParCov&); -#pragma link C++ function o2::vertexing::DCAFitter3::process(const o2::track::TrackParCov&, const o2::track::TrackParCov&, const o2::track::TrackParCov&); - -#endif diff --git a/Common/DCAFitter/src/DCAFitterN.cxx b/Common/DCAFitter/src/DCAFitterN.cxx deleted file mode 100644 index dbc992bcf687d..0000000000000 --- a/Common/DCAFitter/src/DCAFitterN.cxx +++ /dev/null @@ -1,33 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file DCAFitterN.cxx -/// \brief Defintions for N-prongs secondary vertex fit -/// \author ruben.shahoyan@cern.ch - -#include "DCAFitter/DCAFitterN.h" - -namespace o2 -{ -namespace vertexing -{ - -void __dummy_instance__() -{ - DCAFitter2 ft2; - DCAFitter3 ft3; - o2::track::TrackParCov tr; - ft2.process(tr, tr); - ft3.process(tr, tr, tr); -} - -} // namespace vertexing -} // namespace o2 diff --git a/Common/DCAFitter/src/FwdDCAFitterN.cxx b/Common/DCAFitter/src/FwdDCAFitterN.cxx deleted file mode 100644 index 3198c6dfbfd52..0000000000000 --- a/Common/DCAFitter/src/FwdDCAFitterN.cxx +++ /dev/null @@ -1,33 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file DCAFitterN.cxx -/// \brief Defintions for N-prongs secondary vertex fit -/// \author ruben.shahoyan@cern.ch, adapted from central barrel to fwd rapidities by Rita Sadek, rita.sadek@cern.ch - -#include "DCAFitter/FwdDCAFitterN.h" - -namespace o2 -{ -namespace vertexing -{ - -void __test_instance__() -{ - FwdDCAFitter2 ft2; - FwdDCAFitter3 ft3; - o2::track::TrackParCovFwd tr; - ft2.process(tr, tr); - ft3.process(tr, tr, tr); -} - -} // namespace vertexing -} // namespace o2 diff --git a/Common/DCAFitter/test/testDCAFitterN.cxx b/Common/DCAFitter/test/testDCAFitterN.cxx deleted file mode 100644 index 7fef1f6889284..0000000000000 --- a/Common/DCAFitter/test/testDCAFitterN.cxx +++ /dev/null @@ -1,471 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -#define BOOST_TEST_MODULE Test DCAFitterN class -#define BOOST_TEST_MAIN -#define BOOST_TEST_DYN_LINK -#include - -#include "DCAFitter/DCAFitterN.h" -#include "CommonUtils/TreeStreamRedirector.h" -#include -#include -#include -#include -#include -#include - -namespace o2 -{ -namespace vertexing -{ - -using Vec3D = ROOT::Math::SVector; - -template -float checkResults(o2::utils::TreeStreamRedirector& outs, std::string& treeName, FITTER& fitter, - Vec3D& vgen, TLorentzVector& genPar, const std::vector& dtMass) -{ - int nCand = fitter.getNCandidates(); - std::array p; - float distMin = 1e9; - bool absDCA = fitter.getUseAbsDCA(); - bool useWghDCA = fitter.getWeightedFinalPCA(); - for (int ic = 0; ic < nCand; ic++) { - const auto& vtx = fitter.getPCACandidate(ic); - auto df = vgen; - df -= vtx; - - TLorentzVector moth, prong; - for (int i = 0; i < fitter.getNProngs(); i++) { - const auto& trc = fitter.getTrack(i, ic); - trc.getPxPyPzGlo(p); - prong.SetVectM({p[0], p[1], p[2]}, dtMass[i]); - moth += prong; - } - auto nIter = fitter.getNIterations(ic); - auto chi2 = fitter.getChi2AtPCACandidate(ic); - double dst = TMath::Sqrt(df[0] * df[0] + df[1] * df[1] + df[2] * df[2]); - distMin = dst < distMin ? dst : distMin; - auto parentTrack = fitter.createParentTrackParCov(ic); - // float genX - outs << treeName.c_str() << "cand=" << ic << "ncand=" << nCand << "nIter=" << nIter << "chi2=" << chi2 - << "genPart=" << genPar << "recPart=" << moth - << "genX=" << vgen[0] << "genY=" << vgen[1] << "genZ=" << vgen[2] - << "dx=" << df[0] << "dy=" << df[1] << "dz=" << df[2] << "dst=" << dst - << "useAbsDCA=" << absDCA << "useWghDCA=" << useWghDCA << "parent=" << parentTrack << "\n"; - } - return distMin; -} - -TLorentzVector generate(Vec3D& vtx, std::vector& vctr, float bz, - TGenPhaseSpace& genPHS, double parMass, const std::vector& dtMass, std::vector forceQ) -{ - const float errYZ = 1e-2, errSlp = 1e-3, errQPT = 2e-2; - std::array covm = { - errYZ * errYZ, - 0., errYZ * errYZ, - 0, 0., errSlp * errSlp, - 0., 0., 0., errSlp * errSlp, - 0., 0., 0., 0., errQPT * errQPT}; - bool accept = true; - TLorentzVector parent, d0, d1, d2; - do { - accept = true; - double y = gRandom->Rndm() - 0.5; - double pt = 0.1 + gRandom->Rndm() * 3; - double mt = TMath::Sqrt(parMass * parMass + pt * pt); - double pz = mt * TMath::SinH(y); - double phi = gRandom->Rndm() * TMath::Pi() * 2; - double en = mt * TMath::CosH(y); - double rdec = 10.; // radius of the decay - vtx[0] = rdec * TMath::Cos(phi); - vtx[1] = rdec * TMath::Sin(phi); - vtx[2] = rdec * pz / pt; - parent.SetPxPyPzE(pt * TMath::Cos(phi), pt * TMath::Sin(phi), pz, en); - int nd = dtMass.size(); - genPHS.SetDecay(parent, nd, dtMass.data()); - genPHS.Generate(); - vctr.clear(); - float p[4]; - for (int i = 0; i < nd; i++) { - auto* dt = genPHS.GetDecay(i); - if (dt->Pt() < 0.05) { - accept = false; - break; - } - dt->GetXYZT(p); - float s, c, x; - std::array params; - o2::math_utils::sincos(dt->Phi(), s, c); - o2::math_utils::rotateZInv(vtx[0], vtx[1], x, params[0], s, c); - - params[1] = vtx[2]; - params[2] = 0.; // since alpha = phi - params[3] = 1. / TMath::Tan(dt->Theta()); - params[4] = (i % 2 ? -1. : 1.) / dt->Pt(); - covm[14] = errQPT * errQPT * params[4] * params[4]; - // - // randomize - float r1, r2; - gRandom->Rannor(r1, r2); - params[0] += r1 * errYZ; - params[1] += r2 * errYZ; - gRandom->Rannor(r1, r2); - params[2] += r1 * errSlp; - params[3] += r2 * errSlp; - params[4] *= gRandom->Gaus(1., errQPT); - if (forceQ[i] == 0) { - params[4] = 0.; // impose straight track - } - auto& trc = vctr.emplace_back(x, dt->Phi(), params, covm); - float rad = forceQ[i] == 0 ? 600. : TMath::Abs(1. / trc.getCurvature(bz)); - if (!trc.propagateTo(trc.getX() + (gRandom->Rndm() - 0.5) * rad * 0.05, bz) || - !trc.rotate(trc.getAlpha() + (gRandom->Rndm() - 0.5) * 0.2)) { - printf("Failed to randomize "); - trc.print(); - } - } - } while (!accept); - - return parent; -} - -BOOST_AUTO_TEST_CASE(DCAFitterNProngs) -{ - constexpr int NTest = 10000; - o2::utils::TreeStreamRedirector outStream("dcafitterNTest.root"); - - TGenPhaseSpace genPHS; - constexpr double pion = 0.13957; - constexpr double k0 = 0.49761; - constexpr double kch = 0.49368; - constexpr double dch = 1.86965; - std::vector k0dec = {pion, pion}; - std::vector dchdec = {pion, kch, pion}; - std::vector vctracks; - Vec3D vtxGen; - - double bz = 5.0; - // 2 prongs vertices - { - LOG(info) << "Processing 2-prong Helix - Helix case"; - std::vector forceQ{1, 1}; - - o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter - ft.setBz(bz); - ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - ft.setMaxDXYIni(4); // do not consider V0 seeds with tracks XY-distance exceeding this. This is default anyway - ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - - std::string treeName2A = "pr2a", treeName2AW = "pr2aw", treeName2W = "pr2w"; - TStopwatch swA, swAW, swW; - int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - double meanDA = 0, meanDAW = 0, meanDW = 0; - swA.Stop(); - swAW.Stop(); - swW.Stop(); - for (int iev = 0; iev < NTest; iev++) { - auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); - - ft.setUseAbsDCA(true); - swA.Start(false); - int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - swA.Stop(); - LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - if (ncA) { - auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); - meanDA += minD; - nfoundA++; - } - - ft.setUseAbsDCA(true); - ft.setWeightedFinalPCA(true); - swAW.Start(false); - int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - swAW.Stop(); - LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - if (ncAW) { - auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); - meanDAW += minD; - nfoundAW++; - } - - ft.setUseAbsDCA(false); - ft.setWeightedFinalPCA(false); - swW.Start(false); - int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - swW.Stop(); - LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - if (ncW) { - auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); - meanDW += minD; - nfoundW++; - } - } - ft.print(); - meanDA /= nfoundA ? nfoundA : 1; - meanDAW /= nfoundA ? nfoundA : 1; - meanDW /= nfoundW ? nfoundW : 1; - LOG(info) << "Processed " << NTest << " 2-prong vertices Helix : Helix"; - LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - BOOST_CHECK(nfoundA > 0.99 * NTest); - BOOST_CHECK(nfoundAW > 0.99 * NTest); - BOOST_CHECK(nfoundW > 0.99 * NTest); - BOOST_CHECK(meanDA < 0.1); - BOOST_CHECK(meanDAW < 0.1); - BOOST_CHECK(meanDW < 0.1); - } - - // 2 prongs vertices with one of charges set to 0: Helix : Line - { - std::vector forceQ{1, 1}; - LOG(info) << "Processing 2-prong Helix - Line case"; - o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter - ft.setBz(bz); - ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - - std::string treeName2A = "pr2aHL", treeName2AW = "pr2awHL", treeName2W = "pr2wHL"; - TStopwatch swA, swAW, swW; - int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - double meanDA = 0, meanDAW = 0, meanDW = 0; - swA.Stop(); - swAW.Stop(); - swW.Stop(); - for (int iev = 0; iev < NTest; iev++) { - forceQ[iev % 2] = 1; - forceQ[1 - iev % 2] = 0; - auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); - - ft.setUseAbsDCA(true); - swA.Start(false); - int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - swA.Stop(); - LOG(debug) << "fit abs.dist with final weighted DCA " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - if (ncA) { - auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); - meanDA += minD; - nfoundA++; - } - - ft.setUseAbsDCA(true); - ft.setWeightedFinalPCA(true); - swAW.Start(false); - int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - swAW.Stop(); - LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - if (ncAW) { - auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); - meanDAW += minD; - nfoundAW++; - } - - ft.setUseAbsDCA(false); - ft.setWeightedFinalPCA(false); - swW.Start(false); - int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - swW.Stop(); - LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - if (ncW) { - auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); - meanDW += minD; - nfoundW++; - } - } - ft.print(); - meanDA /= nfoundA ? nfoundA : 1; - meanDAW /= nfoundAW ? nfoundAW : 1; - meanDW /= nfoundW ? nfoundW : 1; - LOG(info) << "Processed " << NTest << " 2-prong vertices: Helix : Line"; - LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - BOOST_CHECK(nfoundA > 0.99 * NTest); - BOOST_CHECK(nfoundAW > 0.99 * NTest); - BOOST_CHECK(nfoundW > 0.99 * NTest); - BOOST_CHECK(meanDA < 0.1); - BOOST_CHECK(meanDAW < 0.1); - BOOST_CHECK(meanDW < 0.1); - } - - // 2 prongs vertices with both of charges set to 0: Line : Line - { - std::vector forceQ{0, 0}; - LOG(info) << "Processing 2-prong Line - Line case"; - o2::vertexing::DCAFitterN<2> ft; // 2 prong fitter - ft.setBz(bz); - ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - - std::string treeName2A = "pr2aLL", treeName2AW = "pr2awLL", treeName2W = "pr2wLL"; - TStopwatch swA, swAW, swW; - int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - double meanDA = 0, meanDAW = 0, meanDW = 0; - swA.Stop(); - swAW.Stop(); - swW.Stop(); - for (int iev = 0; iev < NTest; iev++) { - forceQ[0] = forceQ[1] = 0; - auto genParent = generate(vtxGen, vctracks, bz, genPHS, k0, k0dec, forceQ); - - ft.setUseAbsDCA(true); - swA.Start(false); - int ncA = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - swA.Stop(); - LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - if (ncA) { - auto minD = checkResults(outStream, treeName2A, ft, vtxGen, genParent, k0dec); - meanDA += minD; - nfoundA++; - } - - ft.setUseAbsDCA(true); - ft.setWeightedFinalPCA(true); - swAW.Start(false); - int ncAW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - swAW.Stop(); - LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - if (ncAW) { - auto minD = checkResults(outStream, treeName2AW, ft, vtxGen, genParent, k0dec); - meanDAW += minD; - nfoundAW++; - } - - ft.setUseAbsDCA(false); - ft.setWeightedFinalPCA(false); - swW.Start(false); - int ncW = ft.process(vctracks[0], vctracks[1]); // HERE WE FIT THE VERTICES - swW.Stop(); - LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - if (ncW) { - auto minD = checkResults(outStream, treeName2W, ft, vtxGen, genParent, k0dec); - meanDW += minD; - nfoundW++; - } - } - ft.print(); - meanDA /= nfoundA ? nfoundA : 1; - meanDAW /= nfoundAW ? nfoundAW : 1; - meanDW /= nfoundW ? nfoundW : 1; - LOG(info) << "Processed " << NTest << " 2-prong vertices: Line : Line"; - LOG(info) << "2-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - LOG(info) << "2-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - LOG(info) << "2-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - BOOST_CHECK(nfoundA > 0.99 * NTest); - BOOST_CHECK(nfoundAW > 0.99 * NTest); - BOOST_CHECK(nfoundW > 0.99 * NTest); - BOOST_CHECK(meanDA < 0.1); - BOOST_CHECK(meanDAW < 0.1); - BOOST_CHECK(meanDW < 0.1); - } - - // 3 prongs vertices - { - std::vector forceQ{1, 1, 1}; - - o2::vertexing::DCAFitterN<3> ft; // 3 prong fitter - ft.setBz(bz); - ft.setPropagateToPCA(true); // After finding the vertex, propagate tracks to the DCA. This is default anyway - ft.setMaxR(200); // do not consider V0 seeds with 2D circles crossing above this R. This is default anyway - ft.setMaxDZIni(4); // do not consider V0 seeds with tracks Z-distance exceeding this. This is default anyway - ft.setMinParamChange(1e-3); // stop iterations if max correction is below this value. This is default anyway - ft.setMinRelChi2Change(0.9); // stop iterations if chi2 improves by less that this factor - - std::string treeName3A = "pr3a", treeName3AW = "pr3aw", treeName3W = "pr3w"; - TStopwatch swA, swAW, swW; - int nfoundA = 0, nfoundAW = 0, nfoundW = 0; - double meanDA = 0, meanDAW = 0, meanDW = 0; - swA.Stop(); - swAW.Stop(); - swW.Stop(); - for (int iev = 0; iev < NTest; iev++) { - auto genParent = generate(vtxGen, vctracks, bz, genPHS, dch, dchdec, forceQ); - - ft.setUseAbsDCA(true); - swA.Start(false); - int ncA = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES - swA.Stop(); - LOG(debug) << "fit abs.dist " << iev << " NC: " << ncA << " Chi2: " << (ncA ? ft.getChi2AtPCACandidate(0) : -1); - if (ncA) { - auto minD = checkResults(outStream, treeName3A, ft, vtxGen, genParent, dchdec); - meanDA += minD; - nfoundA++; - } - - ft.setUseAbsDCA(true); - ft.setWeightedFinalPCA(true); - swAW.Start(false); - int ncAW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES - swAW.Stop(); - LOG(debug) << "fit abs.dist " << iev << " NC: " << ncAW << " Chi2: " << (ncAW ? ft.getChi2AtPCACandidate(0) : -1); - if (ncAW) { - auto minD = checkResults(outStream, treeName3AW, ft, vtxGen, genParent, dchdec); - meanDAW += minD; - nfoundAW++; - } - - ft.setUseAbsDCA(false); - ft.setWeightedFinalPCA(false); - swW.Start(false); - int ncW = ft.process(vctracks[0], vctracks[1], vctracks[2]); // HERE WE FIT THE VERTICES - swW.Stop(); - LOG(debug) << "fit wgh.dist " << iev << " NC: " << ncW << " Chi2: " << (ncW ? ft.getChi2AtPCACandidate(0) : -1); - if (ncW) { - auto minD = checkResults(outStream, treeName3W, ft, vtxGen, genParent, dchdec); - meanDW += minD; - nfoundW++; - } - } - ft.print(); - meanDA /= nfoundA ? nfoundA : 1; - meanDAW /= nfoundAW ? nfoundAW : 1; - meanDW /= nfoundW ? nfoundW : 1; - LOG(info) << "Processed " << NTest << " 3-prong vertices"; - LOG(info) << "3-prongs with abs.dist minization: eff= " << float(nfoundA) / NTest - << " mean.dist to truth: " << meanDA << " CPU time: " << swA.CpuTime(); - LOG(info) << "3-prongs with abs.dist but wghPCA: eff= " << float(nfoundAW) / NTest - << " mean.dist to truth: " << meanDAW << " CPU time: " << swAW.CpuTime(); - LOG(info) << "3-prongs with wgh.dist minization: eff= " << float(nfoundW) / NTest - << " mean.dist to truth: " << meanDW << " CPU time: " << swW.CpuTime(); - BOOST_CHECK(nfoundA > 0.99 * NTest); - BOOST_CHECK(nfoundAW > 0.99 * NTest); - BOOST_CHECK(nfoundW > 0.99 * NTest); - BOOST_CHECK(meanDA < 0.1); - BOOST_CHECK(meanDAW < 0.1); - BOOST_CHECK(meanDW < 0.1); - } - - outStream.Close(); -} - -} // namespace vertexing -} // namespace o2