]>
Commit | Line | Data |
---|---|---|
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" | |
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}; | |
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*)fPoints;} | |
65 | // | |
66 | void SetBz(Double_t bz) {fBz = bz;} | |
67 | Double_t GetBz() const {return fBz;} | |
1d06ac63 | 68 | Bool_t IsFieldON() const {return TMath::Abs(fBz)>1e-5 && !IsSwitched2Line();} |
6be22b3f | 69 | Bool_t IsTypeCosmics() const {return TestBit(kCosmicsBit);} |
70 | Bool_t IsTypeCollision() const {return !IsTypeCosmics();} | |
71 | Int_t GetCharge() const {return fCharge;} | |
72 | Int_t GetSignQB() const {return fBz<0 ? -fCharge:fCharge;} | |
73 | void GetResiduals(Double_t *res, Int_t ipnt) const; | |
74 | void GetResiduals(Double_t *resPCA, const Double_t* xyz, const Double_t* covI=0) const; | |
75 | Double_t GetPosition( Double_t *xyzPCA, const Double_t* xyz, const Double_t* covI=0) const; | |
76 | Double_t GetPosition( Double_t *xyzPCA, const AliTrackPoint *pntCovInv) const; | |
77 | void GetResiduals(Double_t *xyzPCA, const AliTrackPoint *pntCovInv) const; | |
78 | void GetPosition(Double_t *xyz, Double_t t) const; | |
79 | void GetPosition(Double_t *xyz, Int_t pnt) const; | |
80 | void GetDirCos(Double_t *dircos, Double_t t) const; | |
81 | Double_t GetPCA2PlaneInfo(Double_t *xyz, Double_t *dir=0, Int_t axis=kY, Double_t axval=0) const; | |
66214d86 | 82 | void GetT0Info(Double_t *xyz, Double_t *dir=0) const; |
6be22b3f | 83 | Double_t CalcChi2NDF() const; |
84 | Double_t GetChi2NDF() const {return fChi2NDF;} | |
85 | Double_t GetParPCA(const double *xyz, const double *covI=0) const; | |
1d06ac63 | 86 | Double_t CalcParPCA(Int_t ipnt) const; |
66214d86 | 87 | Bool_t CalcErrorMatrix(); |
6be22b3f | 88 | // |
89 | void GetDResDParamsLine (Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0) const; | |
90 | void GetDResDParamsLine (Double_t *dXYZdP, Int_t ipnt) const; | |
91 | void GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0); | |
92 | void GetDResDParams(Double_t *dXYZdP, Int_t ipnt); | |
93 | // | |
94 | void GetDResDPosLine (Double_t *dXYZdP,/*const Double_t *xyz,*/ const Double_t *covI=0) const; | |
95 | void GetDResDPosLine (Double_t *dXYZdP, Int_t ipnt) const; | |
96 | void GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0); | |
97 | void GetDResDPos(Double_t *dXYZdP, Int_t ipnt); | |
98 | // | |
99 | Double_t* GetPoint(int ip) const; | |
100 | Bool_t Converged() const {return fIter<fMaxIter;} | |
101 | // | |
102 | Double_t Fit(Int_t extQ=0, Double_t extPT=-1,Double_t extPTerr=0); | |
103 | Double_t FitLine(); | |
104 | Double_t FitHelix(Int_t extQ=0, Double_t extPT=-1,Double_t extPTerr=0); | |
105 | Bool_t FitLineCrude(); | |
106 | Bool_t FitHelixCrude(Int_t imposedQ=0); | |
107 | // | |
108 | Int_t GetParAxis() const {return fParAxis;} | |
109 | Int_t GetAxID(Int_t id) const {return fkAxID ? fkAxID[id] : -1;} | |
110 | Int_t GetAxCID(Int_t id) const {return fkAxCID ? fkAxCID[id] : -1;} | |
111 | Int_t GetFirst() const {return fPntFirst;} | |
112 | Int_t GetLast() const {return fPntLast;} | |
113 | // | |
114 | Int_t GetNParams() const {return IsFieldON() ? 5:4;} | |
115 | Bool_t InvertPointsCovMat(); | |
116 | // | |
1d06ac63 | 117 | Int_t* GetElsId() const {return fElsId;} |
6be22b3f | 118 | Double_t* GetElsDR() const {return fElsDR;} |
119 | // | |
120 | Double_t* GetCovI(Int_t ip) const {return fCovI + ip*6;} | |
121 | Double_t* GetCovI() const {return fCovI;} | |
122 | Double_t* GetParams() const {return (Double_t*)&fParams[0];} | |
123 | Double_t GetParam(Int_t ip) const {return fParams[ip];} | |
124 | Double_t* GetTs() const {return (Double_t*)fCurT;} | |
125 | Double_t GetT(Int_t ip) const {return fCurT[ip];} | |
126 | Double_t GetLineOffset(Int_t axis) const; | |
127 | Double_t GetLineSlope(Int_t axis) const; | |
128 | // | |
1d06ac63 | 129 | Bool_t IsSwitched2Line() const {return fSwitch2Line;} |
6be22b3f | 130 | Bool_t IsELossON() const {return TestBit(kELossBit)&&IsFieldON();} |
131 | Bool_t IsFitDone() const {return TestBit(kFitDoneBit);} | |
132 | Bool_t IsCovInv() const {return TestBit(kCovInvBit);} | |
133 | Bool_t IsCovIgnored() const {return TestBit(kIgnoreCovBit);} | |
134 | Int_t GetMaxIterations() const {return fMaxIter;} | |
135 | Int_t GetNIterations() const {return fIter;} | |
136 | Double_t GetEps() const {return fEps;} | |
137 | Double_t GetMass() const {return fMass;} | |
138 | // | |
1d06ac63 | 139 | void Switch2Line(Bool_t v=kTRUE) {fSwitch2Line = v;} |
140 | void SetMaxRforHelix(Double_t r) {fMaxRforHelix = r>0 ? r : 1e9;} | |
6be22b3f | 141 | void SetCharge(Int_t q=1) {fCharge = q<0 ? -1:1;} |
142 | void SetELossON(Bool_t v=kTRUE) {SetBit(kELossBit,v);} | |
143 | void SetTypeCosmics(Bool_t v=kTRUE) {SetBit(kCosmicsBit,v);} | |
144 | void SetTypeCollision(Bool_t v=kTRUE) {SetTypeCosmics(!v);} | |
145 | void SetFitDone(Bool_t v=kTRUE) {SetBit(kFitDoneBit,v);} | |
146 | void SetCovInv(Bool_t v=kTRUE) {SetBit(kCovInvBit,v);} | |
147 | void SetIgnoreCov(Bool_t v=kTRUE) {SetBit(kIgnoreCovBit,v);} | |
148 | void SetParAxis(Int_t ax); | |
149 | void SetMaxIterations(Int_t n=20) {fMaxIter = n<2 ? 2:n;} | |
150 | void SetEps(Double_t eps=1e-6) {fEps = eps<0 ? GetMachinePrec() : eps;} | |
151 | void SetMass(Double_t m=0.13957) {fMass = m<5E-4 ? 5E-4 : m;} | |
152 | void Reset(); | |
153 | void BuildMaterialLUT(Int_t ntri=3000); | |
154 | // | |
155 | virtual void Print(Option_t *opt="") const; | |
156 | // | |
157 | static void GetNormal(Double_t *norm,const Float_t *covMat); | |
158 | // | |
159 | protected: | |
160 | void InitAux(); | |
161 | Double_t GetMachinePrec(); | |
162 | Int_t ChoseParAxis() const; | |
163 | Double_t GetParPCALine(const Double_t *xyz, const Double_t *covI=0) const; | |
164 | Double_t GetParPCAHelix(const Double_t *xyz, const Double_t *covI=0) const; | |
165 | Double_t GetParPCACircle(Double_t x, Double_t y) const; | |
166 | Double_t GetHelixParAtR(Double_t r) const; | |
167 | // | |
168 | void GetDtDPosLine(Double_t *dtpos,/*const Double_t *xyz,*/ const Double_t *covI=0) const; | |
169 | Double_t GetDtDParamsLine(Double_t *dtparam,const Double_t *xyz, const Double_t *covI=0) const; | |
170 | // | |
171 | Double_t GetDRofELoss(Double_t t,Double_t cdip,Double_t rhoL, | |
172 | const Double_t *normS, Double_t &p,Double_t &e) const; | |
173 | static Bool_t IsZero(Double_t v,Double_t threshold = 1e-16) {return TMath::Abs(v)<threshold; } | |
174 | ||
175 | // | |
176 | protected: | |
177 | const AliTrackPointArray *fPoints; // current points | |
178 | AliParamSolver* fParSol; // solver for parametric linearized systems | |
179 | // | |
180 | Double_t fBz; // magnetic field | |
181 | Int_t fCharge; // track charge +1=+, -1=- | |
182 | Int_t fPntFirst; // first point to fit | |
183 | Int_t fPntLast; // last point to fit | |
184 | Int_t fNPBooked; // number of points booked | |
185 | Int_t fParAxis; // parameterization axis | |
186 | Double_t *fCovI; //! inverted cov.matrix for each point | |
187 | Double_t fParams[kMaxParam]; // fitted params | |
188 | Double_t fParamsCov[kMaxParamSq]; // fit cov matrix | |
189 | Double_t fChi2NDF; // fit chi2/NDF | |
190 | Int_t fMaxIter; // max number of iterations | |
191 | Int_t fIter; // real number of iterations | |
192 | Double_t fEps; // precision | |
193 | Double_t fMass; // assumed particle mass for ELoss Calculation | |
1d06ac63 | 194 | Bool_t fSwitch2Line; // decided to switch to line |
195 | Double_t fMaxRforHelix; // above this radius use straight line fit | |
6be22b3f | 196 | // |
197 | const Int_t *fkAxID; // axis IDs | |
198 | const Int_t *fkAxCID; // axis combinations IDs | |
199 | // | |
200 | // internal storage | |
201 | Double_t *fCurT; // track parameter for each point | |
202 | // | |
203 | // storage to account e-loss | |
204 | Int_t fFirstPosT; // id of the first positive t index in fElsId | |
205 | Int_t fNElsPnt; // number of e-loss layers seen by the track | |
206 | Int_t *fElsId; // index of increasing t-ordering in the fCurT | |
207 | Double_t *fElsDR; // delta_Radius for each e-loss layer | |
208 | // | |
209 | static Double_t fgRhoLITS[kMaxLrITS]; // <rho*L> for each material layer | |
210 | static const Double_t fgkRLayITS[kMaxLrITS]; // radii of material layers | |
211 | static const Double_t fgkZSpanITS[kMaxLrITS]; // half Z span of the material layer | |
212 | static const Int_t fgkPassivLrITS[3]; // list of passive layer enums | |
213 | static const Int_t fgkActiveLrITS[6]; // list of active layer enums | |
214 | static const Double_t fgkAlmostZero; // tiny double | |
215 | static const Double_t fgkCQConv; // R = PT/Bz/fgkCQConv with GeV,kGauss,cm | |
216 | static const Int_t fgkAxisID[3][3]; // permutations of axis | |
217 | static const Int_t fgkAxisCID[3][6]; // cov matrix elements for axis selection | |
218 | ||
219 | ClassDef(AliITSTPArrayFit,0); | |
220 | }; | |
221 | ||
222 | //____________________________________________________ | |
223 | inline void AliITSTPArrayFit::GetPosition(Double_t *xyz, Int_t pnt) const | |
224 | { | |
225 | // track position at measured point pnt | |
226 | GetPosition(xyz,fCurT[pnt]); | |
227 | } | |
228 | ||
229 | //____________________________________________________ | |
230 | inline Double_t AliITSTPArrayFit::GetParPCA(const double *xyz, const double *covI) const | |
231 | { | |
232 | // get parameter for the point with least weighted distance to the point | |
233 | if (IsFieldON()) return GetParPCAHelix(xyz,covI); | |
234 | else return GetParPCALine(xyz,covI); | |
235 | } | |
236 | ||
237 | //____________________________________________________ | |
238 | inline Double_t* AliITSTPArrayFit::GetPoint(Int_t ip) const | |
239 | { | |
240 | static double xyz[3]; | |
241 | xyz[kX] = fPoints->GetX()[ip]; | |
242 | xyz[kY] = fPoints->GetY()[ip]; | |
243 | xyz[kZ] = fPoints->GetZ()[ip]; | |
244 | return &xyz[0]; | |
245 | } | |
246 | ||
247 | //____________________________________________________ | |
248 | inline Double_t AliITSTPArrayFit::Fit(Int_t extQ,Double_t extPT,Double_t extPTerr) | |
249 | { | |
250 | if (IsFieldON()) return FitHelix(extQ,extPT,extPTerr); | |
251 | else return FitLine(); | |
252 | } | |
253 | ||
254 | ||
255 | #endif |