ITS/AliITSTPArrayFit.h

   1 /* Copyright(c) 2009-2011, ALICE Experiment at CERN, All rights reserved. *
   2  * See cxx source for full Copyright notice                               */
   3
   4 /* $Id$ */
   5
   6
   7 #ifndef ALIITSTPARRAYFIT_H
   8 #define ALIITSTPARRAYFIT_H
   9
  10 ///////////////////////////////////////////////////////////////////////////////////////////////
  11 //                                                                                           //
  12 // The line is defined by equations (1)                                                      //
  13 // a0*z+a1*x-a0*a1=0 and                                                                     //
  14 // b0*z+b1*y-b0*b1=0                                                                         //
  15 // where x,y,z are NOT the lab axes but z is the lab axis along which the track              //
  16 // has the largest lever arm and x,y are the remaining 2 axis in                             //
  17 // the order of fgkAxisID[z][0], fgkAxisID[z][1]                                             //
  18 // The parameters are fParams[kA0,kB0,kA1,kB1] and the axis chosen as the independent        //
  19 // var. is fParAxis (i.e. if fParAxis==kZ, then a0=ax,b0=bx, a1=ay,b1=by)                    //
  20 //                                                                                           //
  21 //                                                                                           //
  22 // The helix is defined by the equations (2)                                                 //
  23 // X(t) = (dr+R)*cos(phi0) - (R+sum{dRi})*cos(t+phi0) + sum{dRi*cos(phi0+ti)}                //
  24 // Y(t) = (dr+R)*sin(phi0) - (R+sum{dRi})*sin(t+phi0) + sum{dRi*sin(phi0+ti)}                //
  25 // Z(t) = dz - (R+sum{dRi})*t*tg(dip) + sum{dRi*ti}*tg(dip)                                  //
  26 // where dRi is the change of the radius due to the ELoss at parameter ti                    //
  27 //                                                                                           //
  28 // Author: ruben.shahoyan@cern.ch                                                            //
  29 //                                                                                           //
  30 ///////////////////////////////////////////////////////////////////////////////////////////////
  31
  32
  33 #include <TObject.h>
  34 #include <TMath.h>
  35 #include <AliTrackPointArray.h>
  36 class AliSymMatrix;
  37 class AliLog;
  38 class AliParamSolver;
  39
  40
  41 class AliITSTPArrayFit : public TObject
  42 {
  43  public:
  44   enum {kFitDoneBit=BIT(14),kCovInvBit=BIT(15),
  45         kCosmicsBit=BIT(16),kELossBit=BIT(17),
  46         kIgnoreCovBit=BIT(18),
  47         kMask=BIT(24)-1};
  48   enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5,kScl=6,kNCov};
  49   enum {kA0,kB0,kA1,kB1};                // line params
  50   enum {kD0,kPhi0,kR0,kDZ,kDip};         // helix params
  51   enum {kX,kY,kZ};
  52   enum {kMaxParam=6,kMaxParamSq = kMaxParam*(kMaxParam+1)/2};
  53   enum {kLrBeamPime, kLrSPD1,kLrSPD2, kLrShield1, kLrSDD1,kLrSDD2, kLrShield2, kLrSSD1,kLrSSD2,kMaxLrITS};
  54   //
  55  public:
  56   AliITSTPArrayFit();
  57   AliITSTPArrayFit(Int_t npoints);
  58   AliITSTPArrayFit(const AliITSTPArrayFit &fit);
  59   AliITSTPArrayFit& operator= (const AliITSTPArrayFit& src);
  60   virtual ~AliITSTPArrayFit();
  61   //
  62   void          AttachPoints(const AliTrackPointArray* points, Int_t pfirst=-1,Int_t plast=-1);
  63   Bool_t        SetFirstLast(Int_t pfirst=-1,Int_t plast=-1);
  64   AliTrackPointArray* GetPoints()                           const {return (AliTrackPointArray*)fkPoints;}
  65   //
  66   void          SetBz(Double_t bz)                                {fBz = bz;}
  67   Double_t      GetBz()                                     const {return fBz;}
  68   Bool_t        IsHelix()                                   const {return fParAxis<0;}
  69   Bool_t        IsFieldON()                                 const {return TMath::Abs(fBz)>1e-5;}
  70   Bool_t        IsTypeCosmics()                             const {return TestBit(kCosmicsBit);}
  71   Bool_t        IsTypeCollision()                           const {return !IsTypeCosmics();}
  72   Int_t         GetCharge()                                 const {return fCharge;}
  73   Int_t         GetSignQB()                                 const {return fBz<0 ? -fCharge:fCharge;}
  74   void          GetResiduals(Double_t *res, Int_t ipnt)     const;
  75   void          GetResiduals(Double_t *resPCA, const Double_t* xyz, const Double_t* covI=0, Double_t sclCovI=-1)  const;
  76   Double_t      GetPosition( Double_t *xyzPCA, const Double_t* xyz, const Double_t* covI=0, Double_t sclCovI=-1)  const;
  77   Double_t      GetPosition( Double_t *xyzPCA, const AliTrackPoint *pntCovInv,Bool_t useErr=kFALSE) const;
  78   void          GetResiduals(Double_t *xyzPCA, const AliTrackPoint *pntCovInv,Bool_t useErr=kFALSE) const;
  79   void          GetPosition(Double_t *xyz, Double_t t)      const;
  80   void          GetPosition(Double_t *xyz, Int_t pnt)       const;
  81   void          GetDirCos(Double_t *dircos, Double_t t)     const;
  82   Double_t      GetPCA2PlaneInfo(Double_t *xyz, Double_t *dir=0, Int_t axis=kY, Double_t axval=0) const;
  83   void          GetT0Info(Double_t *xyz, Double_t *dir=0)   const;
  84   Double_t      CalcChi2NDF()                               const;
  85   Double_t      GetChi2NDF()                                const {return fChi2NDF;}
  86   Double_t      GetParPCA(const double *xyz, const double *covI=0, Double_t sclCovI=-1)        const;
  87   Double_t      CalcParPCA(Int_t ipnt)                      const;
  88   Bool_t        CalcErrorMatrix();
  89   //
  90   void          GetDResDParamsLine (Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0/*,Double_t sclCovI=-1*/) const;
  91   void          GetDResDParamsLine (Double_t *dXYZdP, Int_t ipnt) const;
  92   void          GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0, Double_t sclCovI=-1);
  93   void          GetDResDParams(Double_t *dXYZdP, Int_t ipnt);
  94   //
  95   void          GetDResDPosLine (Double_t *dXYZdP,/*const Double_t *xyz,*/ const Double_t *covI=0/*,Double_t sclCovI=-1*/) const;
  96   void          GetDResDPosLine (Double_t *dXYZdP, Int_t ipnt) const;
  97   void          GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0, Double_t sclCovI=-1) const;
  98   void          GetDResDPos(Double_t *dXYZdP, Int_t ipnt);
  99   //
 100   Double_t*     GetPoint(int ip)                            const;
 101   Bool_t        Converged()                                 const {return fIter<fMaxIter;}
 102   //
 103   Double_t      Fit(Int_t extQ=0, Double_t extPT=-1,Double_t extPTerr=0);
 104   Double_t      FitLine();
 105   Double_t      FitHelix(Int_t extQ=0, Double_t extPT=-1,Double_t extPTerr=0);
 106   Bool_t        FitLineCrude();
 107   Bool_t        FitHelixCrude(Int_t imposedQ=0);
 108   //
 109   Int_t         GetParAxis()                                const {return fParAxis;}
 110   Int_t         GetAxID(Int_t id)                           const {return fkAxID  ? fkAxID[id] : -1;}
 111   Int_t         GetAxCID(Int_t id)                          const {return fkAxCID ? fkAxCID[id] : -1;}
 112   Int_t         GetFirst()                                  const {return fPntFirst;}
 113   Int_t         GetLast()                                   const {return fPntLast;}
 114   //
 115   Int_t         GetNParams()                                const {return IsFieldON() ? 5:4;}
 116   Bool_t        InvertPointsCovMat();
 117   //
 118   Int_t*        GetElsId() const {return fElsId;}
 119   Double_t*     GetElsDR() const {return fElsDR;}
 120   //
 121   Double_t      GetCovIScale(Int_t ip)                      const {return ip<fNPBooked ? fCovI[ip*kNCov+kScl] : -1.0;}
 122   Double_t*     GetCovI(Int_t ip)                           const {return fCovI + kNCov*ip;}
 123   Double_t*     GetCovI()                                   const {return fCovI;}
 124   Double_t*     GetParams()                                 const {return (Double_t*)&fParams[0];}
 125   Double_t      GetParam(Int_t ip)                          const {return fParams[ip];}
 126   Double_t*     GetTs()                                     const {return (Double_t*)fCurT;}
 127   Double_t      GetT(Int_t ip)                              const {return fCurT[ip];}
 128   Double_t      GetLineOffset(Int_t axis)                   const;
 129   Double_t      GetLineSlope(Int_t axis)                    const;
 130   //
 131   Bool_t        IsSwitched2Line()                           const {return fSwitch2Line;}
 132   Bool_t        IsELossON()                                 const {return TestBit(kELossBit)&&IsFieldON();}
 133   Bool_t        IsFitDone()                                 const {return TestBit(kFitDoneBit);}
 134   Bool_t        IsCovInv()                                  const {return TestBit(kCovInvBit);}
 135   Bool_t        IsCovIgnored()                              const {return TestBit(kIgnoreCovBit);}
 136   Int_t         GetMaxIterations()                          const {return fMaxIter;}
 137   Int_t         GetNIterations()                            const {return fIter;}
 138   Double_t      GetEps()                                    const {return fEps;}
 139   Double_t      GetMass()                                   const {return fMass;}
 140   //
 141   Double_t      GetPt()                                     const;
 142   Double_t      GetP()                                      const;
 143
 144   //
 145   void          Switch2Line(Bool_t v=kTRUE)                       {fSwitch2Line = v;}
 146   void          SetMaxRforHelix(Double_t r)                       {fMaxRforHelix = r>0 ? r : 1e9;}
 147   void          SetCharge(Int_t q=1)                              {fCharge = q<0 ? -1:1;}
 148   void          SetELossON(Bool_t v=kTRUE)                        {SetBit(kELossBit,v);}
 149   void          SetTypeCosmics(Bool_t v=kTRUE)                    {SetBit(kCosmicsBit,v);}
 150   void          SetTypeCollision(Bool_t v=kTRUE)                  {SetTypeCosmics(!v);}
 151   void          SetFitDone(Bool_t v=kTRUE)                        {SetBit(kFitDoneBit,v);}
 152   void          SetCovInv(Bool_t v=kTRUE)                         {SetBit(kCovInvBit,v);}
 153   void          SetIgnoreCov(Bool_t v=kTRUE)                      {SetBit(kIgnoreCovBit,v);}
 154   void          SetParAxis(Int_t ax);
 155   void          SetMaxIterations(Int_t n=20)                      {fMaxIter = n<2 ? 2:n;}
 156   void          SetEps(Double_t eps=1e-6)                         {fEps = eps<0 ? GetMachinePrec() : eps;}
 157   void          SetMass(Double_t m=0.13957)                       {fMass = m<5E-4 ? 5E-4 : m;}
 158   void          Reset();
 159   void          BuildMaterialLUT(Int_t ntri=3000);
 160   void          SetCovIScale(Int_t ip, Double_t scl=-1.0);
 161   void          ResetCovIScale(Double_t scl=-1.0)                 {for (int i=fNPBooked;i--;) SetCovIScale(i,scl);}
 162   //
 163   virtual void  Print(Option_t *opt="")                    const;
 164   //
 165   static void   GetNormal(Double_t *norm,const Float_t *covMat);
 166   //
 167  protected:
 168   void          InitAux();
 169   Int_t         ChoseParAxis()                                              const;
 170   Double_t      GetParPCALine(const Double_t *xyz, const Double_t *covI=0/*, Double_t sclCovI=-1*/)  const;
 171   Double_t      GetParPCAHelix(const Double_t *xyz, const Double_t *covI=0, Double_t sclCovI=-1) const;
 172   Double_t      GetParPCACircle(Double_t x, Double_t y)                     const;
 173   Double_t      GetHelixParAtR(Double_t r)                                  const;
 174   //
 175   void          GetDtDPosLine(Double_t *dtpos,/*const Double_t *xyz,*/  const Double_t *covI=0/*, Double_t sclCovI=-1*/)  const;
 176   Double_t      GetDtDParamsLine(Double_t *dtparam,const Double_t *xyz, const Double_t *covI=0/*, Double_t sclCovI=-1*/)  const;
 177   //
 178   Double_t      GetDRofELoss(Double_t t,Double_t cdip,Double_t rhoL,
 179                              const Double_t *normS, Double_t &p,Double_t &e) const;
 180   static Bool_t IsZero(Double_t v,Double_t threshold = 1e-16)     {return TMath::Abs(v)<threshold; }
 181   static Double_t      GetMachinePrec();
 182   //
 183  protected:
 184   const AliTrackPointArray *fkPoints;               // current points
 185   AliParamSolver* fParSol;                         // solver for parametric linearized systems
 186   //
 187   Double_t  fBz;                                   // magnetic field
 188   Int_t     fCharge;                               // track charge +1=+, -1=-
 189   Int_t     fPntFirst;                             // first point to fit
 190   Int_t     fPntLast;                              // last point to fit
 191   Int_t     fNPBooked;                             // number of points booked
 192   Int_t     fParAxis;                              // parameterization axis
 193   Double_t *fCovI;                                 //! inverted cov.matrix for each point
 194   Double_t  fParams[kMaxParam];                    // fitted params
 195   Double_t  fParamsCov[kMaxParamSq];               // fit cov matrix
 196   Double_t  fChi2NDF;                              // fit chi2/NDF
 197   Int_t     fMaxIter;                              // max number of iterations
 198   Int_t     fIter;                                 // real number of iterations
 199   Double_t  fEps;                                  // precision
 200   Double_t  fMass;                                 // assumed particle mass for ELoss Calculation
 201   Bool_t    fSwitch2Line;                          // decided to switch to line
 202   Double_t  fMaxRforHelix;                         // above this radius use straight line fit
 203   //
 204   const Int_t  *fkAxID;                            // axis IDs
 205   const Int_t  *fkAxCID;                           // axis combinations IDs
 206   //
 207   // internal storage
 208   Double_t *fCurT;                                 // track parameter for each point
 209   //
 210   // storage to account e-loss
 211   Int_t     fFirstPosT;                            // id of the first positive t index in fElsId
 212   Int_t     fNElsPnt;                              // number of e-loss layers seen by the track
 213   Int_t    *fElsId;                                // index of increasing t-ordering in the fCurT
 214   Double_t *fElsDR;                                // delta_Radius for each e-loss layer
 215   //
 216   static       Double_t fgRhoLITS[kMaxLrITS];      // <rho*L> for each material layer
 217   static const Double_t fgkRLayITS[kMaxLrITS];     // radii of material layers
 218   static const Double_t fgkZSpanITS[kMaxLrITS];    // half Z span of the material layer
 219   static const Int_t    fgkPassivLrITS[3];         // list of passive layer enums
 220   static const Int_t    fgkActiveLrITS[6];         // list of active layer enums
 221   static const Double_t fgkAlmostZero;             // tiny double
 222   static const Double_t fgkCQConv;                 // R = PT/Bz/fgkCQConv with GeV,kGauss,cm
 223   static const Int_t    fgkAxisID[3][3];           // permutations of axis
 224   static const Int_t    fgkAxisCID[3][6];          // cov matrix elements for axis selection
 225
 226   ClassDef(AliITSTPArrayFit,0);
 227 };
 228
 229 //____________________________________________________
 230 inline void AliITSTPArrayFit::GetPosition(Double_t *xyz, Int_t pnt) const
 231 {
 232   // track position at measured point pnt
 233   GetPosition(xyz,fCurT[pnt]);
 234 }
 235
 236 //____________________________________________________
 237 inline Double_t AliITSTPArrayFit::GetParPCA(const double *xyz, const double *covI, Double_t sclCovI) const
 238 {
 239   // get parameter for the point with least weighted distance to the point
 240   if (IsFieldON()) return GetParPCAHelix(xyz,covI,sclCovI);
 241   else             return GetParPCALine(xyz,covI/*,sclCovI*/);
 242 }
 243
 244 //____________________________________________________
 245 inline Double_t* AliITSTPArrayFit::GetPoint(Int_t ip) const
 246 {
 247   // get point xyz
 248   static double xyz[3];
 249   xyz[kX] = fkPoints->GetX()[ip];
 250   xyz[kY] = fkPoints->GetY()[ip];
 251   xyz[kZ] = fkPoints->GetZ()[ip];
 252   return &xyz[0];
 253 }
 254
 255 //____________________________________________________
 256 inline Double_t AliITSTPArrayFit::Fit(Int_t extQ,Double_t extPT,Double_t extPTerr)
 257 {
 258   // perform the fit
 259   if (IsFieldON()) return FitHelix(extQ,extPT,extPTerr);
 260   else             return FitLine();
 261 }
 262
 263 //____________________________________________________
 264 inline void  AliITSTPArrayFit::SetCovIScale(Int_t ip, Double_t scl)
 265 {
 266   // rescale inverted error matrix of specific point
 267   if (ip>=fNPBooked) return;
 268   if (TMath::Abs(scl-GetCovIScale(ip))<1e-7) ResetBit(kFitDoneBit);
 269   fCovI[ip*kNCov+kScl] = scl;
 270 }
 271
 272 #endif
 273