]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSTPArrayFit.h
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated survey information for ladder 2 of SDD layer 3 (M. Poghosyan)
[u/mrichter/AliRoot.git] / ITS / AliITSTPArrayFit.h
CommitLineData
6be22b3f 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"
36class AliSymMatrix;
37class AliLog;
38class AliParamSolver;
39
40
41class 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};
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);
24391cd5 64 AliTrackPointArray* GetPoints() const {return (AliTrackPointArray*)fkPoints;}
6be22b3f 65 //
66 void SetBz(Double_t bz) {fBz = bz;}
67 Double_t GetBz() const {return fBz;}
ef24eb3b 68 Bool_t IsHelix() const {return fParAxis<0;}
69 Bool_t IsFieldON() const {return TMath::Abs(fBz)>1e-5;}
6be22b3f 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) const;
76 Double_t GetPosition( Double_t *xyzPCA, const Double_t* xyz, const Double_t* covI=0) const;
ef24eb3b 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;
6be22b3f 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;
66214d86 83 void GetT0Info(Double_t *xyz, Double_t *dir=0) const;
6be22b3f 84 Double_t CalcChi2NDF() const;
85 Double_t GetChi2NDF() const {return fChi2NDF;}
86 Double_t GetParPCA(const double *xyz, const double *covI=0) const;
1d06ac63 87 Double_t CalcParPCA(Int_t ipnt) const;
66214d86 88 Bool_t CalcErrorMatrix();
6be22b3f 89 //
90 void GetDResDParamsLine (Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0) 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);
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) 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);
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 //
1d06ac63 118 Int_t* GetElsId() const {return fElsId;}
6be22b3f 119 Double_t* GetElsDR() const {return fElsDR;}
120 //
121 Double_t* GetCovI(Int_t ip) const {return fCovI + ip*6;}
122 Double_t* GetCovI() const {return fCovI;}
123 Double_t* GetParams() const {return (Double_t*)&fParams[0];}
124 Double_t GetParam(Int_t ip) const {return fParams[ip];}
125 Double_t* GetTs() const {return (Double_t*)fCurT;}
126 Double_t GetT(Int_t ip) const {return fCurT[ip];}
127 Double_t GetLineOffset(Int_t axis) const;
128 Double_t GetLineSlope(Int_t axis) const;
129 //
1d06ac63 130 Bool_t IsSwitched2Line() const {return fSwitch2Line;}
6be22b3f 131 Bool_t IsELossON() const {return TestBit(kELossBit)&&IsFieldON();}
132 Bool_t IsFitDone() const {return TestBit(kFitDoneBit);}
133 Bool_t IsCovInv() const {return TestBit(kCovInvBit);}
134 Bool_t IsCovIgnored() const {return TestBit(kIgnoreCovBit);}
135 Int_t GetMaxIterations() const {return fMaxIter;}
136 Int_t GetNIterations() const {return fIter;}
137 Double_t GetEps() const {return fEps;}
138 Double_t GetMass() const {return fMass;}
ef24eb3b 139 //
140 Double_t GetPt() const;
141 Double_t GetP() const;
142
6be22b3f 143 //
1d06ac63 144 void Switch2Line(Bool_t v=kTRUE) {fSwitch2Line = v;}
145 void SetMaxRforHelix(Double_t r) {fMaxRforHelix = r>0 ? r : 1e9;}
6be22b3f 146 void SetCharge(Int_t q=1) {fCharge = q<0 ? -1:1;}
147 void SetELossON(Bool_t v=kTRUE) {SetBit(kELossBit,v);}
148 void SetTypeCosmics(Bool_t v=kTRUE) {SetBit(kCosmicsBit,v);}
149 void SetTypeCollision(Bool_t v=kTRUE) {SetTypeCosmics(!v);}
150 void SetFitDone(Bool_t v=kTRUE) {SetBit(kFitDoneBit,v);}
151 void SetCovInv(Bool_t v=kTRUE) {SetBit(kCovInvBit,v);}
152 void SetIgnoreCov(Bool_t v=kTRUE) {SetBit(kIgnoreCovBit,v);}
153 void SetParAxis(Int_t ax);
154 void SetMaxIterations(Int_t n=20) {fMaxIter = n<2 ? 2:n;}
155 void SetEps(Double_t eps=1e-6) {fEps = eps<0 ? GetMachinePrec() : eps;}
156 void SetMass(Double_t m=0.13957) {fMass = m<5E-4 ? 5E-4 : m;}
157 void Reset();
158 void BuildMaterialLUT(Int_t ntri=3000);
159 //
160 virtual void Print(Option_t *opt="") const;
161 //
162 static void GetNormal(Double_t *norm,const Float_t *covMat);
163 //
164 protected:
165 void InitAux();
6be22b3f 166 Int_t ChoseParAxis() const;
167 Double_t GetParPCALine(const Double_t *xyz, const Double_t *covI=0) const;
168 Double_t GetParPCAHelix(const Double_t *xyz, const Double_t *covI=0) const;
169 Double_t GetParPCACircle(Double_t x, Double_t y) const;
170 Double_t GetHelixParAtR(Double_t r) const;
171 //
172 void GetDtDPosLine(Double_t *dtpos,/*const Double_t *xyz,*/ const Double_t *covI=0) const;
173 Double_t GetDtDParamsLine(Double_t *dtparam,const Double_t *xyz, const Double_t *covI=0) const;
174 //
175 Double_t GetDRofELoss(Double_t t,Double_t cdip,Double_t rhoL,
176 const Double_t *normS, Double_t &p,Double_t &e) const;
177 static Bool_t IsZero(Double_t v,Double_t threshold = 1e-16) {return TMath::Abs(v)<threshold; }
24391cd5 178 static Double_t GetMachinePrec();
6be22b3f 179 //
180 protected:
24391cd5 181 const AliTrackPointArray *fkPoints; // current points
6be22b3f 182 AliParamSolver* fParSol; // solver for parametric linearized systems
183 //
184 Double_t fBz; // magnetic field
185 Int_t fCharge; // track charge +1=+, -1=-
186 Int_t fPntFirst; // first point to fit
187 Int_t fPntLast; // last point to fit
188 Int_t fNPBooked; // number of points booked
189 Int_t fParAxis; // parameterization axis
190 Double_t *fCovI; //! inverted cov.matrix for each point
191 Double_t fParams[kMaxParam]; // fitted params
192 Double_t fParamsCov[kMaxParamSq]; // fit cov matrix
193 Double_t fChi2NDF; // fit chi2/NDF
194 Int_t fMaxIter; // max number of iterations
195 Int_t fIter; // real number of iterations
196 Double_t fEps; // precision
197 Double_t fMass; // assumed particle mass for ELoss Calculation
1d06ac63 198 Bool_t fSwitch2Line; // decided to switch to line
199 Double_t fMaxRforHelix; // above this radius use straight line fit
6be22b3f 200 //
201 const Int_t *fkAxID; // axis IDs
202 const Int_t *fkAxCID; // axis combinations IDs
203 //
204 // internal storage
205 Double_t *fCurT; // track parameter for each point
206 //
207 // storage to account e-loss
208 Int_t fFirstPosT; // id of the first positive t index in fElsId
209 Int_t fNElsPnt; // number of e-loss layers seen by the track
210 Int_t *fElsId; // index of increasing t-ordering in the fCurT
211 Double_t *fElsDR; // delta_Radius for each e-loss layer
212 //
213 static Double_t fgRhoLITS[kMaxLrITS]; // <rho*L> for each material layer
214 static const Double_t fgkRLayITS[kMaxLrITS]; // radii of material layers
215 static const Double_t fgkZSpanITS[kMaxLrITS]; // half Z span of the material layer
216 static const Int_t fgkPassivLrITS[3]; // list of passive layer enums
217 static const Int_t fgkActiveLrITS[6]; // list of active layer enums
218 static const Double_t fgkAlmostZero; // tiny double
219 static const Double_t fgkCQConv; // R = PT/Bz/fgkCQConv with GeV,kGauss,cm
220 static const Int_t fgkAxisID[3][3]; // permutations of axis
221 static const Int_t fgkAxisCID[3][6]; // cov matrix elements for axis selection
222
223 ClassDef(AliITSTPArrayFit,0);
224};
225
226//____________________________________________________
227inline void AliITSTPArrayFit::GetPosition(Double_t *xyz, Int_t pnt) const
228{
229 // track position at measured point pnt
230 GetPosition(xyz,fCurT[pnt]);
231}
232
233//____________________________________________________
234inline Double_t AliITSTPArrayFit::GetParPCA(const double *xyz, const double *covI) const
235{
236 // get parameter for the point with least weighted distance to the point
237 if (IsFieldON()) return GetParPCAHelix(xyz,covI);
238 else return GetParPCALine(xyz,covI);
239}
240
241//____________________________________________________
242inline Double_t* AliITSTPArrayFit::GetPoint(Int_t ip) const
243{
244 static double xyz[3];
24391cd5 245 xyz[kX] = fkPoints->GetX()[ip];
246 xyz[kY] = fkPoints->GetY()[ip];
247 xyz[kZ] = fkPoints->GetZ()[ip];
6be22b3f 248 return &xyz[0];
249}
250
251//____________________________________________________
252inline Double_t AliITSTPArrayFit::Fit(Int_t extQ,Double_t extPT,Double_t extPTerr)
253{
254 if (IsFieldON()) return FitHelix(extQ,extPT,extPTerr);
255 else return FitLine();
256}
257
258
259#endif
24391cd5 260