]>
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> | |
b80c197e | 35 | #include <AliTrackPointArray.h> |
6be22b3f | 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}; | |
8102b2c9 | 48 | enum {kXX=0,kXY=1,kXZ=2,kYX=kXY,kYY=3,kYZ=4,kZX=kXZ,kZY=kYZ,kZZ=5,kScl=6,kNCov}; |
6be22b3f | 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; | |
8102b2c9 | 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; | |
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;} | |
8102b2c9 | 86 | Double_t GetParPCA(const double *xyz, const double *covI=0, Double_t sclCovI=-1) const; |
1d06ac63 | 87 | Double_t CalcParPCA(Int_t ipnt) const; |
66214d86 | 88 | Bool_t CalcErrorMatrix(); |
6be22b3f | 89 | // |
8102b2c9 | 90 | void GetDResDParamsLine (Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0/*,Double_t sclCovI=-1*/) const; |
6be22b3f | 91 | void GetDResDParamsLine (Double_t *dXYZdP, Int_t ipnt) const; |
8102b2c9 | 92 | void GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0, Double_t sclCovI=-1); |
6be22b3f | 93 | void GetDResDParams(Double_t *dXYZdP, Int_t ipnt); |
94 | // | |
8102b2c9 | 95 | void GetDResDPosLine (Double_t *dXYZdP,/*const Double_t *xyz,*/ const Double_t *covI=0/*,Double_t sclCovI=-1*/) const; |
6be22b3f | 96 | void GetDResDPosLine (Double_t *dXYZdP, Int_t ipnt) const; |
b80c197e | 97 | void GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI=0, Double_t sclCovI=-1) const; |
6be22b3f | 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 | // | |
8102b2c9 | 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;} | |
6be22b3f | 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 | // | |
1d06ac63 | 131 | Bool_t IsSwitched2Line() const {return fSwitch2Line;} |
6be22b3f | 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;} | |
ef24eb3b | 140 | // |
141 | Double_t GetPt() const; | |
142 | Double_t GetP() const; | |
143 | ||
6be22b3f | 144 | // |
1d06ac63 | 145 | void Switch2Line(Bool_t v=kTRUE) {fSwitch2Line = v;} |
146 | void SetMaxRforHelix(Double_t r) {fMaxRforHelix = r>0 ? r : 1e9;} | |
6be22b3f | 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); | |
8102b2c9 | 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);} | |
6be22b3f | 162 | // |
8102b2c9 | 163 | virtual void Print(Option_t *opt="") const; |
6be22b3f | 164 | // |
165 | static void GetNormal(Double_t *norm,const Float_t *covMat); | |
166 | // | |
167 | protected: | |
168 | void InitAux(); | |
6be22b3f | 169 | Int_t ChoseParAxis() const; |
8102b2c9 | 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; | |
6be22b3f | 172 | Double_t GetParPCACircle(Double_t x, Double_t y) const; |
173 | Double_t GetHelixParAtR(Double_t r) const; | |
174 | // | |
8102b2c9 | 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; | |
6be22b3f | 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; } | |
24391cd5 | 181 | static Double_t GetMachinePrec(); |
6be22b3f | 182 | // |
183 | protected: | |
24391cd5 | 184 | const AliTrackPointArray *fkPoints; // current points |
6be22b3f | 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 | |
1d06ac63 | 201 | Bool_t fSwitch2Line; // decided to switch to line |
202 | Double_t fMaxRforHelix; // above this radius use straight line fit | |
6be22b3f | 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 | //____________________________________________________ | |
8102b2c9 | 237 | inline Double_t AliITSTPArrayFit::GetParPCA(const double *xyz, const double *covI, Double_t sclCovI) const |
6be22b3f | 238 | { |
239 | // get parameter for the point with least weighted distance to the point | |
8102b2c9 | 240 | if (IsFieldON()) return GetParPCAHelix(xyz,covI,sclCovI); |
241 | else return GetParPCALine(xyz,covI/*,sclCovI*/); | |
6be22b3f | 242 | } |
243 | ||
244 | //____________________________________________________ | |
245 | inline Double_t* AliITSTPArrayFit::GetPoint(Int_t ip) const | |
246 | { | |
8102b2c9 | 247 | // get point xyz |
6be22b3f | 248 | static double xyz[3]; |
24391cd5 | 249 | xyz[kX] = fkPoints->GetX()[ip]; |
250 | xyz[kY] = fkPoints->GetY()[ip]; | |
251 | xyz[kZ] = fkPoints->GetZ()[ip]; | |
6be22b3f | 252 | return &xyz[0]; |
253 | } | |
254 | ||
255 | //____________________________________________________ | |
256 | inline Double_t AliITSTPArrayFit::Fit(Int_t extQ,Double_t extPT,Double_t extPTerr) | |
257 | { | |
8102b2c9 | 258 | // perform the fit |
6be22b3f | 259 | if (IsFieldON()) return FitHelix(extQ,extPT,extPTerr); |
260 | else return FitLine(); | |
261 | } | |
262 | ||
8102b2c9 | 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 | } | |
6be22b3f | 271 | |
272 | #endif | |
24391cd5 | 273 |