]>
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); | |
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 | //____________________________________________________ | |
227 | inline 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 | //____________________________________________________ | |
234 | inline 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 | //____________________________________________________ | |
242 | inline 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 | //____________________________________________________ | |
252 | inline 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 |