Line data Source code
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 0 : AliTrackPointArray* GetPoints() const {return (AliTrackPointArray*)fkPoints;}
65 : //
66 0 : void SetBz(Double_t bz) {fBz = bz;}
67 0 : Double_t GetBz() const {return fBz;}
68 0 : Bool_t IsHelix() const {return fParAxis<0;}
69 0 : Bool_t IsFieldON() const {return TMath::Abs(fBz)>1e-5;}
70 0 : Bool_t IsTypeCosmics() const {return TestBit(kCosmicsBit);}
71 0 : Bool_t IsTypeCollision() const {return !IsTypeCosmics();}
72 0 : Int_t GetCharge() const {return fCharge;}
73 0 : 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 0 : 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 0 : 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 0 : Int_t GetParAxis() const {return fParAxis;}
110 0 : Int_t GetAxID(Int_t id) const {return fkAxID ? fkAxID[id] : -1;}
111 0 : Int_t GetAxCID(Int_t id) const {return fkAxCID ? fkAxCID[id] : -1;}
112 0 : Int_t GetFirst() const {return fPntFirst;}
113 0 : Int_t GetLast() const {return fPntLast;}
114 : //
115 0 : Int_t GetNParams() const {return IsFieldON() ? 5:4;}
116 : Bool_t InvertPointsCovMat();
117 : //
118 0 : Int_t* GetElsId() const {return fElsId;}
119 0 : Double_t* GetElsDR() const {return fElsDR;}
120 : //
121 0 : Double_t GetCovIScale(Int_t ip) const {return ip<fNPBooked ? fCovI[ip*kNCov+kScl] : -1.0;}
122 0 : Double_t* GetCovI(Int_t ip) const {return fCovI + kNCov*ip;}
123 0 : Double_t* GetCovI() const {return fCovI;}
124 0 : Double_t* GetParams() const {return (Double_t*)&fParams[0];}
125 0 : Double_t GetParam(Int_t ip) const {return fParams[ip];}
126 0 : Double_t* GetTs() const {return (Double_t*)fCurT;}
127 0 : 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 0 : Bool_t IsSwitched2Line() const {return fSwitch2Line;}
132 0 : Bool_t IsELossON() const {return TestBit(kELossBit)&&IsFieldON();}
133 0 : Bool_t IsFitDone() const {return TestBit(kFitDoneBit);}
134 0 : Bool_t IsCovInv() const {return TestBit(kCovInvBit);}
135 0 : Bool_t IsCovIgnored() const {return TestBit(kIgnoreCovBit);}
136 0 : Int_t GetMaxIterations() const {return fMaxIter;}
137 0 : Int_t GetNIterations() const {return fIter;}
138 0 : Double_t GetEps() const {return fEps;}
139 0 : Double_t GetMass() const {return fMass;}
140 : //
141 : Double_t GetPt() const;
142 : Double_t GetP() const;
143 :
144 : //
145 0 : void Switch2Line(Bool_t v=kTRUE) {fSwitch2Line = v;}
146 0 : void SetMaxRforHelix(Double_t r) {fMaxRforHelix = r>0 ? r : 1e9;}
147 0 : void SetCharge(Int_t q=1) {fCharge = q<0 ? -1:1;}
148 0 : void SetELossON(Bool_t v=kTRUE) {SetBit(kELossBit,v);}
149 0 : void SetTypeCosmics(Bool_t v=kTRUE) {SetBit(kCosmicsBit,v);}
150 0 : void SetTypeCollision(Bool_t v=kTRUE) {SetTypeCosmics(!v);}
151 0 : void SetFitDone(Bool_t v=kTRUE) {SetBit(kFitDoneBit,v);}
152 0 : void SetCovInv(Bool_t v=kTRUE) {SetBit(kCovInvBit,v);}
153 0 : void SetIgnoreCov(Bool_t v=kTRUE) {SetBit(kIgnoreCovBit,v);}
154 : void SetParAxis(Int_t ax);
155 0 : void SetMaxIterations(Int_t n=20) {fMaxIter = n<2 ? 2:n;}
156 0 : void SetEps(Double_t eps=1e-6) {fEps = eps<0 ? GetMachinePrec() : eps;}
157 0 : 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 0 : 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 0 : 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 116 : 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 0 : GetPosition(xyz,fCurT[pnt]);
234 0 : }
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 0 : if (IsFieldON()) return GetParPCAHelix(xyz,covI,sclCovI);
241 0 : else return GetParPCALine(xyz,covI/*,sclCovI*/);
242 0 : }
243 :
244 : //____________________________________________________
245 : inline Double_t* AliITSTPArrayFit::GetPoint(Int_t ip) const
246 : {
247 : // get point xyz
248 : static double xyz[3];
249 0 : xyz[kX] = fkPoints->GetX()[ip];
250 0 : xyz[kY] = fkPoints->GetY()[ip];
251 0 : xyz[kZ] = fkPoints->GetZ()[ip];
252 0 : 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 0 : if (IsFieldON()) return FitHelix(extQ,extPT,extPTerr);
260 0 : else return FitLine();
261 0 : }
262 :
263 : //____________________________________________________
264 : inline void AliITSTPArrayFit::SetCovIScale(Int_t ip, Double_t scl)
265 : {
266 : // rescale inverted error matrix of specific point
267 0 : if (ip>=fNPBooked) return;
268 0 : if (TMath::Abs(scl-GetCovIScale(ip))<1e-7) ResetBit(kFitDoneBit);
269 0 : fCovI[ip*kNCov+kScl] = scl;
270 0 : }
271 :
272 : #endif
273 :
|
