From 7e57b62becf35f6f51c7d2e749caba5ee2235a23 Mon Sep 17 00:00:00 2001 From: shahoian Date: Tue, 2 Mar 2021 11:23:18 +0100 Subject: [PATCH] Remove non-templated DCAFitter (overridden by DCAFitterN) --- Analysis/Tasks/validation.cxx | 1 - Detectors/Base/CMakeLists.txt | 11 - .../Base/include/DetectorsBase/DCAFitter.h | 370 ------------ Detectors/Base/src/DCAFitter.cxx | 571 ------------------ Detectors/Base/src/DetectorsBaseLinkDef.h | 2 - 5 files changed, 955 deletions(-) delete mode 100644 Detectors/Base/include/DetectorsBase/DCAFitter.h delete mode 100644 Detectors/Base/src/DCAFitter.cxx diff --git a/Analysis/Tasks/validation.cxx b/Analysis/Tasks/validation.cxx index 1a6f6dafdee25..dfe0b84567841 100644 --- a/Analysis/Tasks/validation.cxx +++ b/Analysis/Tasks/validation.cxx @@ -10,7 +10,6 @@ #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "Framework/AnalysisDataModel.h" -#include "DetectorsBase/DCAFitter.h" #include "ReconstructionDataFormats/Track.h" #include diff --git a/Detectors/Base/CMakeLists.txt b/Detectors/Base/CMakeLists.txt index dda6c34f843a1..20d1db9eb08a0 100644 --- a/Detectors/Base/CMakeLists.txt +++ b/Detectors/Base/CMakeLists.txt @@ -16,7 +16,6 @@ o2_add_library(DetectorsBase src/MatLayerCyl.cxx src/MatLayerCylSet.cxx src/Ray.cxx - src/DCAFitter.cxx src/BaseDPLDigitizer.cxx src/CTFCoderBase.cxx PUBLIC_LINK_LIBRARIES FairRoot::Base @@ -36,7 +35,6 @@ o2_add_library(DetectorsBase o2_target_root_dictionary(DetectorsBase HEADERS include/DetectorsBase/Detector.h - include/DetectorsBase/DCAFitter.h include/DetectorsBase/GeometryManager.h include/DetectorsBase/MaterialManager.h include/DetectorsBase/Propagator.h @@ -57,15 +55,6 @@ if(BUILD_SIMULATION) VMCWORKDIR=${CMAKE_BINARY_DIR}/stage/${CMAKE_INSTALL_DATADIR}) endif() -o2_add_test( - DCAFitter - SOURCES test/testDCAFitter.cxx - COMPONENT_NAME DetectorsBase - PUBLIC_LINK_LIBRARIES O2::DetectorsBase - LABELS detectorsbase - ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage - VMCWORKDIR=${CMAKE_BINARY_DIR}/stage/${CMAKE_INSTALL_DATADIR}) - o2_add_test_root_macro(test/buildMatBudLUT.C PUBLIC_LINK_LIBRARIES O2::DetectorsBase LABELS detectorsbase) diff --git a/Detectors/Base/include/DetectorsBase/DCAFitter.h b/Detectors/Base/include/DetectorsBase/DCAFitter.h deleted file mode 100644 index 976bcafde041c..0000000000000 --- a/Detectors/Base/include/DetectorsBase/DCAFitter.h +++ /dev/null @@ -1,370 +0,0 @@ -// Copyright CERN and copyright holders of ALICE O2. This software is -// distributed under the terms of the GNU General Public License v3 (GPL -// Version 3), copied verbatim in the file "COPYING". -// -// See http://alice-o2.web.cern.ch/license for full licensing information. -// -// 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 DCAFitter.h -/// \brief Defintions for DCA fitter class -/// \author ruben.shahoyan@cern.ch - -#ifndef _ALICEO2_DCA_FITTER_ -#define _ALICEO2_DCA_FITTER_ - -#include -#include - -//#define _ADAPT_FOR_ALIROOT_ // to make it compatible with AliRoot AliExternalTrackParam - -#ifdef _ADAPT_FOR_ALIROOT_ -#include "AliExternalTrackParam.h" -typedef AliExternalTrackParam Track; -typedef float ftype_t; // type precision for standard calculation (prefer float) -typedef double dtype_t; // type precision for calculation with risk of round-off errors -#define CONSTRDEF \ - { \ - } -#else -#define CONSTRDEF = default - -#include "ReconstructionDataFormats/Track.h" - -namespace o2 -{ -namespace base -{ -#endif - -class DCAFitter -{ - public: -#ifdef _ADAPT_FOR_ALIROOT_ -#define getSigmaY2 GetSigmaY2 -#define getSigmaZ2 GetSigmaZ2 -#define getSigmaZY GetSigmaZY -#define getAlpha GetAlpha -#define getX GetX -#define getY GetY -#define getZ GetZ -#define getSnp GetSnp -#define getTgl GetTgl -#define getParam GetParameter -#define getCurvature GetC -#define propagateTo PropagateTo -#define propagateParamTo PropagateParamOnlyTo - // -#else - using Track = o2::track::TrackParCov; - using ftype_t = float; // type precision for standard calculation (prefer float) - using dtype_t = double; // type precision for calculation with risk of round-off errors -#endif - - // ---> Auxiliary structs used by DCA finder - - //---------------------------------------------------- - ///< Inverse cov matrix of the point defined by the track - struct TrackCovI { - ftype_t sxx, syy, syz, szz; - - TrackCovI(const Track& trc) { set(trc); } - TrackCovI() CONSTRDEF; - void set(const Track& trc) - { - // we assign Y error to X for DCA calculation of 2 points - // (otherwise for quazi-collinear tracks the X will not be constrained) - ftype_t cyy = trc.getSigmaY2(), czz = trc.getSigmaZ2(), cyz = trc.getSigmaZY(), cxx = cyy; - ftype_t detYZ = cyy * czz - cyz * cyz; - if (detYZ > 0.) { - sxx = 1. / cxx; - syy = czz / detYZ; - syz = -cyz / detYZ; - szz = cyy / detYZ; - } else { - syy = 0.0; // failure - } - } - }; - - //---------------------------------------------------- - ///< Derivative (up to 2) of the TrackParam position over its running param X - struct TrackDeriv2 { - ftype_t dydx, dzdx, d2ydx2, d2zdx2; - - TrackDeriv2() CONSTRDEF; - TrackDeriv2(const Track& trc, ftype_t bz) { set(trc, bz); } - void set(const Track& trc, ftype_t bz) - { - ftype_t snp = trc.getSnp(), csp = TMath::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 - } - }; - - //---------------------------------------------------- - //< precalculated track radius, center, alpha sin,cos and their combinations - struct TrcAuxPar { - dtype_t c, s; // cos ans sin of track alpha - dtype_t cc, cs, ss; // products - dtype_t r, xCen, yCen; // helix radius and center in lab - - TrcAuxPar() CONSTRDEF; - TrcAuxPar(const Track& trc, ftype_t bz) { set(trc, bz); } - - void set(const Track& trc, ftype_t bz) - { - c = TMath::Cos(trc.getAlpha()); - s = TMath::Sin(trc.getAlpha()); - setRCen(trc, bz); - cc = c * c; - ss = s * s; - cs = c * s; - } - - void setRCen(const Track& tr, ftype_t bz); - - void glo2loc(ftype_t vX, ftype_t vY, ftype_t& vXL, ftype_t& vYL) const - { - // rotate XY in global frame to the frame of track with angle A - vXL = vX * c + vY * s; - vYL = -vX * s + vY * c; - } - - void loc2glo(ftype_t vXL, ftype_t vYL, ftype_t& vX, ftype_t& vY) const - { - // rotate XY in local alpha frame to global frame - vX = vXL * c - vYL * s; - vY = vXL * s + vYL * c; - } - void glo2loc(dtype_t vX, dtype_t vY, dtype_t& vXL, dtype_t& vYL) const - { - // rotate XY in global frame to the frame of track with angle A - vXL = vX * c + vY * s; - vYL = -vX * s + vY * c; - } - - void loc2glo(dtype_t vXL, dtype_t vYL, dtype_t& vX, dtype_t& vY) const - { - // rotate XY in local alpha frame to global frame - vX = vXL * c - vYL * s; - vY = vXL * s + vYL * c; - } - }; - - //---------------------------------------------------- - //< coefficients of the track-point contribution to the PCA (Vx,Vy,Vz) to 2 points in lab frame represented via local points coordinates as - //< Vx = mXX0*x0+mXY0*y0+mXZ0*z0 + mXX1*x1+mXY1*y1+mXZ1*z1 - //< Vy = mYX0*x0+mYY0*y0+mYZ0*z0 + mYX1*x1+mYY1*y1+mYZ1*z1 - //< Vz = mZX0*x0+mZY0*y0+mZZ0*z0 + mZX1*x1+mZY1*y1+mZZ1*z1 - //< where {x0,y0,z0} and {x1,y1,z1} are track positions in their local frames - struct TrackCoefVtx { - ftype_t mXX, mXY, mXZ, mYX, mYY, mYZ, mZX, mZY, mZZ; - TrackCoefVtx() CONSTRDEF; - }; - - //---------------------------------------------------- - ///< particular point on track trajectory (in track proper alpha-frame) - struct Triplet { - ftype_t x, y, z; - Triplet(Track& trc) { set(trc); } - Triplet(ftype_t px = 0, ftype_t py = 0, ftype_t pz = 0) : x(px), y(py), z(pz) {} - void set(const Track& trc) - { - x = trc.getX(); - y = trc.getY(); - z = trc.getZ(); - } - }; - - //---------------------------------------------------- - //< crossing coordinates of 2 circles - struct CrossInfo { - ftype_t xDCA[2]; - ftype_t yDCA[2]; - int nDCA; - - CrossInfo() CONSTRDEF; - CrossInfo(const TrcAuxPar& trc0, const TrcAuxPar& trc1) { set(trc0, trc1); } - void set(const TrcAuxPar& trc0, const TrcAuxPar& trc1); - - void notTouchingXY(ftype_t dist, ftype_t xDist, ftype_t yDist, const TrcAuxPar& trcA, ftype_t 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.xCen, trcA.yCen ), xDist = trcB.xCen = trcA.xCen ... - // There are 2 special cases: - // (a) small circle is inside the large one: provide rBSign as -trcB.r - // (b) circle are side by side: provide rBSign as trcB.r - nDCA = 1; - auto t2d = (dist + trcA.r - rBSign) / dist; - xDCA[0] = trcA.xCen + 0.5 * (xDist * t2d); - yDCA[0] = trcA.yCen + 0.5 * (yDist * t2d); - } - }; - - struct Derivatives { - dtype_t dChidx0, dChidx1; // 1st derivatives of chi2 vs tracks local parameters X - dtype_t dChidx0dx0, dChidx1dx1, dChidx0dx1; // 2nd derivatives of chi2 vs tracks local parameters X - }; - - // <--- Auxiliary structs used by DCA finder - - //=============================================================================== - - DCAFitter() CONSTRDEF; - - int getMaxIter() const - { - return mMaxIter; - } - ftype_t getMaxR() const { return TMath::Sqrt(mMaxR2); } - ftype_t getMaxChi2() const { return mMaxChi2; } - ftype_t getMinParamChange() const { return mMinParamChange; } - ftype_t getBz() const { return mBz; } - bool getUseAbsDCA() const { return mUseAbsDCA; } - - void setMaxIter(int n = 20) { mMaxIter = n > 2 ? n : 2; } - void setMaxR(ftype_t r = 200.) { mMaxR2 = r * r; } - void setMaxChi2(ftype_t chi2 = 999.) { mMaxChi2 = chi2; } - void setBz(ftype_t bz) { mBz = bz; } - void setMinParamChange(ftype_t x = 1e-3) { mMinParamChange = x > 1e-4 ? x : 1.e-4; } - void setMinRelChi2Change(ftype_t r = 0.9) { mMinRelChi2Change = r > 0.1 ? r : 999.; } - void setUseAbsDCA(bool v) { mUseAbsDCA = v; } - - void clear() - { - mNCandidates = 0; - mCrossIDCur = 0; - mCrossIDAlt = -1; - } - - DCAFitter(ftype_t bz, ftype_t minRelChiChange = 0.9, ftype_t minXChange = 1e-3, ftype_t maxChi = 999, int n = 20, ftype_t maxR = 200.) - { - setMaxIter(n); - setMaxR(maxR); - setMaxChi2(maxChi); - setMinParamChange(minXChange); - setMinRelChi2Change(minRelChiChange); - setBz(bz); - setUseAbsDCA(false); // by default use weighted DCA definition (much slower) - } - - ///< number of validated V0 candidates (at most 2 are possible) - int getNCandidates() const { return mNCandidates; } - - ///< return PCA candidate (no check for its validity) - const Triplet& getPCACandidate(int cand) const { return mPCA[cand]; } - - ///< return Chi2 at PCA candidate (no check for its validity) - ftype_t getChi2AtPCACandidate(int cand) const { return mChi2[cand]; } - - ///< 1st track params propagated to V0 candidate (no check for the candidate validity) - const Track& getTrack0(int cand) const { return mCandTr0[cand]; } - - ///< 2nd track params propagated to V0 candidate (no check for the candidate validity) - const Track& getTrack1(int cand) const { return mCandTr1[cand]; } - - ///< calculate parameters tracks at PCA - int process(const Track& trc0, const Track& trc1) - { - // find dca of 2 tracks - TrcAuxPar trc0Aux(trc0, mBz), trc1Aux(trc1, mBz); - return process(trc0, trc0Aux, trc1, trc1Aux); - } - - ///< calculate parameters tracks at PCA, using precalculated aux info // = TrcAuxPar(track) // - int process(const Track& trc0, const TrcAuxPar& trc0Aux, - const Track& trc1, const TrcAuxPar& trc1Aux); - - ///< minimizer for abs distance definition of DCA, starting with cached tracks - bool processCandidateDCA(const TrcAuxPar& trc0Aux, const TrcAuxPar& trc1Aux); - - ///< minimizer for weighted distance definition of DCA (chi2), starting with cached tracks - bool processCandidateChi2(const TrcAuxPar& trc0Aux, const TrcAuxPar& trc1Aux); - - ///< minimize w/o preliminary propagation to XY crossing points - int processAsIs(const Track& trc0, const Track& trc1); - - ///< calculate squared distance between 2 tracks - static ftype_t getDistance2(const Track& trc0, const Track& trc1); - - ///< calculate half sum of squared distances between 2 tracks and vertex - static ftype_t getDistance2(ftype_t x, ftype_t y, ftype_t z, const Track& trc0, const Track& trc1); - - protected: - void calcPCACoefs(const TrcAuxPar& trc0Aux, const TrackCovI& trcEI0, - const TrcAuxPar& trc1Aux, const TrackCovI& trcEI1, - TrackCoefVtx& trCFVT0, TrackCoefVtx& trCFVT1) const; - - ///< PCA with weighted DCA definition - void calcPCA(const Track& trc0, const TrackCoefVtx& trCFVT0, - const Track& trc1, const TrackCoefVtx& trCFVT1, - Triplet& v) const; - - ///< PCA with abs DCA definition - void calcPCA(const Track& trc0, const TrcAuxPar& trc0Aux, - const Track& trc1, const TrcAuxPar& trc1Aux, - Triplet& v) const; - - ///< chi2 (weighted distance) - ftype_t calcChi2(const Triplet& pca, - const Track& trc0, const TrcAuxPar& trc0Aux, const TrackCovI& trcEI0, - const Track& trc1, const TrcAuxPar& trc1Aux, const TrackCovI& trcEI1) const; - - ///< DCA (abs distance) - ftype_t calcDCA(const Track& tPnt0, const TrcAuxPar& trc0Aux, const Track& tPnt1, const TrcAuxPar& trc1Aux) const; - ftype_t calcDCA(const Triplet& pca, const Track& trc0, const TrcAuxPar& trc0Aux, const Track& trc1, const TrcAuxPar& trc1Aux) const; - - Triplet calcResid(const Track& trc, const TrcAuxPar& alpCS, const Triplet& vtx) const - { - ftype_t vlX, vlY; // Vertex XY in track local frame - alpCS.glo2loc(vtx.x, vtx.y, vlX, vlY); - return Triplet(trc.getX() - vlX, trc.getY() - vlY, trc.getZ() - vtx.z); - } - - void chi2Deriv(const Track& trc0, const TrackDeriv2& tDer0, const TrcAuxPar& trc0Aux, const TrackCovI& trcEI0, const TrackCoefVtx& trCFVT0, - const Track& trc1, const TrackDeriv2& tDer1, const TrcAuxPar& trc1Aux, const TrackCovI& trcEI1, const TrackCoefVtx& trCFVT1, - Derivatives& deriv) const; - - void DCADeriv(const Track& trc0, const TrackDeriv2& tDer0, const TrcAuxPar& trc0Aux, - const Track& trc1, const TrackDeriv2& tDer1, const TrcAuxPar& trc1Aux, - Derivatives& deriv) const; - - bool closerToAlternative(ftype_t x, ftype_t y) const; - - private: - bool mUseAbsDCA; // ignore track errors (minimize abs DCA to vertex) - int mMaxIter; // max iterations - ftype_t mMaxR2; // max radius to consider - ftype_t mMaxChi2; // max chi2 to accept - ftype_t mMinRelChi2Change; // stop iterations if relative chi2 change is less than requested - ftype_t mMinParamChange; // stop iterations when both X params change by less than this value - - ftype_t mBz; // mag field for simple propagation - - // Working arrays - CrossInfo mCrossings; //! analystical XY crossings (max 2) of the seeds - int mCrossIDCur; //! XY crossing being tested - int mCrossIDAlt; //! XY crossing alternative to the one being tested. Abandon fit if it converges to it - - int mNCandidates; //! number of consdered candidates - Triplet mPCA[2]; //! PCA for 2 possible cases - ftype_t mChi2[2]; //! Chi2 at PCA candidate - TrackCovI mTrcEI0[2], mTrcEI1[2]; //! errors for each track candidate - TrackCoefVtx mTrCFVT0[2], mTrCFVT1[2]; //! coefficients of PCA vs track points for each track - Track mCandTr0[2], mCandTr1[2]; //! Tracks at PCA, max 2 candidates. Note: Errors are at seed XY point - - ClassDefNV(DCAFitter, 1); -}; - -#ifndef _ADAPT_FOR_ALIROOT_ -} -} -#endif - -#endif // _O2_DCA_FITTER_ diff --git a/Detectors/Base/src/DCAFitter.cxx b/Detectors/Base/src/DCAFitter.cxx deleted file mode 100644 index 6fd2cf12200b4..0000000000000 --- a/Detectors/Base/src/DCAFitter.cxx +++ /dev/null @@ -1,571 +0,0 @@ -// Copyright CERN and copyright holders of ALICE O2. This software is -// distributed under the terms of the GNU General Public License v3 (GPL -// Version 3), copied verbatim in the file "COPYING". -// -// See http://alice-o2.web.cern.ch/license for full licensing information. -// -// 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 DCAFitter.cxx -/// \brief Implementations for DCA fitter class - -#include "DetectorsBase/DCAFitter.h" - -#ifndef _ADAPT_FOR_ALIROOT_ -using namespace o2::base; -using ftype_t = DCAFitter::ftype_t; -using dtype_t = DCAFitter::dtype_t; - -#else -#include "DCAFitter.h" -#endif - -//_____________________________________________________________________________________ -void DCAFitter::CrossInfo::set(const TrcAuxPar& trc0, const TrcAuxPar& trc1) -{ - // calculate up to 2 crossings between 2 circles - nDCA = 0; - const auto& trcA = trc0.r > trc1.r ? trc0 : trc1; // designate the largest circle as A - const auto& trcB = trc0.r > trc1.r ? trc1 : trc0; - ftype_t xDist = trcB.xCen - trcA.xCen, yDist = trcB.yCen - trcA.yCen; - ftype_t dist2 = xDist * xDist + yDist * yDist, dist = TMath::Sqrt(dist2), rsum = trcA.r + trcB.r; - if (TMath::Abs(dist) < 1e-12) { - return; // 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) - notTouchingXY(dist, xDist, yDist, trcA, trcB.r); - } else if (dist + trcB.r < trcA.r) { // 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.r); - } 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 (TMath::Abs(xDist) < TMath::Abs(yDist)) { - ftype_t a = (trcA.r * trcA.r - trcB.r * trcB.r + dist2) / (2. * yDist), b = -xDist / yDist, ab = a * b, bb = b * b; - ftype_t det = ab * ab - (1. + bb) * (a * a - trcA.r * trcA.r); - if (det > 0.) { - det = TMath::Sqrt(det); - xDCA[0] = (-ab + det) / (1. + b * b); - yDCA[0] = a + b * xDCA[0] + trcA.yCen; - xDCA[0] += trcA.xCen; - xDCA[1] = (-ab - det) / (1. + b * b); - yDCA[1] = a + b * xDCA[1] + trcA.yCen; - xDCA[1] += trcA.xCen; - 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.r); - } - } else { - ftype_t a = (trcA.r * trcA.r - trcB.r * trcB.r + dist2) / (2. * xDist), b = -yDist / xDist, ab = a * b, bb = b * b; - ftype_t det = ab * ab - (1. + bb) * (a * a - trcA.r * trcA.r); - if (det > 0.) { - det = TMath::Sqrt(det); - yDCA[0] = (-ab + det) / (1. + bb); - xDCA[0] = a + b * yDCA[0] + trcA.xCen; - yDCA[0] += trcA.yCen; - yDCA[1] = (-ab - det) / (1. + bb); - xDCA[1] = a + b * yDCA[1] + trcA.xCen; - yDCA[1] += trcA.yCen; - 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.r); - } - } - } -} - -//_____________________________________________________________________________________ -void DCAFitter::TrcAuxPar::setRCen(const Track& tr, ftype_t bz) -{ - // set track radius and circle coordinates in global frame, no check for bz==0! - r = 1. / tr.getCurvature(bz); - ftype_t sn = tr.getSnp(), cs = TMath::Sqrt((1. - sn) * (1. + sn)); - loc2glo(tr.getX() - sn * r, tr.getY() + cs * r, xCen, yCen); - r = TMath::Abs(r); -} - -//___________________________________________________________________ -int DCAFitter::process(const Track& trc0, const TrcAuxPar& trc0Aux, - const Track& trc1, const TrcAuxPar& trc1Aux) -{ - // find dca of 2 tracks with aux info preclaculated - mNCandidates = 0; - mCrossings.set(trc0Aux, trc1Aux); // find at most 2 candidates of 2 circles crossing - if (!mCrossings.nDCA) { - return 0; // no crossing - } - for (int ic = 0; ic < mCrossings.nDCA; ic++) { - // both XY crossings may eventually converge to the same point. To stop asap this redundant step - // we abandon fit if it appears to be closer to alternative crossing - mCrossIDCur = ic; - mCrossIDAlt = mCrossings.nDCA == 2 ? 1 - ic : -1; - // check if radius is acceptable - if (mCrossings.xDCA[ic] * mCrossings.xDCA[ic] + mCrossings.yDCA[ic] * mCrossings.yDCA[ic] > mMaxR2) { - continue; - } - // find dca starting from proximity of transverse point xv,yv - ftype_t xl, yl; - mCandTr0[mNCandidates] = trc0; - mCandTr1[mNCandidates] = trc1; - trc0Aux.glo2loc(mCrossings.xDCA[ic], mCrossings.yDCA[ic], xl, yl); - if (!mCandTr0[mNCandidates].propagateTo(xl, mBz)) { - continue; - } - trc1Aux.glo2loc(mCrossings.xDCA[ic], mCrossings.yDCA[ic], xl, yl); - if (!mCandTr1[mNCandidates].propagateTo(xl, mBz)) { - continue; - } - if (mUseAbsDCA ? processCandidateDCA(trc0Aux, trc1Aux) : processCandidateChi2(trc0Aux, trc1Aux)) { - mNCandidates++; // candidate validated - } - } - return mNCandidates; -} - -//___________________________________________________________________ -bool DCAFitter::processCandidateChi2(const TrcAuxPar& trc0Aux, const TrcAuxPar& trc1Aux) -{ - // find best chi2 (weighted DCA) of 2 tracks already propagated to their approximate vicinity - mChi2[mNCandidates] = 1e9; - Track &trc0 = mCandTr0[mNCandidates], &trc1 = mCandTr1[mNCandidates]; - TrackCovI &trcEI0 = mTrcEI0[mNCandidates], &trcEI1 = mTrcEI1[mNCandidates]; - // get error matrices at initial point - trcEI0.set(trc0); - trcEI1.set(trc1); - TrackCoefVtx &trCFVT0 = mTrCFVT0[mNCandidates], &trCFVT1 = mTrCFVT1[mNCandidates]; // get coefficients of PCA vs track points - calcPCACoefs(trc0Aux, trcEI0, trc1Aux, trcEI1, trCFVT0, trCFVT1); - TrackDeriv2 tDer0, tDer1; // their derivatives over track param X - Derivatives deriv; // chi2 1st and 2nd derivatives over tracks params - - ftype_t chi2 = 0, chi2Prev = 0; - int iter = 0; - do { - tDer0.set(trc0, mBz); // tracks derivatives - tDer1.set(trc1, mBz); // over their X params - chi2Deriv(trc0, tDer0, trc0Aux, trcEI0, trCFVT0, trc1, tDer1, trc1Aux, trcEI1, trCFVT1, deriv); - - // do Newton-Rapson iteration with corrections = - dchi2/d{x0,x1} * [ d^2chi2/d{x0,x1}^2 ]^-1 - dtype_t detDer2 = deriv.dChidx0dx0 * deriv.dChidx1dx1 - deriv.dChidx0dx1 * deriv.dChidx0dx1; - dtype_t detDer2I = 1. / detDer2; - dtype_t dX0 = -(deriv.dChidx0 * deriv.dChidx1dx1 - deriv.dChidx1 * deriv.dChidx0dx1) * detDer2I; - dtype_t dX1 = -(deriv.dChidx1 * deriv.dChidx0dx0 - deriv.dChidx0 * deriv.dChidx0dx1) * detDer2I; - if (!trc0.propagateParamTo(trc0.getX() + dX0, mBz) || !trc1.propagateParamTo(trc1.getX() + dX1, mBz)) { - return false; - } - Triplet pca; - calcPCA(trc0, trCFVT0, trc1, trCFVT1, pca); - // if there are 2 XY crossings, fit with both of them used as a starting point may converge to - // the same point. To stop asap this redundant step we abandon fit if it appears to be closer - // to alternative crossing - if (mCrossIDAlt >= 0 && closerToAlternative(pca.x, pca.y)) { - return false; - } - - chi2 = calcChi2(pca, trc0, trc0Aux, trcEI0, trc1, trc1Aux, trcEI1); - - if ((TMath::Abs(dX0) < mMinParamChange && TMath::Abs(dX1) < mMinParamChange) || - (iter && chi2 / chi2Prev > mMinRelChi2Change)) { - break; // converged - } - chi2Prev = chi2; - } while (++iter < mMaxIter); - // - if (chi2 < mMaxChi2) { - auto& pca = mPCA[mNCandidates]; - calcPCA(trc0, trCFVT0, trc1, trCFVT1, pca); - mChi2[mNCandidates] = calcChi2(pca, trc0, trc0Aux, trcEI0, trc1, trc1Aux, trcEI1); - return true; - } - return false; -} - -//___________________________________________________________________ -int DCAFitter::processAsIs(const Track& trc0, const Track& trc1) -{ - // find dca of 2 tracks w/o preliminary propagation to XY crossing points - clear(); - TrcAuxPar trc0Aux(trc0, mBz), trc1Aux(trc1, mBz); - mCandTr0[mNCandidates] = trc0; - mCandTr1[mNCandidates] = trc1; - if (mUseAbsDCA ? processCandidateDCA(trc0Aux, trc1Aux) : processCandidateChi2(trc0Aux, trc1Aux)) { - mNCandidates++; // candidate validated - } - return mNCandidates; -} - -//___________________________________________________________________ -bool DCAFitter::processCandidateDCA(const TrcAuxPar& trc0Aux, const TrcAuxPar& trc1Aux) -{ - // find DCA of 2 tracks already propagated to its approximate vicinity W/O APPLYING ANY ERRORS - // i.e. the absolute distance is minimized - Track &trc0 = mCandTr0[mNCandidates], &trc1 = mCandTr1[mNCandidates]; - - TrackDeriv2 tDer0, tDer1; // their derivatives over track param X - Derivatives deriv; // chi2 1st and 2nd derivatives over tracks params - - ftype_t chi2 = 0, chi2Prev = 0; - int iter = 0; - do { - tDer0.set(trc0, mBz); // tracks derivatives - tDer1.set(trc1, mBz); // over their X params - DCADeriv(trc0, tDer0, trc0Aux, trc1, tDer1, trc1Aux, deriv); - - // do Newton-Rapson iteration with corrections = - dchi2/d{x0,x1} * [ d^2chi2/d{x0,x1}^2 ]^-1 - dtype_t detDer2 = deriv.dChidx0dx0 * deriv.dChidx1dx1 - deriv.dChidx0dx1 * deriv.dChidx0dx1; - dtype_t detDer2I = 1. / detDer2; - dtype_t dX0 = -(deriv.dChidx0 * deriv.dChidx1dx1 - deriv.dChidx1 * deriv.dChidx0dx1) * detDer2I; - dtype_t dX1 = -(deriv.dChidx1 * deriv.dChidx0dx0 - deriv.dChidx0 * deriv.dChidx0dx1) * detDer2I; - if (!trc0.propagateParamTo(trc0.getX() + dX0, mBz) || !trc1.propagateParamTo(trc1.getX() + dX1, mBz)) { - return false; - } - - if (mCrossIDAlt >= 0) { - Triplet pca; - calcPCA(trc0, trc0Aux, trc1, trc1Aux, pca); - // if there are 2 XY crossings, fit with both of them used as a starting point may converge to - // the same point. To stop asap this redundant step we abandon fit if it appears to be closer - // to alternative crossing - if (closerToAlternative(pca.x, pca.y)) { - return false; - } - } - chi2 = calcDCA(trc0, trc0Aux, trc1, trc1Aux); - if ((TMath::Abs(dX0) < mMinParamChange && TMath::Abs(dX1) < mMinParamChange) || - (iter && chi2 / chi2Prev > mMinRelChi2Change)) { - break; // converged - } - chi2Prev = chi2; - } while (++iter < mMaxIter); - // - if (chi2 < mMaxChi2) { - auto& pca = mPCA[mNCandidates]; - calcPCA(trc0, trc0Aux, trc1, trc1Aux, pca); - mChi2[mNCandidates] = calcDCA(trc0, trc0Aux, trc1, trc1Aux); - return true; - } - return false; -} - -//___________________________________________________________________ -void DCAFitter::calcPCA(const Track& trc0, const TrackCoefVtx& trCFVT0, - const Track& trc1, const TrackCoefVtx& trCFVT1, - Triplet& pca) const -{ - // calculate the PCA (Vx,Vy,Vz) to 2 points in lab frame represented via local points coordinates - pca.x = trCFVT0.mXX * trc0.getX() + trCFVT0.mXY * trc0.getY() + trCFVT0.mXZ * trc0.getZ() + trCFVT1.mXX * trc1.getX() + trCFVT1.mXY * trc1.getY() + trCFVT1.mXZ * trc1.getZ(); - pca.y = trCFVT0.mYX * trc0.getX() + trCFVT0.mYY * trc0.getY() + trCFVT0.mYZ * trc0.getZ() + trCFVT1.mYX * trc1.getX() + trCFVT1.mYY * trc1.getY() + trCFVT1.mYZ * trc1.getZ(); - pca.z = trCFVT0.mZX * trc0.getX() + trCFVT0.mZY * trc0.getY() + trCFVT0.mZZ * trc0.getZ() + trCFVT1.mZX * trc1.getX() + trCFVT1.mZY * trc1.getY() + trCFVT1.mZZ * trc1.getZ(); -} - -//___________________________________________________________________ -void DCAFitter::calcPCA(const Track& trc0, const TrcAuxPar& trc0Aux, - const Track& trc1, const TrcAuxPar& trc1Aux, - Triplet& pca) const -{ - // calculate the PCA (Vx,Vy,Vz) to 2 points in lab frame represented via local points coordinates - // w/o accounting for the errors of the points - ftype_t xg, yg; - trc0Aux.loc2glo(trc0.getX(), trc0.getY(), xg, yg); - trc1Aux.loc2glo(trc1.getX(), trc1.getY(), pca.x, pca.y); - pca.x += xg; - pca.y += yg; - pca.z = 0.5 * (trc0.getZ() + trc1.getZ()); - pca.x *= 0.5; - pca.y *= 0.5; -} - -//___________________________________________________________________ -void DCAFitter::chi2Deriv(const Track& trc0, const TrackDeriv2& tDer0, const TrcAuxPar& trc0Aux, const TrackCovI& trcEI0, const TrackCoefVtx& trCFVT0, - const Track& trc1, const TrackDeriv2& tDer1, const TrcAuxPar& trc1Aux, const TrackCovI& trcEI1, const TrackCoefVtx& trCFVT1, - Derivatives& deriv) const -{ - // calculate 1st and 2nd derivatives of wighted DCA (chi2) over track parameters X - Triplet vtx; - calcPCA(trc0, trCFVT0, trc1, trCFVT1, vtx); // calculate PCA for current track-points positions - // calculate residuals - auto res0 = calcResid(trc0, trc0Aux, vtx); - auto res1 = calcResid(trc1, trc1Aux, vtx); - - // res0.x,dy0,dz0 = x0 - Xl0, y0 - Yl0, ... , whith x0,y0,z0: track coords in its frame, and Xl0, Yl0, Zl0 : vertex rotated to same frame - - // aux params to minimize multiplications - ftype_t dx0s = res0.x * trcEI0.sxx; - ftype_t dy0sz0t = res0.y * trcEI0.syy + res0.z * trcEI0.syz; - ftype_t dy0tz0s = res0.y * trcEI0.syz + res0.z * trcEI0.szz; - - ftype_t dx1s = res1.x * trcEI1.sxx; - ftype_t dy1sz1t = res1.y * trcEI1.syy + res1.z * trcEI1.syz; - ftype_t dy1tz1s = res1.y * trcEI1.syz + res1.z * trcEI1.szz; - - //------------------------------- - // At the moment keep this in double. TODO: check if ftype_t is ok - dtype_t xx0DtXYXZx0 = trCFVT0.mXX + trCFVT0.mXY * tDer0.dydx + trCFVT0.mXZ * tDer0.dzdx; - dtype_t xx1DtXYXZx1 = trCFVT1.mXX + trCFVT1.mXY * tDer1.dydx + trCFVT1.mXZ * tDer1.dzdx; - dtype_t yz0DtYYYZx0 = trCFVT0.mYX + trCFVT0.mYY * tDer0.dydx + trCFVT0.mYZ * tDer0.dzdx; - dtype_t yz1DtYYYZx1 = trCFVT1.mYX + trCFVT1.mYY * tDer1.dydx + trCFVT1.mYZ * tDer1.dzdx; - - dtype_t DtXYXZx02 = trCFVT0.mXY * tDer0.d2ydx2 + trCFVT0.mXZ * tDer0.d2zdx2; - dtype_t DtYYYZx02 = trCFVT0.mYY * tDer0.d2ydx2 + trCFVT0.mYZ * tDer0.d2zdx2; - dtype_t DtXYXZx12 = trCFVT1.mXY * tDer1.d2ydx2 + trCFVT1.mXZ * tDer1.d2zdx2; - dtype_t DtYYYZx12 = trCFVT1.mYY * tDer1.d2ydx2 + trCFVT1.mYZ * tDer1.d2zdx2; - - dtype_t FDdx0Dx0 = 1. - (trc0Aux.c * xx0DtXYXZx0 + trc0Aux.s * yz0DtYYYZx0); - dtype_t FDdy0Dx0 = trc0Aux.s * xx0DtXYXZx0 - trc0Aux.c * yz0DtYYYZx0 + tDer0.dydx; - dtype_t FDdz0Dx0 = -trCFVT0.mZX - trCFVT0.mZY * tDer0.dydx + tDer0.dzdx * (1. - trCFVT0.mZZ); - - dtype_t FDdx1Dx1 = 1. - (trc1Aux.c * xx1DtXYXZx1 + trc1Aux.s * yz1DtYYYZx1); - dtype_t FDdy1Dx1 = trc1Aux.s * xx1DtXYXZx1 - trc1Aux.c * yz1DtYYYZx1 + tDer1.dydx; - dtype_t FDdz1Dx1 = -trCFVT1.mZX - trCFVT1.mZY * tDer1.dydx + tDer1.dzdx * (1. - trCFVT1.mZZ); - - dtype_t FDdx0Dx1 = -(trc0Aux.c * xx1DtXYXZx1 + trc0Aux.s * yz1DtYYYZx1); - dtype_t FDdy0Dx1 = trc0Aux.s * xx1DtXYXZx1 - trc0Aux.c * yz1DtYYYZx1; - dtype_t FDdz0Dx1 = FDdz1Dx1 - tDer1.dzdx; // -trCFVT1.mZX - trCFVT1.mZY*tDer1.dydx - trCFVT1.mZZ*tDer1.dzdx; - - dtype_t FDdx1Dx0 = -(trc1Aux.c * xx0DtXYXZx0 + trc1Aux.s * yz0DtYYYZx0); - dtype_t FDdy1Dx0 = trc1Aux.s * xx0DtXYXZx0 - trc1Aux.c * yz0DtYYYZx0; - dtype_t FDdz1Dx0 = FDdz0Dx0 - tDer0.dzdx; // -trCFVT0.mZX - trCFVT0.mZY*tDer0.dydx - trCFVT0.mZZ*tDer0.dzdx; - - dtype_t FDdx0Dx0x0 = -(trc0Aux.c * DtXYXZx02 + trc0Aux.s * DtYYYZx02); - dtype_t FDdy0Dx0x0 = tDer0.d2ydx2 + trc0Aux.s * DtXYXZx02 - trc0Aux.c * DtYYYZx02; - dtype_t FDdz0Dx0x0 = -trCFVT0.mZY * tDer0.d2ydx2 + tDer0.d2zdx2 * (1. - trCFVT0.mZZ); - - dtype_t FDdx1Dx1x1 = -(trc1Aux.c * DtXYXZx12 + trc1Aux.s * DtYYYZx12); - dtype_t FDdy1Dx1x1 = tDer1.d2ydx2 + trc1Aux.s * DtXYXZx12 - trc1Aux.c * DtYYYZx12; - dtype_t FDdz1Dx1x1 = -trCFVT1.mZY * tDer1.d2ydx2 + tDer1.d2zdx2 * (1. - trCFVT1.mZZ); - - dtype_t FDdx0Dx1x1 = -(trc0Aux.c * DtXYXZx12 + trc0Aux.s * DtYYYZx12); - dtype_t FDdx1Dx0x0 = -(trc1Aux.c * DtXYXZx02 + trc1Aux.s * DtYYYZx02); - - dtype_t FDdy1Dx0x0 = trc1Aux.s * DtXYXZx02 - trc1Aux.c * DtYYYZx02; - dtype_t FDdy0Dx1x1 = trc0Aux.s * DtXYXZx12 - trc0Aux.c * DtYYYZx12; - - dtype_t FDdz0Dx1x1 = FDdz1Dx1x1 - tDer1.d2zdx2; // -(trCFVT1.mZY*tDer1.d2ydx2 + trCFVT1.mZZ*tDer1.d2zdx2); - dtype_t FDdz1Dx0x0 = FDdz0Dx0x0 - tDer0.d2zdx2; // -(trCFVT0.mZY*tDer0.d2ydx2 + trCFVT0.mZZ*tDer0.d2zdx2); - - dtype_t FD00YYYZ = FDdy0Dx0 * trcEI0.syy + FDdz0Dx0 * trcEI0.syz; - dtype_t FD11YYYZ = FDdy1Dx1 * trcEI1.syy + FDdz1Dx1 * trcEI1.syz; - dtype_t FD10YYYZ = FDdy1Dx0 * trcEI1.syy + FDdz1Dx0 * trcEI1.syz; - dtype_t FD01YYYZ = FDdy0Dx1 * trcEI0.syy + FDdz0Dx1 * trcEI0.syz; - dtype_t FD00YZZZ = FDdy0Dx0 * trcEI0.syz + FDdz0Dx0 * trcEI0.szz; - dtype_t FD11YZZZ = FDdy1Dx1 * trcEI1.syz + FDdz1Dx1 * trcEI1.szz; - - dtype_t FD10YZZZ = FDdy1Dx0 * trcEI1.syz + FDdz1Dx0 * trcEI1.szz; - dtype_t FD01YZZZ = FDdy0Dx1 * trcEI0.syz + FDdz0Dx1 * trcEI0.szz; - - // 1st derivatives over track params x - deriv.dChidx0 = dx0s * FDdx0Dx0 + dx1s * FDdx1Dx0 + dy0sz0t * FDdy0Dx0 + dy1sz1t * FDdy1Dx0 + dy0tz0s * FDdz0Dx0 + dy1tz1s * FDdz1Dx0; - deriv.dChidx1 = dx1s * FDdx1Dx1 + dx0s * FDdx0Dx1 + dy1sz1t * FDdy1Dx1 + dy0sz0t * FDdy0Dx1 + dy1tz1s * FDdz1Dx1 + dy0tz0s * FDdz0Dx1; - - // 2nd derivative over track params x - deriv.dChidx0dx0 = dx0s * FDdx0Dx0x0 + dx1s * FDdx1Dx0x0 + FDdy0Dx0x0 * dy0sz0t + FDdy1Dx0x0 * dy1sz1t + FDdz0Dx0x0 * dy0tz0s + FDdz1Dx0x0 * dy1tz1s + - FDdx0Dx0 * FDdx0Dx0 * trcEI0.sxx + FDdx1Dx0 * FDdx1Dx0 * trcEI1.sxx + - FDdy0Dx0 * FD00YYYZ + FDdy1Dx0 * FD10YYYZ + FDdz0Dx0 * FD00YZZZ + FDdz1Dx0 * FD10YZZZ; - - deriv.dChidx1dx1 = dx1s * FDdx1Dx1x1 + dx0s * FDdx0Dx1x1 + FDdy1Dx1x1 * dy1sz1t + FDdy0Dx1x1 * dy0sz0t + FDdz1Dx1x1 * dy1tz1s + FDdz0Dx1x1 * dy0tz0s + - FDdx1Dx1 * FDdx1Dx1 * trcEI1.sxx + FDdx0Dx1 * FDdx0Dx1 * trcEI0.sxx + - FDdy1Dx1 * FD11YYYZ + FDdy0Dx1 * FD01YYYZ + FDdz1Dx1 * FD11YZZZ + FDdz0Dx1 * FD01YZZZ; - - // N.B.: cross-derivatice - // FDdx0Dx0x1 -> 0, FDdy0Dx0x1 -> 0, FDdz0Dx0x1 -> 0, FDdx1Dx1x0 -> 0; - // FDdy1Dx1x0 -> 0, FDdz1Dx1x0 -> 0, FDdx0Dx1x0 -> 0, FDdy0Dx1x0 -> 0; - // FDdz0Dx1x0 -> 0, FDdx1Dx0x1 -> 0, FDdy1Dx0x1 -> 0, FDdz1Dx0x1 -> 0; - - deriv.dChidx0dx1 = - // this part is = 0 due to the N.B. - // dx0s*FDdx0Dx0x1 + dx1s*FDdx1Dx0x1 + dy0s*FDdy0Dx0x1 + dy1s*FDdy1Dx0x1 + dz0s*FDdz0Dx0x1 + dz1s*FDdz1Dx0x1 + - // dy0t*FDdz0Dx0x1 + dy1t*FDdz1Dx0x1 +dz0t*FDdy0Dx0x1 + dz1t*FDdy1Dx0x1 + - FDdx0Dx0 * FDdx0Dx1 * trcEI0.sxx + FDdx1Dx0 * FDdx1Dx1 * trcEI1.sxx + - FDdy0Dx0 * (FDdy0Dx1 * trcEI0.syy + FDdz0Dx1 * trcEI0.syz) + - FDdy1Dx0 * (FDdy1Dx1 * trcEI1.syy + FDdz1Dx1 * trcEI1.syz) + - FDdz0Dx0 * (FDdy0Dx1 * trcEI0.syz + FDdz0Dx1 * trcEI0.szz) + - FDdz1Dx0 * (FDdy1Dx1 * trcEI1.syz + FDdz1Dx1 * trcEI1.szz); - // -} - -//___________________________________________________________________ -void DCAFitter::DCADeriv(const Track& trc0, const TrackDeriv2& tDer0, const TrcAuxPar& trc0Aux, - const Track& trc1, const TrackDeriv2& tDer1, const TrcAuxPar& trc1Aux, - Derivatives& deriv) const -{ - // DCA derivative calculation with Chi2 defined as the absolule distance (no errors applied) - ftype_t cosDA = trc0Aux.c * trc1Aux.c + trc0Aux.s * trc1Aux.s; // cos(A0-A1) - ftype_t sinDA = trc0Aux.s * trc1Aux.c - trc0Aux.c * trc1Aux.s; // sin(A0-A1) - ftype_t dx01 = trc0.getX() - (trc1.getX() * cosDA + trc1.getY() * sinDA); - ftype_t dx10 = trc1.getX() - (trc0.getX() * cosDA - trc0.getY() * sinDA); - ftype_t dy01 = trc0.getY() - (trc1.getY() * cosDA - trc1.getX() * sinDA); - ftype_t dy10 = trc1.getY() - (trc0.getY() * cosDA + trc0.getX() * sinDA); - ftype_t dz = trc0.getZ() - trc1.getZ(); - - // 1st derivatives over track params x - deriv.dChidx0 = 0.5 * (tDer0.dydx * dy01 + dx01 + dz * tDer0.dzdx); - deriv.dChidx1 = 0.5 * (tDer1.dydx * dy10 + dx10 - dz * tDer1.dzdx); - - // 2nd derivative over track params x - deriv.dChidx0dx0 = 0.5 * (1. + tDer0.dydx * tDer0.dydx + tDer0.dzdx * tDer0.dzdx + dz * tDer0.d2zdx2 + tDer0.d2ydx2 * dy01); - deriv.dChidx1dx1 = 0.5 * (1. + tDer1.dydx * tDer1.dydx + tDer1.dzdx * tDer1.dzdx - dz * tDer1.d2zdx2 + tDer1.d2ydx2 * dy10); - deriv.dChidx0dx1 = 0.5 * (sinDA * (tDer0.dydx - tDer1.dydx) - cosDA * (1. + tDer0.dydx * tDer1.dydx) - tDer0.dzdx * tDer1.dzdx); -} - -//___________________________________________________________________ -void DCAFitter::calcPCACoefs(const TrcAuxPar& trc0Aux, const TrackCovI& trcEI0, - const TrcAuxPar& trc1Aux, const TrackCovI& trcEI1, - TrackCoefVtx& trCFVT0, TrackCoefVtx& trCFVT1) const -{ - // calculate coefficients of the PCA (Vx,Vy,Vz) to 2 points in lab frame represented via local points coordinates as - // Vx = mXX0*x0+mXY0*y0+mXZ0*z0 + mXX1*x1+mXY1*y1+mXZ1*z1 - // Vy = mYX0*x0+mYY0*y0+mYZ0*z0 + mYX1*x1+mYY1*y1+mYZ1*z1 - // Vz = mZX0*x0+mZY0*y0+mZZ0*z0 + mZX1*x1+mZY1*y1+mZZ1*z1 - // where {x0,y0,z0} and {x1,y1,z1} are track positions in their local frames - // - // we find the PCA of 2 tracks poins weighted by their errors, i.e. minimizing - // chi2 = .... - // these are the coefficients of dChi2/d{Vx,Vy,vZ} = 0 - - // At the moment keep this in dtype_t. TODO: check if ftype_t is ok - dtype_t axx = trc0Aux.cc * trcEI0.sxx + trc1Aux.cc * trcEI1.sxx + trc0Aux.ss * trcEI0.syy + trc1Aux.ss * trcEI1.syy; - dtype_t axy = trc0Aux.cs * (trcEI0.sxx - trcEI0.syy) + trc1Aux.cs * (trcEI1.sxx - trcEI1.syy); - dtype_t axz = -(trc0Aux.s * trcEI0.syz + trc1Aux.s * trcEI1.syz); - dtype_t ayy = trc0Aux.ss * trcEI0.sxx + trc1Aux.ss * trcEI1.sxx + trc0Aux.cc * trcEI0.syy + trc1Aux.cc * trcEI1.syy; // = (trcEI0.sxx + trcEI1.sxx + trcEI0.syy + trcEI1.syy) - axx - dtype_t ayz = trc0Aux.c * trcEI0.syz + trc1Aux.c * trcEI1.syz; - dtype_t azz = trcEI0.szz + trcEI1.szz; - // - // define some aux variables - dtype_t axxyy = axx * ayy, axxzz = axx * azz, axxyz = axx * ayz, - axyxy = axy * axy, axyxz = axy * axz, axyyz = axy * ayz, axyzz = axy * azz, - axzxz = axz * axz, axzyy = axz * ayy, axzyz = axz * ayz, - ayzyz = ayz * ayz, ayyzz = ayy * azz; - dtype_t dAxxyyAxyxy = axxyy - axyxy, dAxyyzAxzyy = axyyz - axzyy, dAxyxzAxxyz = axyxz - axxyz; - dtype_t dAxzyzAxyzz = axzyz - axyzz, dAyyzzAyzyz = ayyzz - ayzyz, dAxxzzAxzxz = axxzz - axzxz; - dtype_t det = -dAxyyzAxzyy * axz + dAxyxzAxxyz * ayz + dAxxyyAxyxy * azz, detI = 1. / det; - - dtype_t dfxPCS0 = dAyyzzAyzyz * trc0Aux.c + dAxzyzAxyzz * trc0Aux.s, dfxQCS0 = dAxzyzAxyzz * trc0Aux.c - dAyyzzAyzyz * trc0Aux.s; - dtype_t dfyPCS0 = dAxzyzAxyzz * trc0Aux.c + dAxxzzAxzxz * trc0Aux.s, dfyQCS0 = dAxxzzAxzxz * trc0Aux.c - dAxzyzAxyzz * trc0Aux.s; - dtype_t dfzPCS0 = dAxyyzAxzyy * trc0Aux.c + dAxyxzAxxyz * trc0Aux.s, dfzQCS0 = dAxyxzAxxyz * trc0Aux.c - dAxyyzAxzyy * trc0Aux.s; - - dtype_t dfxPCS1 = dAyyzzAyzyz * trc1Aux.c + dAxzyzAxyzz * trc1Aux.s, dfxQCS1 = dAxzyzAxyzz * trc1Aux.c - dAyyzzAyzyz * trc1Aux.s; - dtype_t dfyPCS1 = dAxzyzAxyzz * trc1Aux.c + dAxxzzAxzxz * trc1Aux.s, dfyQCS1 = dAxxzzAxzxz * trc1Aux.c - dAxzyzAxyzz * trc1Aux.s; - dtype_t dfzPCS1 = dAxyyzAxzyy * trc1Aux.c + dAxyxzAxxyz * trc1Aux.s, dfzQCS1 = dAxyxzAxxyz * trc1Aux.c - dAxyyzAxzyy * trc1Aux.s; - // - trCFVT0.mXX = detI * (dfxPCS0 * trcEI0.sxx); - trCFVT0.mXY = detI * (dfxQCS0 * trcEI0.syy + dAxyyzAxzyy * trcEI0.syz); - trCFVT0.mXZ = detI * (dfxQCS0 * trcEI0.syz + dAxyyzAxzyy * trcEI0.szz); - - trCFVT0.mYX = detI * (dfyPCS0 * trcEI0.sxx); - trCFVT0.mYY = detI * (dfyQCS0 * trcEI0.syy + dAxyxzAxxyz * trcEI0.syz); - trCFVT0.mYZ = detI * (dfyQCS0 * trcEI0.syz + dAxyxzAxxyz * trcEI0.szz); - - trCFVT0.mZX = detI * (dfzPCS0 * trcEI0.sxx); - trCFVT0.mZY = detI * (dfzQCS0 * trcEI0.syy + dAxxyyAxyxy * trcEI0.syz); - trCFVT0.mZZ = detI * (dfzQCS0 * trcEI0.syz + dAxxyyAxyxy * trcEI0.szz); - - trCFVT1.mXX = detI * (dfxPCS1 * trcEI1.sxx); - trCFVT1.mXY = detI * (dfxQCS1 * trcEI1.syy + dAxyyzAxzyy * trcEI1.syz); - trCFVT1.mXZ = detI * (dfxQCS1 * trcEI1.syz + dAxyyzAxzyy * trcEI1.szz); - - trCFVT1.mYX = detI * (dfyPCS1 * trcEI1.sxx); - trCFVT1.mYY = detI * (dfyQCS1 * trcEI1.syy + dAxyxzAxxyz * trcEI1.syz); - trCFVT1.mYZ = detI * (dfyQCS1 * trcEI1.syz + dAxyxzAxxyz * trcEI1.szz); - - trCFVT1.mZX = detI * (dfzPCS1 * trcEI1.sxx); - trCFVT1.mZY = detI * (dfzQCS1 * trcEI1.syy + dAxxyyAxyxy * trcEI1.syz); - trCFVT1.mZZ = detI * (dfzQCS1 * trcEI1.syz + dAxxyyAxyxy * trcEI1.szz); -} - -//___________________________________________________________________ -ftype_t DCAFitter::calcChi2(const Triplet& pca, - const Track& trc0, const TrcAuxPar& trc0Aux, const TrackCovI& trcEI0, - const Track& trc1, const TrcAuxPar& trc1Aux, const TrackCovI& trcEI1) const -{ - ftype_t chi2 = 0; - ftype_t xl, yl, dx, dy, dz; - trc0Aux.glo2loc(pca.x, pca.y, xl, yl); - dx = trc0.getX() - xl; - dy = trc0.getY() - yl; - dz = trc0.getZ() - pca.z; - chi2 += dx * dx * trcEI0.sxx + dy * dy * trcEI0.syy + dz * dz * trcEI0.szz + 2. * dy * dz * trcEI0.syz; - trc1Aux.glo2loc(pca.x, pca.y, xl, yl); - dx = trc1.getX() - xl; - dy = trc1.getY() - yl; - dz = trc1.getZ() - pca.z; - chi2 += dx * dx * trcEI1.sxx + dy * dy * trcEI1.syy + dz * dz * trcEI1.szz + 2. * dy * dz * trcEI1.syz; - return chi2; -} - -//___________________________________________________________________ -ftype_t DCAFitter::calcDCA(const Triplet& pca, - const Track& trc0, const TrcAuxPar& trc0Aux, const Track& trc1, const TrcAuxPar& trc1Aux) const -{ - // calculate distance (non-weighted) of closest approach of 2 points in their local frame - // (long way, see alternative getDCA w/o explicit vertex calculation) - ftype_t chi2 = 0; - ftype_t xl, yl, dx, dy, dz; - trc0Aux.glo2loc(pca.x, pca.y, xl, yl); - dx = trc0.getX() - xl; - dy = trc0.getY() - yl; - dz = trc0.getZ() - pca.z; - chi2 += dx * dx + dy * dy + dz * dz; - trc1Aux.glo2loc(pca.x, pca.y, xl, yl); - dx = trc1.getX() - xl; - dy = trc1.getY() - yl; - dz = trc1.getZ() - pca.z; - chi2 += dz * dz + dy * dy + dz * dz; - return chi2; -} - -//___________________________________________________________________ -ftype_t DCAFitter::calcDCA(const Track& trc0, const TrcAuxPar& trc0Aux, const Track& trc1, const TrcAuxPar& trc1Aux) const -{ - // calculate distance (non-weighted) of closest approach of 2 points in their local frame - ftype_t chi2 = 0; - dtype_t cosDA = trc0Aux.c * trc1Aux.c + trc0Aux.s * trc1Aux.s; // cos(A0-A1) - dtype_t sinDA = trc0Aux.s * trc1Aux.c - trc0Aux.c * trc1Aux.s; // sin(A0-A1) - dtype_t dx = trc0.getX() - trc1.getX(), dy = trc0.getY() - trc1.getY(), dz = trc0.getZ() - trc1.getZ(); - chi2 = 0.5 * (dx * dx + dy * dy + dz * dz) + (1. - cosDA) * (trc0.getX() * trc1.getX() + trc0.getY() * trc1.getY()) + - sinDA * (trc0.getY() * trc1.getX() - trc1.getY() * trc0.getX()); - - return chi2; -} - -//___________________________________________________________________ -bool DCAFitter::closerToAlternative(ftype_t x, ftype_t y) const -{ - // check if the point x,y is closer to the seeding XY point being tested or to alternative see (if any) - ftype_t dxCur = x - mCrossings.xDCA[mCrossIDCur], dyCur = y - mCrossings.yDCA[mCrossIDCur]; - ftype_t dxAlt = x - mCrossings.xDCA[mCrossIDAlt], dyAlt = y - mCrossings.yDCA[mCrossIDAlt]; - return dxCur * dxCur + dyCur * dyCur > dxAlt * dxAlt + dyAlt * dyAlt; -} - -//___________________________________________________________________ -ftype_t DCAFitter::getDistance2(const Track& trc0, const Track& trc1) -{ - // calculate un-weighted distance^2 between 2 tracks - dtype_t cosalp0 = TMath::Cos(trc0.getAlpha()), sinalp0 = TMath::Sin(trc0.getAlpha()); - dtype_t cosalp1 = TMath::Cos(trc1.getAlpha()), sinalp1 = TMath::Sin(trc1.getAlpha()); - dtype_t x0 = trc0.getX() * cosalp0 - trc0.getY() * sinalp0; - dtype_t y0 = trc0.getX() * sinalp0 + trc0.getY() * cosalp0; - dtype_t x1 = trc1.getX() * cosalp1 - trc1.getY() * sinalp1; - dtype_t y1 = trc1.getX() * sinalp1 + trc1.getY() * cosalp1; - - dtype_t dx = x0 - x1, dy = y0 - y1, dz = trc0.getZ() - trc1.getZ(); - return dx * dx + dy * dy + dz * dz; -} - -//___________________________________________________________________ -ftype_t DCAFitter::getDistance2(ftype_t x, ftype_t y, ftype_t z, const Track& trc0, const Track& trc1) -{ - // calculate un-weighted 1/2 distance^2 between 2 tracks and vertex coordinates - dtype_t cosalp0 = TMath::Cos(trc0.getAlpha()), sinalp0 = TMath::Sin(trc0.getAlpha()); - dtype_t cosalp1 = TMath::Cos(trc1.getAlpha()), sinalp1 = TMath::Sin(trc1.getAlpha()); - dtype_t x0 = trc0.getX() * cosalp0 - trc0.getY() * sinalp0; - dtype_t y0 = trc0.getX() * sinalp0 + trc0.getY() * cosalp0; - dtype_t x1 = trc1.getX() * cosalp1 - trc1.getY() * sinalp1; - dtype_t y1 = trc1.getX() * sinalp1 + trc1.getY() * cosalp1; - - dtype_t dx0 = x0 - x, dy0 = y0 - y, dz0 = trc0.getZ() - z; - dtype_t dx1 = x1 - x, dy1 = y1 - y, dz1 = trc1.getZ() - z; - return 0.5 * (dx0 * dx0 + dy0 * dy0 + dz0 * dz0 + dx1 * dx1 + dy1 * dy1 + dz1 * dz1); -} diff --git a/Detectors/Base/src/DetectorsBaseLinkDef.h b/Detectors/Base/src/DetectorsBaseLinkDef.h index 5e7edcfdc77f3..735d73a64d71a 100644 --- a/Detectors/Base/src/DetectorsBaseLinkDef.h +++ b/Detectors/Base/src/DetectorsBaseLinkDef.h @@ -27,8 +27,6 @@ #pragma link C++ class o2::base::MatLayerCyl + ; #pragma link C++ class o2::base::MatLayerCylSet + ; -#pragma link C++ class o2::base::DCAFitter + ; - #pragma link C++ class o2::ctf::CTFCoderBase + ; #endif