Fix for the case of non-existent calibration files
[u/mrichter/AliRoot.git] / ITS / AliITSTPArrayFit.cxx
1 /**************************************************************************
2  * Copyright(c) 2009-2011, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 /* $Id$ */
16 ///////////////////////////////////////////////////////////////////////////////////////////////
17 //                                                                                           //
18 // The line is defined by equations (1)                                                      //
19 // a0*z+a1*x-a0*a1=0 and                                                                     //
20 // b0*z+b1*y-b0*b1=0                                                                         //
21 // where x,y,z are NOT the lab axes but z is the lab axis along which the track              //
22 // has the largest lever arm and x,y are the remaining 2 axis in                             //
23 // the order of fgkAxisID[z][0], fgkAxisID[z][1]                                             //
24 // The parameters are fParams[kA0,kB0,kA1,kB1] and the axis chosen as the independent        //
25 // var. is fParAxis (i.e. if fParAxis==kZ, then a0=ax,b0=bx, a1=ay,b1=by)                    //
26 //                                                                                           //
27 //                                                                                           //
28 // The helix is defined by the equations (2)                                                 //
29 // X(t) = (dr+R)*cos(phi0) - (R+sum{dRi})*cos(t+phi0) + sum{dRi*cos(phi0+ti)}                //
30 // Y(t) = (dr+R)*sin(phi0) - (R+sum{dRi})*sin(t+phi0) + sum{dRi*sin(phi0+ti)}                //
31 // Z(t) = dz - (R+sum{dRi})*t*tg(dip) + sum{dRi*ti}*tg(dip)                                  //
32 // where dRi is the change of the radius due to the ELoss at parameter ti                    //
33 //                                                                                           //
34 // Author: ruben.shahoyan@cern.ch                                                            //
35 //                                                                                           //
36 ///////////////////////////////////////////////////////////////////////////////////////////////
37
38 #include "AliITSTPArrayFit.h"
39 #include "AliExternalTrackParam.h"
40 #include "AliSymMatrix.h"
41 #include "AliLog.h"
42 #include "AliParamSolver.h"
43 #include "AliGeomManager.h"
44 #include "AliITSgeomTGeo.h"
45 #include "AliTracker.h"
46 #include <TRandom.h>
47 #include <TArrayD.h>
48
49 ClassImp(AliITSTPArrayFit)
50
51 const Int_t  AliITSTPArrayFit::fgkAxisID[3][3] = { 
52   {AliITSTPArrayFit::kY,AliITSTPArrayFit::kZ,AliITSTPArrayFit::kX},
53   {AliITSTPArrayFit::kZ,AliITSTPArrayFit::kX,AliITSTPArrayFit::kY},
54   {AliITSTPArrayFit::kX,AliITSTPArrayFit::kY,AliITSTPArrayFit::kZ} };
55
56 const Int_t  AliITSTPArrayFit::fgkAxisCID[3][6] = { 
57   {AliITSTPArrayFit::kYY,AliITSTPArrayFit::kYZ,AliITSTPArrayFit::kXY,
58    AliITSTPArrayFit::kZZ,AliITSTPArrayFit::kXZ,AliITSTPArrayFit::kXX},
59   //
60   {AliITSTPArrayFit::kZZ,AliITSTPArrayFit::kXZ,AliITSTPArrayFit::kYZ,
61    AliITSTPArrayFit::kXX,AliITSTPArrayFit::kYX,AliITSTPArrayFit::kYY},
62   //
63   {AliITSTPArrayFit::kXX,AliITSTPArrayFit::kXY,AliITSTPArrayFit::kXZ,
64    AliITSTPArrayFit::kYY,AliITSTPArrayFit::kYZ,AliITSTPArrayFit::kZZ}
65 };
66 //
67
68 const Double_t AliITSTPArrayFit::fgkAlmostZero = 1E-55;
69 const Double_t AliITSTPArrayFit::fgkCQConv = 0.299792458e-3;// R = PT/Bz/fgkCQConv with GeV,kGauss,cm
70 const Double_t AliITSTPArrayFit::fgkZSpanITS[AliITSTPArrayFit::kMaxLrITS] = {
71   36. ,14.1,14.1,  38., 22.2,29.7, 51.   ,43.1,48.9};
72
73 const Double_t AliITSTPArrayFit::fgkRLayITS[AliITSTPArrayFit::kMaxLrITS] = {
74   2.94, 3.9,7.6, 11.04, 15.0,23.9, 29.44 ,38.0,43.0};
75
76 const Int_t    AliITSTPArrayFit::fgkPassivLrITS[3] = 
77   {AliITSTPArrayFit::kLrBeamPime,AliITSTPArrayFit::kLrShield1,AliITSTPArrayFit::kLrShield2};
78
79 const Int_t    AliITSTPArrayFit::fgkActiveLrITS[6] = 
80   {AliITSTPArrayFit::kLrSPD1,AliITSTPArrayFit::kLrSPD2,
81    AliITSTPArrayFit::kLrSDD1,AliITSTPArrayFit::kLrSDD2,
82    AliITSTPArrayFit::kLrSSD1,AliITSTPArrayFit::kLrSSD2};
83
84 Double_t AliITSTPArrayFit::fgRhoLITS[AliITSTPArrayFit::kMaxLrITS] = {
85   1.48e-01, 2.48e-01,2.57e-01, 1.34e-01, 3.34e-01,3.50e-01, 2.22e-01, 2.38e-01,2.25e-01};
86
87 //____________________________________________________
88 AliITSTPArrayFit::AliITSTPArrayFit() :
89   fkPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1),
90   fPntLast(-1),fNPBooked(0),fParAxis(-1),fCovI(0),fChi2NDF(0),
91   fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fSwitch2Line(kFALSE),fMaxRforHelix(6.e5),
92   fkAxID(0),fkAxCID(0),fCurT(0),
93   fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
94 {
95   // default constructor
96   for (int i=kMaxParam;i--;)   fParams[i] = 0;
97   for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0;
98   SetMass();
99 }
100
101 //____________________________________________________
102 AliITSTPArrayFit::AliITSTPArrayFit(Int_t np) :
103   fkPoints(0),fParSol(0),fBz(0),fCharge(0),fPntFirst(-1),
104   fPntLast(-1),fNPBooked(np),fParAxis(-1),fCovI(0),fChi2NDF(0),
105   fMaxIter(20),fIter(0),fEps(1e-6),fMass(0),fSwitch2Line(kFALSE),fMaxRforHelix(6.e5),
106   fkAxID(0),fkAxCID(0),fCurT(0),fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
107 {
108   // constructor with booking of np points
109   for (int i=kMaxParam;i--;)   fParams[i] = 0;
110   for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0;
111   InitAux();
112   SetEps();
113   SetMass();
114   SetMaxIterations();
115 }
116
117 //____________________________________________________
118 AliITSTPArrayFit::AliITSTPArrayFit(const AliITSTPArrayFit &src) : 
119   TObject(src),fkPoints(src.fkPoints),fParSol(0),fBz(src.fBz),
120   fCharge(src.fCharge),fPntFirst(src.fPntFirst),fPntLast(src.fPntLast),fNPBooked(src.fNPBooked),
121   fParAxis(src.fParAxis),fCovI(0),fChi2NDF(0),fMaxIter(20),fIter(0),fEps(0),fMass(src.fMass),
122   fSwitch2Line(src.fSwitch2Line),fMaxRforHelix(src.fMaxRforHelix),fkAxID(0),fkAxCID(0),fCurT(0),
123   fFirstPosT(0),fNElsPnt(0),fElsId(0),fElsDR(0)
124 {
125   // copy constructor
126   InitAux();
127   memcpy(fCovI,src.fCovI,fNPBooked*kNCov*sizeof(Double_t));
128   for (int i=kMaxParam;i--;)   fParams[i] = src.fParams[i];
129   for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i];
130   memcpy(fCurT,src.fCurT,fNPBooked*sizeof(Double_t));
131   SetEps(src.fEps);
132   SetMaxIterations(src.fMaxIter);
133   //  
134 }
135
136 //____________________________________________________
137 AliITSTPArrayFit &AliITSTPArrayFit::operator =(const AliITSTPArrayFit& src)
138 {
139   // assignment operator
140   if (this==&src) return *this;
141   ((TObject*)this)->operator=(src);
142   fkPoints   = src.fkPoints;
143   if (!fParSol) fParSol = new AliParamSolver(*src.fParSol);
144   else *fParSol = *src.fParSol; 
145   fBz       = src.fBz; 
146   fCharge   = src.fCharge;
147   fNPBooked = src.fNPBooked;
148   fPntFirst = src.fPntFirst;
149   fPntLast  = src.fPntLast;
150   InitAux();
151   memcpy(fCovI,src.fCovI,fNPBooked*kNCov*sizeof(Double_t));
152   for (int i=kMaxParam;i--;)   fParams[i] = src.fParams[i];
153   for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i];
154   SetParAxis(src.fParAxis);
155   fNElsPnt   = src.fNElsPnt;
156   fFirstPosT = src.fFirstPosT;
157   memcpy(fCurT  ,src.fCurT  ,fNPBooked*sizeof(Double_t));
158   memcpy(fElsId ,src.fElsId ,fNPBooked*sizeof(Int_t));
159   memcpy(fElsDR ,src.fElsDR ,fNPBooked*sizeof(Double_t));
160   memcpy(fCurT  ,src.fCurT  ,fNPBooked*sizeof(Double_t));
161   SetEps(src.fEps);
162   SetMaxIterations(src.fMaxIter);
163   //
164   return *this;
165   //
166 }
167
168 //____________________________________________________
169 AliITSTPArrayFit::~AliITSTPArrayFit()
170 {
171   // destructor
172   delete   fParSol;
173   delete[] fCovI;
174   delete[] fCurT;
175   delete[] fElsId;
176   delete[] fElsDR;
177 }
178
179 //____________________________________________________
180 void AliITSTPArrayFit::Reset()
181 {
182   // reset to process new track
183   if (fParSol) fParSol->Clear();
184   fkPoints=0; 
185   fNElsPnt = 0;
186   fFirstPosT = 0;
187   //  fBz = 0;
188   fCharge = 0;
189   fIter = 0;
190   fPntFirst=fPntLast=-1; 
191   SetParAxis(-1);
192   fSwitch2Line = kFALSE;
193   ResetBit(kFitDoneBit|kCovInvBit);
194 }
195
196 //____________________________________________________
197 void AliITSTPArrayFit::AttachPoints(const AliTrackPointArray* points, Int_t pfirst,Int_t plast) 
198 {
199   // create from piece of AliTrackPointArray
200   Reset();
201   fkPoints = points;
202   int np = points->GetNPoints();
203   if (fNPBooked<np) {
204     fNPBooked = np;
205     InitAux();
206   }
207   fPntFirst = pfirst<0 ? 0 : pfirst;
208   fPntLast  = plast<fPntFirst ? np-1 : plast;
209   //
210   for (int i=kMaxParam;i--;)   fParams[i] = 0;
211   for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0;
212   //
213   InvertPointsCovMat();
214   ResetCovIScale();
215   //
216 }
217
218 //____________________________________________________
219 Bool_t AliITSTPArrayFit::SetFirstLast(Int_t pfirst,Int_t plast) 
220 {
221   // set first and last point to fit
222   const AliTrackPointArray* pnts = fkPoints;
223   if (!pnts) {AliError("TrackPointArray is not attached yet"); return kFALSE;}
224   AttachPoints(pnts,pfirst,plast);
225   return kTRUE;
226   //
227 }
228
229 //____________________________________________________
230 Bool_t AliITSTPArrayFit::InvertPointsCovMat()
231 {
232   // invert the cov.matrices of the points
233   for (int i=fPntFirst;i<=fPntLast;i++) {
234     //
235     float *cov = (float*)fkPoints->GetCov() + i*6; // pointer on cov.matrix
236     //
237     Double_t t0 = cov[kYY]*cov[kZZ] - cov[kYZ]*cov[kYZ];
238     Double_t t1 = cov[kXY]*cov[kZZ] - cov[kXZ]*cov[kYZ];
239     Double_t t2 = cov[kXY]*cov[kYZ] - cov[kXZ]*cov[kYY];
240     Double_t det = cov[kXX]*t0 - cov[kXY]*t1 + cov[kXZ]*t2;
241     if (IsZero(det,1e-18)) { // one of errors is 0, fix this
242       double norm[3];
243       TGeoHMatrix hcov;
244       TGeoRotation hrot,hrotI;
245       GetNormal(norm,cov);
246       double phi = TMath::ATan2(norm[1],norm[0]);
247       hrot.SetAngles(-phi*TMath::RadToDeg(),0.,0.);
248       hrotI.SetAngles(phi*TMath::RadToDeg(),0.,0.);
249       //      
250       Double_t hcovel[9];
251       hcovel[0] = cov[kXX];
252       hcovel[1] = cov[kXY];
253       hcovel[2] = cov[kXZ];
254       hcovel[3] = cov[kXY];
255       hcovel[4] = cov[kYY];
256       hcovel[5] = cov[kYZ];
257       hcovel[6] = cov[kXZ];
258       hcovel[7] = cov[kYZ];
259       hcovel[8] = cov[kZZ];
260       hcov.SetRotation(hcovel);
261       //
262       Double_t *hcovscl = hcov.GetRotationMatrix(); 
263       //      printf(">> %+e %+e %+e\n   %+e %+e %+e\n   %+e %+e %+e\n\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[3],hcovscl[4],hcovscl[5],hcovscl[6],hcovscl[7],hcovscl[8]);
264       //      printf("Rot by %+.e (%+.e %+.e %+.e)\n",phi, norm[0],norm[1],norm[2]);
265       hcov.Multiply(&hrotI);
266       hcov.MultiplyLeft(&hrot);
267       //      printf("|| %+e %+e %+e\n   %+e %+e %+e\n   %+e %+e %+e\n\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[3],hcovscl[4],hcovscl[5],hcovscl[6],hcovscl[7],hcovscl[8]);
268       if (hcovscl[0]<1e-8) hcovscl[0] = 1e-8;
269       if (hcovscl[4]<1e-8) hcovscl[4] = 1e-8;
270       if (hcovscl[8]<1e-8) hcovscl[8] = 1e-8;
271       //      printf("** %+e %+e %+e\n   %+e %+e %+e\n   %+e %+e %+e\n\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[3],hcovscl[4],hcovscl[5],hcovscl[6],hcovscl[7],hcovscl[8]);
272       hcov.Multiply(&hrot);
273       hcov.MultiplyLeft(&hrotI);
274       //      printf("^^ %+e %+e %+e\n   %+e %+e %+e\n   %+e %+e %+e\n\n",hcovscl[0],hcovscl[1],hcovscl[2],hcovscl[3],hcovscl[4],hcovscl[5],hcovscl[6],hcovscl[7],hcovscl[8]);
275       cov[kXX] = hcovscl[0];
276       cov[kXY] = hcovscl[1];
277       cov[kXZ] = hcovscl[2];
278       cov[kYY] = hcovscl[4];
279       cov[kYZ] = hcovscl[5];
280       cov[kZZ] = hcovscl[8];
281     }
282     t0 = cov[kYY]*cov[kZZ] - cov[kYZ]*cov[kYZ];
283     t1 = cov[kXY]*cov[kZZ] - cov[kXZ]*cov[kYZ];
284     t2 = cov[kXY]*cov[kYZ] - cov[kXZ]*cov[kYY];
285     det = cov[kXX]*t0 - cov[kXY]*t1 + cov[kXZ]*t2;
286     //
287     AliDebug(2,Form("%+.4e %+.4e %+.4e -> %+.4e",t0,t1,t2,det));
288     if (IsZero(det,fgkAlmostZero)) {
289       AliInfo(Form("Cov.Matrix for point %d is singular",i));
290       return kFALSE;
291     }
292     //
293     Double_t *covI = GetCovI(i);
294     covI[kXX] =  t0/det;
295     covI[kXY] = -t1/det;
296     covI[kXZ] =  t2/det;
297     covI[kYY] =  (cov[kXX]*cov[kZZ] - cov[kXZ]*cov[kXZ])/det;
298     covI[kYZ] =  (cov[kXY]*cov[kXZ] - cov[kXX]*cov[kYZ])/det;
299     covI[kZZ] =  (cov[kXX]*cov[kYY] - cov[kXY]*cov[kXY])/det;
300     //
301   }
302   SetCovInv();
303   return kTRUE;
304 }
305
306 //____________________________________________________
307 void AliITSTPArrayFit::InitAux()
308 {
309   // init auxiliary space
310   if (fCovI) delete[] fCovI;
311   if (fCurT) delete[] fCurT;
312   //
313   fCovI   = new Double_t[kNCov*fNPBooked];
314   fCurT   = new Double_t[fNPBooked+kMaxLrITS];
315   fElsId  = new Int_t[fNPBooked+kMaxLrITS];
316   fElsDR  = new Double_t[fNPBooked+kMaxLrITS];
317   memset(fElsDR,0,(fNPBooked+kMaxLrITS)*sizeof(Double_t));
318   memset(fCovI,0,fNPBooked*kNCov*sizeof(Double_t));
319   ResetCovIScale();
320   //
321 }
322
323 //____________________________________________________
324 Bool_t AliITSTPArrayFit::FitLineCrude()
325 {
326   // perform linear fit w/o accounting the errors
327   // fit is done in the parameterization
328   // x = res[0] + res[1]*z
329   // y = res[2] + res[3]*z
330   // where x,y,z are NOT the lab axes but z is the lab axis along which the track 
331   // has the largest lever arm and x,y are the remaining 2 axis in 
332   // the order of fgkAxisID[z][0], fgkAxisID[z][1]
333   //
334   int np = fPntLast - fPntFirst + 1;
335   if (np<2) {
336     AliError("At least 2 points are needed for straight line fit");
337     return kFALSE;
338   }
339   //
340   if (fParAxis<0) SetParAxis(ChoseParAxis());
341   Double_t sZ=0,sZZ=0,sY=0,sYZ=0,sX=0,sXZ=0,det=0;
342   //
343   const float *coord[3] = {fkPoints->GetX(),fkPoints->GetY(),fkPoints->GetZ()};
344   const Float_t *varZ = coord[ fParAxis  ];
345   const Float_t *varX = coord[ fkAxID[kX] ];
346   const Float_t *varY = coord[ fkAxID[kY] ];
347   //
348   for (int i=fPntFirst;i<=fPntLast;i++) {
349     sZ  += varZ[i];
350     sZZ += varZ[i]*varZ[i];
351     //
352     sX  += varX[i];
353     sXZ += varX[i]*varZ[i];
354     //
355     sY  += varY[i];
356     sYZ += varY[i]*varZ[i];
357   }
358   det = sZZ*np-sZ*sZ;
359   if (TMath::Abs(det)<fgkAlmostZero) return kFALSE;
360   fParams[0] = (sX*sZZ-sZ*sXZ)/det;
361   fParams[1] = (sXZ*np-sZ*sX)/det;
362   //
363   fParams[2] = (sY*sZZ-sZ*sYZ)/det;
364   fParams[3] = (sYZ*np-sZ*sY)/det;
365   //
366   return kTRUE;
367   //
368 }
369
370 //____________________________________________________
371 void AliITSTPArrayFit::SetParAxis(Int_t ax)
372 {
373   // select the axis which will be used as a parameter for the line: longest baseline
374   if (ax>kZ) {
375     AliInfo(Form("Wrong axis choice: %d",ax));
376     fParAxis = -1;
377   }
378   fParAxis = ax;
379   if (ax>=0) {
380     fkAxID  = fgkAxisID[ax];
381     fkAxCID = fgkAxisCID[ax];
382   }
383   else {
384     fkAxID = fkAxCID = 0;
385   }
386   //
387 }
388
389 //____________________________________________________
390 Int_t AliITSTPArrayFit::ChoseParAxis() const
391 {
392   // select the variable with largest base as a parameter
393   Double_t cmn[3]={1.e9,1.e9,1.e9},cmx[3]={-1.e9,-1.e9,-1.e9};
394   //
395   const float *coord[3] = {fkPoints->GetX(),fkPoints->GetY(),fkPoints->GetZ()};
396   for (int i=fPntFirst;i<=fPntLast;i++) {
397     for (int j=3;j--;) {
398       Double_t val = coord[j][i];
399       if (cmn[j]>val) cmn[j] = val;
400       if (cmx[j]<val) cmx[j] = val;
401     }
402   }
403   //
404   int axis = kZ;
405   if (cmx[axis]-cmn[axis] < cmx[kX]-cmn[kX]) axis = kX;
406   if (cmx[axis]-cmn[axis] < cmx[kY]-cmn[kY]) axis = kY;
407   return axis;
408   //
409 }
410
411 //____________________________________________________
412 Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const Double_t *xyz, const Double_t *covI, Double_t sclCovI) const
413 {
414   // calculate the position of the track at PCA to xyz
415   Double_t t = GetParPCA(xyz,covI,sclCovI);
416   GetPosition(xyzPCA,t);
417   return t;
418 }
419
420 //____________________________________________________
421 Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const AliTrackPoint *pntCovInv, Bool_t useErr) const
422 {
423   // calculate the position of the track at PCA to pntCovInv
424   // NOTE: the covariance matrix of the point must be inverted
425   double t;
426   double xyz[3] = {pntCovInv->GetX(),pntCovInv->GetY(),pntCovInv->GetZ()};
427   if (useErr) {
428     Double_t covI[6];;
429     for (int i=6;i--;) covI[i] = pntCovInv->GetCov()[i];
430     t = GetParPCA(xyz,covI);
431   }
432   else t = GetParPCA(xyz);
433   GetPosition(xyzPCA,t);
434   return t;
435 }
436
437 //____________________________________________________
438 void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const AliTrackPoint *pntCovInv, Bool_t useErr) const
439 {
440   // calculate the residuals  of the track at PCA to pntCovInv
441   // NOTE: the covariance matrix of the point must be inverted
442   GetPosition(resPCA,pntCovInv,useErr);
443   resPCA[0] -= pntCovInv->GetX();
444   resPCA[1] -= pntCovInv->GetY();
445   resPCA[2] -= pntCovInv->GetZ();
446 }
447
448 //____________________________________________________
449 void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const Double_t *xyz, const Double_t *covI, Double_t sclCovI) const
450 {
451   // calculate the residuals of the track at PCA to xyz
452   GetPosition(resPCA,xyz,covI,sclCovI);
453   resPCA[kX] -= xyz[kX];
454   resPCA[kY] -= xyz[kY];
455   resPCA[kZ] -= xyz[kZ];
456 }
457
458 //____________________________________________________
459 Double_t AliITSTPArrayFit::GetParPCALine(const Double_t *xyz, const Double_t *covI/*, Double_t sclCovI*/) const
460 {
461   // get parameter for the point with least weighted distance to the point
462   //
463   Double_t rhs,denom;
464   Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ];
465   Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ];
466   Double_t dz =             -xyz[ fkAxID[kZ] ];
467   //
468   if (covI) {
469     Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ];
470     Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ];
471     Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ];
472     rhs   = tx*dx + ty*dy + tz*dz;
473     denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]);
474   }
475   else {
476     rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz;
477     denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1);
478   }
479   //
480   return rhs/denom;
481   //
482 }
483
484 //____________________________________________________
485 void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, /*const Double_t *xyz,*/ const Double_t *covI/*,Double_t sclCovI*/) const
486 {
487   // calculate detivative of the PCA residuals vs point position and fill in user provide
488   // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp}
489   //
490   Double_t dTdP[3];
491   GetDtDPosLine(dTdP, /*xyz,*/ covI/*,sclCovI*/); // derivative of the t-param over point position
492   //
493   for (int i=3;i--;) {
494     int var = fkAxID[i];
495     Double_t *curd = dXYZdP + var*3;   // d/dCoord_i
496     curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[var];
497     curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[var];
498     curd[ fkAxID[kZ] ] = dTdP[var];
499     curd[    var     ]-= 1.;
500   }
501   //
502 }
503
504 //____________________________________________________
505 void AliITSTPArrayFit::GetDResDParamsLine(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI/*,Double_t sclCovI*/) const
506 {
507   // calculate detivative of the PCA residuals vs line parameters and fill in user provide
508   // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn}
509   //
510   Double_t dTdP[4];
511   Double_t t = GetDtDParamsLine(dTdP, xyz, covI /*,sclCovI*/); // derivative of the t-param over line params
512   //
513   Double_t *curd = dXYZdP + kA0*3; // d/dA0
514   curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kA0] + 1.;
515   curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kA0];
516   curd[ fkAxID[kZ] ] = dTdP[kA0];
517   //
518   curd = dXYZdP + kB0*3; // d/dB0
519   curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kB0] + t;
520   curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kB0];
521   curd[ fkAxID[kZ] ] = dTdP[kB0];
522   //
523   curd = dXYZdP + kA1*3; // d/dA1
524   curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kA1];
525   curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kA1] + 1.;
526   curd[ fkAxID[kZ] ] = dTdP[kA1];
527   //
528   curd = dXYZdP + kB1*3; // d/dB1
529   curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kB1];
530   curd[ fkAxID[kY] ] = fParams[kB1]*dTdP[kB1] + t;
531   curd[ fkAxID[kZ] ] = dTdP[kB1];
532   //
533 }
534
535 //____________________________________________________
536 Double_t AliITSTPArrayFit::GetDtDParamsLine(Double_t *dtparam,const Double_t *xyz, const Double_t *covI) const
537 {
538   // get t-param detivative over the parameters for the point with least weighted distance to the point
539   //
540   Double_t rhs,denom;
541   Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ];
542   Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ];
543   Double_t dz =             -xyz[ fkAxID[kZ] ];
544   Double_t rhsDA0,rhsDA1,rhsDB0,rhsDB1,denDB0,denDB1;
545   //
546   if (covI) {
547     Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ];
548     Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ];
549     Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ];
550     rhs = tx*dx + ty*dy + tz*dz;
551     denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]);
552     //
553     rhsDA0 = tx;
554     rhsDA1 = ty;
555     rhsDB0 = covI[ fkAxCID[kXX] ]*dx + covI[ fkAxCID[kXY] ]*dy + covI[ fkAxCID[kXZ] ]*dz;
556     rhsDB1 = covI[ fkAxCID[kXY] ]*dx + covI[ fkAxCID[kYY] ]*dy + covI[ fkAxCID[kYZ] ]*dz;
557     //
558     denDB0 = -(tx + tx);
559     denDB1 = -(ty + ty);
560   }
561   else {
562     rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz;
563     denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1);
564     //
565     rhsDA0 = fParams[kB0];
566     rhsDB0 = dx;
567     rhsDA1 = fParams[kB1];
568     rhsDB1 = dy;
569     //
570     denDB0 = -(fParams[kB0]+fParams[kB0]);
571     denDB1 = -(fParams[kB1]+fParams[kB1]);
572     //
573   }
574   //
575   Double_t denom2 = denom*denom;
576   dtparam[kA0] = rhsDA0/denom;    // denom does not depend on A0,A1
577   dtparam[kA1] = rhsDA1/denom;
578   dtparam[kB0] = rhsDB0/denom - rhs/denom2 * denDB0;
579   dtparam[kB1] = rhsDB1/denom - rhs/denom2 * denDB1;
580   //
581   return rhs/denom;
582 }
583
584 //____________________________________________________
585 void AliITSTPArrayFit::GetDtDPosLine(Double_t *dtpos,/*const Double_t *xyz,*/ const Double_t *covI) const
586 {
587   // get t-param detivative over the parameters for the point with least weighted distance to the point
588   //
589   //  Double_t rhs;
590   //  Double_t dx = fParams[kA0]-xyz[ fkAxID[kX] ];
591   //  Double_t dy = fParams[kA1]-xyz[ fkAxID[kY] ];
592   //  Double_t dz =             -xyz[ fkAxID[kZ] ];
593   Double_t denom;
594   Double_t rhsDX,rhsDY,rhsDZ;
595   //
596   if (covI) {
597     Double_t tx = fParams[kB0]*covI[ fkAxCID[kXX] ] + fParams[kB1]*covI[ fkAxCID[kXY] ] + covI[ fkAxCID[kXZ] ];
598     Double_t ty = fParams[kB0]*covI[ fkAxCID[kXY] ] + fParams[kB1]*covI[ fkAxCID[kYY] ] + covI[ fkAxCID[kYZ] ];
599     Double_t tz = fParams[kB0]*covI[ fkAxCID[kXZ] ] + fParams[kB1]*covI[ fkAxCID[kYZ] ] + covI[ fkAxCID[kZZ] ];
600     // rhs = tx*dx + ty*dy + tz*dz;
601     denom = -(fParams[kB0]*(covI[ fkAxCID[kXZ] ] + tx) + fParams[kB1]*(covI[ fkAxCID[kYZ] ] + ty) + covI[ fkAxCID[kZZ] ]);
602     //
603     rhsDX = -tx;
604     rhsDY = -ty;
605     rhsDZ = -tz;
606   }
607   else {
608     // rhs = fParams[kB0]*dx + fParams[kB1]*dy + dz;
609     denom = -(fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1);
610     //
611     rhsDX = -fParams[kB0];
612     rhsDY = -fParams[kB1];
613     rhsDZ = -1;
614     //
615   }
616   //
617   dtpos[ fkAxID[kX] ] = rhsDX/denom;
618   dtpos[ fkAxID[kY] ] = rhsDY/denom;
619   dtpos[ fkAxID[kZ] ] = rhsDZ/denom;
620   //
621   //  return rhs/denom;
622 }
623
624 //____________________________________________________
625 void AliITSTPArrayFit::GetDResDParamsLine(Double_t *dXYZdP, Int_t ipnt) const
626 {
627   // calculate detivative of the PCA residuals vs line parameters and fill in user provide
628   // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn}
629   //
630   if (ipnt<fPntFirst || ipnt>fPntLast) {
631     AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
632     return;
633   }
634   GetDResDParamsLine(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt));
635 }
636
637 //____________________________________________________
638 void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, Int_t ipnt) const
639 {
640   // calculate detivative of the PCA residuals vs point position and fill in user provide
641   // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp}
642   //
643   if (ipnt<fPntFirst || ipnt>fPntLast) {
644     AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
645     return;
646   }
647   GetDResDPosLine(dXYZdP,IsCovIgnored() ? 0 : GetCovI(ipnt)/*,GetCovIScale(ipnt)*/);
648 }
649
650 //____________________________________________________
651 void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, Int_t ipnt)
652 {
653   // calculate detivative of the PCA residuals vs track parameters and fill in user provide
654   // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn}
655   //
656   if (ipnt<fPntFirst || ipnt>fPntLast) {
657     AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
658     return;
659   }
660   GetDResDParams(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt),GetCovIScale(ipnt));
661 }
662
663 //____________________________________________________
664 void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, Int_t ipnt)
665 {
666   // calculate detivative of the PCA residuals vs point position and fill in user provide
667   // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp} 
668   //
669   if (ipnt<fPntFirst || ipnt>fPntLast) {
670     AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
671     return;
672   }
673   GetDResDPos(dXYZdP, GetPoint(ipnt), IsCovIgnored() ? 0 : GetCovI(ipnt), GetCovIScale(ipnt));
674 }
675
676 //____________________________________________________
677 void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI, Double_t sclCovI) 
678 {
679   // get residual detivatives over the track parameters for the point with least weighted distance to the point
680   //
681   if (!IsHelix()) { // for the straight line calculate analytically
682     GetDResDParamsLine(dXYZdP, xyz, covI /*,sclCovI*/);
683     return;
684   }
685   //
686   // calculate derivative numerically
687   const Double_t delta = 0.01;
688   Double_t xyzVar[4][3];
689   //
690   for (int ipar = 5;ipar--;) {
691     double sav = fParams[ipar];
692     fParams[ipar] -= delta;
693     GetPosition(xyzVar[0],xyz,covI,sclCovI);
694     fParams[ipar] += delta/2;
695     GetPosition(xyzVar[1],xyz,covI,sclCovI);
696     fParams[ipar] += delta;
697     GetPosition(xyzVar[2],xyz,covI,sclCovI);
698     fParams[ipar] += delta/2;
699     GetPosition(xyzVar[3],xyz,covI,sclCovI);
700     fParams[ipar] = sav;  // restore
701     //
702     double *curd = dXYZdP + 3*ipar;
703     for (int i=3;i--;) curd[i] = (8.*(xyzVar[2][i]-xyzVar[1][i]) - (xyzVar[3][i]-xyzVar[0][i]))/6./delta;
704   }
705   //
706 }
707
708
709 //____________________________________________________
710 void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI,Double_t sclCovI) const
711 {
712   // get residuals detivative over the point position for the point with least weighted distance to the point
713   //
714
715   if (!IsHelix()) { // for the straight line calculate analytically
716     GetDResDPosLine(dXYZdP, /*xyz,*/ covI /*,sclCovI*/);
717     return;
718   }
719   //
720   // calculate derivative numerically
721   const Double_t delta = 0.005;
722   Double_t xyzVar[4][3];
723   Double_t xyzv[3] = {xyz[0],xyz[1],xyz[2]};
724   //
725   for (int ipar = 3;ipar--;) {
726     double sav = xyzv[ipar];
727     xyzv[ipar] -= delta;
728     GetPosition(xyzVar[0],xyzv,covI,sclCovI);
729     xyzv[ipar] += delta/2;
730     GetPosition(xyzVar[1],xyzv,covI,sclCovI);
731     xyzv[ipar] += delta;
732     GetPosition(xyzVar[2],xyzv,covI,sclCovI);
733     xyzv[ipar] += delta/2;
734     GetPosition(xyzVar[3],xyzv,covI,sclCovI);
735     xyzv[ipar] = sav;  // restore
736     //
737     double *curd = dXYZdP + 3*ipar;
738     for (int i=3;i--;) curd[i] = (8.*(xyzVar[2][i]-xyzVar[1][i]) - (xyzVar[3][i]-xyzVar[0][i]))/6./delta;
739     curd[ipar] -= 1.;
740   }
741   //
742 }
743
744 //________________________________________________________________________________________________________
745 Double_t AliITSTPArrayFit::GetParPCAHelix(const Double_t* xyz, const Double_t* covI,Double_t sclCovI) const
746 {
747   // find track parameter t (eq.2) corresponding to point of closest approach to xyz
748   //
749   Double_t phi  = GetParPCACircle(xyz[kX],xyz[kY]); 
750   Double_t cs = TMath::Cos(fParams[kPhi0]);
751   Double_t sn = TMath::Sin(fParams[kPhi0]);
752   Double_t xc = (fParams[kD0]+fParams[kR0])*cs;
753   Double_t yc = (fParams[kD0]+fParams[kR0])*sn;
754   Double_t dchi2,ddchi2;
755   //
756   Double_t dzD  = -fParams[kR0]*fParams[kDip];
757   Double_t dphi = 0;
758   //
759   double rEps = 1e-5/TMath::Abs(fParams[kR0]); // dphi corresponding to 0.1 micron
760   if (rEps>fEps) rEps = fEps;
761   //
762   int it=0;
763   do {
764     cs = TMath::Cos(phi + fParams[kPhi0]);
765     sn = TMath::Sin(phi + fParams[kPhi0]);
766     //
767     Double_t dxD  =  fParams[kR0]*sn;
768     Double_t dyD  = -fParams[kR0]*cs;
769     Double_t dxDD = -dyD;
770     Double_t dyDD =  dxD;
771     //
772     Double_t dx   = xc - fParams[kR0]*cs - xyz[kX];
773     Double_t dy   = yc - fParams[kR0]*sn - xyz[kY];
774     Double_t dz   = fParams[kDZ] + dzD*phi- xyz[kZ];
775     //
776     if (covI) {
777       Double_t tx = dx*covI[kXX] + dy*covI[kXY] + dz*covI[kXZ];
778       Double_t ty = dx*covI[kXY] + dy*covI[kYY] + dz*covI[kYZ];
779       Double_t tz = dx*covI[kXZ] + dy*covI[kYZ] + dz*covI[kZZ];
780       //
781       Double_t ttx = dxD*covI[kXX] + dyD*covI[kXY] + dzD*covI[kXZ];
782       Double_t tty = dxD*covI[kXY] + dyD*covI[kYY] + dzD*covI[kYZ];
783       Double_t ttz = dxD*covI[kXZ] + dyD*covI[kYZ] + dzD*covI[kZZ];
784       //
785       // chi2   = dx*tx + dy*ty + dz*tz;
786       dchi2  = dxD*tx  + dyD*ty  + dzD*tz;
787       ddchi2 = dxDD*tx + dyDD*ty           + dxD *ttx + dyD *tty + dzD *ttz;
788       //
789       if (sclCovI>0) {dchi2 *= sclCovI; ddchi2 *= sclCovI;}
790     }
791     else {
792       // chi2   = dx*dx + dy*dy + dz*dz;
793       dchi2  = dxD*dx  + dyD*dy  + dzD*dz;
794       ddchi2 = dxDD*dx + dyDD*dy +         + dxD*dxD + dyD*dyD + dzD*dzD;
795     }
796     //
797     if (TMath::Abs(ddchi2)<fgkAlmostZero || TMath::Abs(dphi=dchi2/ddchi2)<rEps) break;
798     phi -= dphi;    
799   } while(++it<fMaxIter);
800
801   //
802   return phi;
803 }
804
805 //________________________________________________________________________________________________________
806 Double_t AliITSTPArrayFit::GetParPCACircle(Double_t x,Double_t y)  const
807 {
808   // find track parameter t (eq.2) corresponding to point on the circle with closest approach to x,y
809   //
810   Double_t r = fParams[kD0]+fParams[kR0];
811   Double_t t = TMath::ATan2( r*TMath::Sin(fParams[kPhi0])-y, r*TMath::Cos(fParams[kPhi0])-x ) - fParams[kPhi0]; 
812   if (fParams[kR0] < 0) t += TMath::Pi();
813   if (t > TMath::Pi())  t -= TMath::Pi()*2;
814   if (t <-TMath::Pi())  t += TMath::Pi()*2;
815   return t;
816 }
817
818 //________________________________________________________________________________________________________
819 Double_t AliITSTPArrayFit::GetHelixParAtR(Double_t r)  const
820 {
821   // find helix parameter t (eq.2) corresponding to point on the circle of radius t
822   //
823   double gam = 1. - (r-fParams[kD0])*(r+fParams[kD0])/fParams[kR0]/(fParams[kD0]+fParams[kR0])/2.;
824   return (TMath::Abs(gam)>1) ?  -1e9 : TMath::ACos(gam);
825 }
826
827 //________________________________________________________________________________________________________
828 Double_t AliITSTPArrayFit::CalcChi2NDF() const
829 {
830   // calculate fit chi2/ndf
831   Double_t chi2 = 0;
832   Double_t dr[3]; // residuals
833   //if (!IsFitDone()) return -1;
834   for (int ipnt=fPntFirst;ipnt<=fPntLast;ipnt++) {
835     GetResiduals(dr,ipnt);
836     Double_t* covI = GetCovI(ipnt);
837     Double_t chi2p = dr[kX]*(dr[kX]*covI[ kXX ]+dr[kY]*covI[ kXY ]+dr[kZ]*covI[ kXZ ])
838       +              dr[kY]*(dr[kX]*covI[ kXY ]+dr[kY]*covI[ kYY ]+dr[kZ]*covI[ kYZ ])
839       +              dr[kZ]*(dr[kX]*covI[ kXZ ]+dr[kY]*covI[ kYZ ]+dr[kZ]*covI[ kZZ ]);
840     if (covI[kScl]>0) chi2p *= covI[kScl]; // rescaling was requested for this point's errors
841     chi2 += chi2p;
842   }
843   int ndf = (fPntLast-fPntFirst+1)*3 - GetNParams();
844   chi2 /= ndf;
845   return chi2;
846 }
847
848 //________________________________________________________________________________________________________
849 void AliITSTPArrayFit::GetResiduals(Double_t *res,Int_t ipnt) const
850 {
851   // calculate residuals at point
852   if (ipnt<fPntFirst || ipnt>fPntLast) {
853     AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast));
854     return;
855   }
856   GetPosition(res,fCurT[ipnt]);
857   res[kX] -= fkPoints->GetX()[ipnt];
858   res[kY] -= fkPoints->GetY()[ipnt];
859   res[kZ] -= fkPoints->GetZ()[ipnt];
860 }
861
862 //________________________________________________________________________________________________________
863 void AliITSTPArrayFit::GetPosition(Double_t *xyz, Double_t t) const
864 {
865   // calculate track position for parameter value t
866   if (IsHelix()) {
867     //
868     Double_t rrho = fParams[kD0]+fParams[kR0];
869     Double_t xc = rrho*TMath::Cos(fParams[kPhi0]);
870     Double_t yc = rrho*TMath::Sin(fParams[kPhi0]);
871     Double_t r  = fParams[kR0];
872     Double_t ze = 0;
873     //
874     if (IsELossON()) {
875       if (t>0) {
876         for (int i=fFirstPosT;i<fNElsPnt;i++) { // along the track direction
877           int indE = fElsId[i];
878           if ( t<fCurT[indE] ) break;       // does not reach this layer on  its way to t 
879           xc += fElsDR[indE] * TMath::Cos(fParams[kPhi0] + fCurT[indE]);
880           yc += fElsDR[indE] * TMath::Sin(fParams[kPhi0] + fCurT[indE]);
881           ze += fElsDR[indE] * fCurT[indE];
882           r  += fElsDR[indE];
883           //printf("ELoss@ %+.2e r:%+.3e got %+.3e\n",fCurT[indE],r,fElsDR[indE]);
884         }
885       } else {
886         for (int i=fFirstPosT;i--;) { // against the track direction
887           int indE = fElsId[i];
888           if ( t>=fCurT[indE] ) break;       // does not reach this layer on  its way to t 
889           xc += fElsDR[indE] * TMath::Cos(fParams[kPhi0] + fCurT[indE]);
890           yc += fElsDR[indE] * TMath::Sin(fParams[kPhi0] + fCurT[indE]);
891           ze += fElsDR[indE] * fCurT[indE];
892           r  += fElsDR[indE];
893           //printf("ELoss@ %+.2e r:%+.3e got %+.3e\n",fCurT[indE],r,fElsDR[indE]);
894         }     
895       }
896     }
897     //
898     xyz[kZ] = fParams[kDZ] - fParams[kDip]*(t*r - ze);
899     //
900     t += fParams[kPhi0];    
901     xyz[kX] = xc - r*TMath::Cos(t);
902     xyz[kY] = yc - r*TMath::Sin(t);
903     //    printf("t: %+.3e xyz:%+.2e %+.2e %+.2e | R %+.6e -> %+.6e | sign %d\n",t-fParams[kPhi0],xyz[0],xyz[1],xyz[2],fParams[kR0],r,GetSignQB());
904   }
905   else {
906     xyz[ fkAxID[kX] ] = fParams[kA0] + fParams[kB0]*t;
907     xyz[ fkAxID[kY] ] = fParams[kA1] + fParams[kB1]*t;
908     xyz[ fParAxis   ] = t;
909   }
910 }
911
912 //________________________________________________________________________________________________________
913 void AliITSTPArrayFit::GetDirCos(Double_t *dircos, Double_t t) const
914 {
915   // calculate track direction cosines for parameter value t
916   if (IsHelix()) {
917     dircos[kZ] = -fParams[kDip];
918     t += fParams[kPhi0];    
919     dircos[kX] = TMath::Sin(t);
920     dircos[kY] =-TMath::Cos(t);
921     double gam = TMath::Sign(1/TMath::Sqrt(dircos[kZ]*dircos[kZ]+dircos[kY]*dircos[kY]+dircos[kX]*dircos[kX]),fParams[kR0]);
922     for (int i=3;i--;) dircos[i] *= gam;
923     if (GetSignQB()>0) for (int i=3;i--;) dircos[i] = -dircos[i]; // positive tracks move along decreasing t
924   }
925   else {
926     double gam = 1/TMath::Sqrt( fParams[kB0]*fParams[kB0] + fParams[kB1]*fParams[kB1] + 1.);
927     dircos[ fkAxID[kX] ] = fParams[kB0]*gam;
928     dircos[ fkAxID[kY] ] = fParams[kB1]*gam;
929     dircos[ fParAxis   ] = gam;
930     // decide direction
931     if (IsTypeCollision()) {
932       static double xyzF[3],xyzL[3];
933       GetPosition(xyzF,fPntFirst);
934       GetPosition(xyzL,fPntLast);
935       double dif = fCurT[fPntLast] - fCurT[fPntFirst];
936       double dr = (xyzL[kX]-xyzF[kX])*(xyzL[kX]+xyzF[kX]) + (xyzL[kY]-xyzF[kY])*(xyzL[kY]+xyzF[kY]);
937       if (dr*dif<0) for (int i=3;i--;) dircos[i] = -dircos[i]; // with increasing t the tracks comes closer to origin
938     }
939     else if (dircos[kY]>0) for (int i=3;i--;) dircos[i] = -dircos[i];  // cosmic tracks have negative angle to Y axis
940   }
941   //
942 }
943
944 //________________________________________________________________________________________________________
945 Double_t AliITSTPArrayFit::GetMachinePrec()
946 {
947   // estimate machine precision
948   Double_t eps=1.0,a;
949   do { a = 1. + (eps=eps/2.0); } while(a>1.);
950   return TMath::Abs(2.*eps);
951 }
952
953 //________________________________________________________________________________________________________
954 Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ)
955 {
956   // crude estimate of helix parameters, w/o errors and Eloss.
957   // Fast Riemann fit: Comp.Phy.Comm.131 (2000) 95
958   //
959   // if charge is not imposed (extQ==0) then it will be determined from the collision type
960   //
961   static TArrayD arrU,arrV,arrW;
962   double *parrW,*parrU,*parrV;
963   Bool_t eloss = IsELossON();
964   //
965   int np = fPntLast - fPntFirst + 1;
966   if (np<3) { AliError("At least 3 points are needed for helix fit"); return kFALSE; }
967   //
968   const float *x=fkPoints->GetX(),*y=fkPoints->GetY(),*z=fkPoints->GetZ(),*cov=fkPoints->GetCov();
969   //
970   if (fPntLast>=arrU.GetSize()) {
971     arrU.Set(2*fPntLast);
972     arrV.Set(2*fPntLast);
973     arrW.Set(2*fPntLast);
974   }
975   parrU = arrU.GetArray();
976   parrV = arrV.GetArray();
977   parrW = arrW.GetArray();
978   //
979   double uav=0,vav=0,wav=0,muu=0,muv=0,muw=0,mvv=0,mvw=0,mww=0;
980   int minRId = fPntFirst;
981   //  
982   // get points span
983   double xmn=1e9,xmx=-1e9, ymn=1e9,ymx=-1e9;
984   for (int i=fPntFirst;i<=fPntLast;i++) {
985     parrW[i] = x[i]*x[i]+y[i]*y[i];
986     if (parrW[i]<parrW[minRId]) minRId = i; // point closest to origin
987     if (xmn>x[i]) xmn = x[i];
988     if (xmx<x[i]) xmx = x[i];
989     if (ymn>y[i]) ymn = y[i];
990     if (ymx<y[i]) ymx = y[i];
991   }
992   int minRId1 = minRId>fPntFirst ? fPntFirst:fPntFirst+1;
993   for (int i=fPntFirst;i<=fPntLast;i++) if (parrW[i]<parrW[minRId1] && i!=minRId) minRId1 = i; 
994   //
995   double xshift = (xmx+xmn)/2 + 10*(ymx-ymn); // shift origin to have uniform weights
996   double yshift = (ymx+ymn)/2 - 10*(xmx-xmn);
997   //  printf("X: %+e %+e Y: %+e %+e | shift: %+e %+e\n",xmn,xmx,ymn,ymx,xshift,yshift);
998   //
999   for (int i=fPntFirst;i<=fPntLast;i++) {
1000     double xs = x[i] - xshift;
1001     double ys = y[i] - yshift;
1002     double w = xs*xs + ys*ys;
1003     double scl = 1./(1.+w);
1004     int i0 = i-fPntFirst;
1005     wav += parrW[i0] = w*scl;
1006     uav += parrU[i0] = xs*scl;
1007     vav += parrV[i0] = ys*scl;
1008   }
1009   uav /= np;    vav /= np;   wav /= np;
1010   //
1011   for (int i=fPntFirst;i<=fPntLast;i++) {
1012     //
1013     // point next to closest
1014     int i0 = i-fPntFirst;
1015     if (parrW[i0]<parrW[minRId1-fPntFirst] && i!=minRId) minRId1 = i; 
1016     double u = parrU[i0] - uav;
1017     double v = parrV[i0] - vav;
1018     double w = parrW[i0] - wav;
1019     muu += u*u;
1020     muv += u*v;
1021     muw += u*w;
1022     mvv += v*v;
1023     mvw += v*w;
1024     mww += w*w;
1025   } 
1026   muu/=np; muv/=np; muw/=np; mvv/=np; mvw/=np; mww/=np;
1027   //
1028   // find eigenvalues:
1029   double trace3 = (muu + mvv + mww)/3.;
1030   double muut = muu-trace3;
1031   double mvvt = mvv-trace3;
1032   double mwwt = mww-trace3;
1033   double q = (muut*(mvvt*mwwt-mvw*mvw) - muv*(muv*mwwt-mvw*muw) + muw*(muv*mvw-mvvt*muw))/2;
1034   double p = (muut*muut+mvvt*mvvt+mwwt*mwwt+2*(muv*muv+muw*muw+mvw*mvw))/6;
1035   double dfpp = p*p*p-q*q;
1036   dfpp = dfpp>0 ? TMath::Sqrt(dfpp)/q : 0;
1037   double ph = TMath::ATan( dfpp )/3.;
1038   if (ph<0) ph += TMath::Pi()/3;
1039   p = p>0 ? TMath::Sqrt(p) : 0;
1040   const double kSqrt3 = 1.73205080;
1041   double snp = TMath::Sin(ph);
1042   double csp = TMath::Cos(ph);
1043   //  double eg1 = trace3 + 2*p*csp;
1044   double eg2 = trace3 - p*(csp+kSqrt3*snp); // smallest one
1045   //  double eg3 = trace3 - p*(csp-kSqrt3*snp);
1046   // eigenvector for min.eigenvalue
1047   muut = muu-eg2;
1048   mvvt = mvv-eg2;
1049   mwwt = mww-eg2;
1050   double n0 = muv*mvw-muw*mvvt;
1051   double n1 = muv*muw-mvw*muut;
1052   double n2 = muut*mvvt-muv*muv;
1053   // normalize to largest one
1054   double nrm = TMath::Abs(n0);
1055   if (nrm<TMath::Abs(n1)) nrm = TMath::Abs(n1);
1056   if (nrm<TMath::Abs(n2)) nrm = TMath::Abs(n2);
1057   n0/=nrm; n1/=nrm; n2/=nrm;
1058   //
1059   double cpar = -(uav*n0 + vav*n1 + wav*n2);
1060   double xc = -n0/(cpar+n2)/2 + xshift;
1061   double yc = -n1/(cpar+n2)/2 + yshift;
1062   double rad = TMath::Sqrt(n0*n0+n1*n1-4*cpar*(cpar+n2))/2./TMath::Abs(cpar+n2);
1063   //
1064   //  printf("Rad: %+e xc: %+e yc: %+e | X0: %+e Y0: %+e | X1: %+e Y1: %+e\n",rad,xc,yc, x[minRId],y[minRId],x[minRId1],y[minRId1]);
1065
1066   // linear circle fit --------------------------------------------------- <<<
1067   //
1068   // decide sign(Q*B) and fill cicrle parameters ------------------------- >>>
1069   int sqb;
1070   if (extQ) {
1071     SetCharge(extQ); 
1072     sqb = fBz<0 ? -GetCharge():GetCharge();
1073   }
1074   else { 
1075     // determine the charge from the collision type and field sign
1076     // the negative Q*B will have positive Vc x dir product Z component
1077     // with Vc={-xc,-yc} : vector from circle center to the origin
1078     // and V0 - track direction vector (take {0,-1,1} for cosmics)
1079     // If Bz is not provided, assume positive Bz
1080     if ( IsTypeCosmics() ) sqb = xc>0 ? -1:1;
1081     else {
1082       // track direction vector as a - diference between the closest and the next to closest to origin points
1083       double v0X = x[minRId1] - x[minRId];
1084       double v0Y = y[minRId1] - y[minRId];
1085       sqb = (yc*v0X - xc*v0Y)>0 ? -1:1;
1086     }
1087     SetCharge( fBz<0 ? -sqb : sqb);
1088   }
1089   //
1090   Double_t phi = TMath::ATan2(yc,xc);
1091   if (sqb<0) phi += TMath::Pi();
1092   if      (phi > TMath::Pi()) phi -= 2.*TMath::Pi();
1093   else if (phi <-TMath::Pi()) phi += 2.*TMath::Pi();
1094   fParams[kPhi0] = phi;  
1095   fParams[kR0]   = sqb<0 ? -rad:rad;  
1096   fParams[kD0] = xc*TMath::Cos(phi) + yc*TMath::Sin(phi) - fParams[kR0];
1097   //
1098   // decide sign(Q*B) and fill cicrle parameters ------------------------- <<<
1099   //
1100   // find z-offset and dip + the parameter t of closest approach to hits - >>>
1101   //
1102   UInt_t hitLrPos=0;  // pattern of hit layers at pos
1103   UInt_t hitLrNeg=0;  // and negative t's
1104
1105   Double_t ss=0,st=0,sz=0,stt=0,szt=0;
1106   for (int i=fPntFirst;i<=fPntLast;i++) {
1107     //
1108     Double_t ze2 = cov[i*6 + kZZ];
1109     Double_t t = TMath::ATan2(yc-y[i],xc-x[i]) - fParams[kPhi0]; // angle at measured z
1110     if (fParams[kR0]<0)  t += TMath::Pi();
1111     if      (t > TMath::Pi()) t -= TMath::Pi()*2;
1112     else if (t <-TMath::Pi()) t += TMath::Pi()*2;
1113     if (ze2<fgkAlmostZero) ze2 = 1E-8;
1114     ze2 = 1./ze2;
1115     ss += ze2;
1116     st += t*ze2;
1117     stt+= t*t*ze2;
1118     sz += z[i]*ze2;
1119     szt+= z[i]*t*ze2;
1120     //
1121     fCurT[i] = t; // parameter of the closest approach to the point
1122     //    printf("%d %+e %+e %+e %+e\n",i,x[i],y[i],z[i],t);
1123     if (eloss) {
1124       double r = TMath::Sqrt(x[i]*x[i]+y[i]*y[i]);
1125       int lr;
1126       for (lr=kMaxLrITS;lr--;) if ( IsZero(r-fgkRLayITS[ lr ],1.) ) break;
1127       if (lr<kMaxLrITS) {
1128         if (t>0) hitLrPos |= (1<<lr);  // set bit of the layer
1129         else     hitLrNeg |= (1<<lr);  // set bit of the layer
1130       }
1131     }
1132   }
1133   double det = ss*stt - st*st;
1134   if (TMath::Abs(det)<fgkAlmostZero) { // no Z dependence
1135     fParams[kDZ]  = sz/ss;
1136     fParams[kDip] = 0;
1137   }
1138   else {
1139     fParams[kDZ]  =  (sz*stt-st*szt)/det;
1140     fParams[kDip] = -(ss*szt-st*sz)/det/fParams[kR0];
1141   }
1142   //
1143   // find z-offset and dip + the parameter t of closest approach to hits - <<<
1144   //
1145   // fill info needed to account for ELoss ------------------------------- >>>
1146   if (eloss) {
1147     fNElsPnt = fPntLast - fPntFirst + 1;
1148     //
1149     // to account for the energy loss in the passive volumes, calculate the relevant t-parameters 
1150     double* tcur = fCurT + fPntFirst;
1151     double* ecur = fElsDR+ fPntFirst;
1152     //
1153     for (int ilp=3;ilp--;) {
1154       int id = fgkPassivLrITS[ilp];
1155       double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1156       if (tp<0) continue; // does not hit this radius
1157       //
1158       tcur[fNElsPnt] = GetSignQB()>0 ? -tp : tp;
1159       ecur[fNElsPnt] = fgRhoLITS[ id ];
1160       fNElsPnt++;
1161       //      printf("Passive  on lr %d  %+e\n",ilp,tcur[fNElsPnt-1]);
1162       //
1163       if (IsTypeCosmics() && !IsZero(tp)) { // 2 crossings for cosmics
1164         tcur[fNElsPnt] = -tcur[fNElsPnt-1];
1165         ecur[fNElsPnt] =  ecur[fNElsPnt-1];
1166         fNElsPnt++;
1167         //printf("Passive* on lr %d  %+e\n",ilp,-tcur[fNElsPnt-1]);
1168       }
1169       //
1170     }
1171     // check if some active layers did not miss the hit, treat them as passive
1172     for (int ilp=6;ilp--;) {
1173       int id = fgkActiveLrITS[ilp];
1174       double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1175       if (tp<0) continue; // does not hit this radius
1176       //
1177       if ( (GetSignQB()>0||IsTypeCosmics()) && !(hitLrNeg & (1<<id)) ) {
1178         tcur[fNElsPnt] = -tp;
1179         ecur[fNElsPnt] = fgRhoLITS[ id ];
1180         fNElsPnt++;
1181         //printf("Missed  on lr %d  %+e\n",ilp,-tp);
1182       }
1183       //
1184       if ( (GetSignQB()<0||IsTypeCosmics()) && !(hitLrPos & (1<<id)) ) {
1185         tcur[fNElsPnt] = tp;
1186         ecur[fNElsPnt] = fgRhoLITS[ id ];
1187         fNElsPnt++;
1188         //printf("Missed* on lr %d  %e\n",ilp,tp);
1189       }
1190     }
1191     //
1192     TMath::Sort(fNElsPnt,fCurT+fPntFirst,fElsId,kFALSE);    // index e-loss points in increasing order
1193     // find the position of smallest positive t-param
1194     for (fFirstPosT=0;fFirstPosT<fNElsPnt;fFirstPosT++) if (fCurT[ fElsId[ fFirstPosT ] ]>0) break;
1195     //
1196     Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]);
1197     Double_t ptot = TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum and energy
1198     Double_t etot = TMath::Sqrt(ptot*ptot + fMass*fMass);      // in the point of closest approach to beam
1199     Double_t normS[3];
1200     //
1201     // Positive t-params: along the track direction for negative track, against for positive
1202     //   Double_t pcur = ptot, ecurr = etot;
1203     for (int ip=fFirstPosT;ip<fNElsPnt;ip++) {
1204       int tID = fElsId[ip];
1205       Double_t t = fCurT[ tID ];
1206       //
1207       if (tID>fPntLast) { // this is not a hit layer but passive layer
1208         double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1209                                   xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1210         normS[0] = -TMath::Cos(php);  // normal to the cylinder at intersection point
1211         normS[1] = -TMath::Sin(php);
1212         normS[2] = 0;
1213       }
1214       else GetNormal(normS,fkPoints->GetCov()+tID*6);   // vector normal to hit module
1215       fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1216     }
1217     //
1218     // negaive t-params: against the track direction for negative track, along for positive
1219     //    pcur  = ptot;
1220     //    ecurr = etot;
1221     for (int ip=fFirstPosT;ip--;) {
1222       int tID = fElsId[ip];
1223       Double_t t = fCurT[ tID ];
1224       //
1225       if (tID>fPntLast) { // this is not a hit layer but passive layer
1226         double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1227                                   xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1228         normS[0] = -TMath::Cos(php);  // normal to the cylinder at intersection point
1229         normS[1] = -TMath::Sin(php);
1230         normS[2] = 0;   
1231       }
1232       else GetNormal(normS,fkPoints->GetCov()+tID*6);   // vector normal to hit module
1233       //
1234       fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1235     }
1236   }
1237   // fill info needed to account for ELoss ------------------------------- <<<
1238   //
1239   return kTRUE;
1240 }
1241
1242 /*
1243 //________________________________________________________________________________________________________
1244 Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ)
1245 {
1246   // crude estimate of helix parameters, w/o errors and Eloss.
1247   // 1st fit the circle (R,xc,yc) by minimizing
1248   // chi2 = sum{ (bx*xi + by*yi + xi^2+yi^2 + rho)^2 } vs bx,by,rho
1249   // with bx = -2*xc, by = -2*yc , rho = xc^2+yc^2 - R2
1250   //
1251   // if charge is not imposed (extQ==0) then it will be determined from the collision type
1252   //
1253   Bool_t eloss = IsELossON();
1254   //
1255   int np = fPntLast - fPntFirst + 1;
1256   if (np<3) { AliError("At least 3 points are needed for helix fit"); return kFALSE; }
1257   //
1258   const float *x=fkPoints->GetX(),*y=fkPoints->GetY(),*z=fkPoints->GetZ(),*cov=fkPoints->GetCov();
1259   //
1260   // linear circle fit --------------------------------------------------- >>>
1261   Double_t sxx=0,sxy=0,syy=0,sx=0,sy=0,rhs0=0,rhs1=0,rhs2=0,minR=1E9;
1262   int minRId = 0;
1263   for (int i=fPntFirst;i<=fPntLast;i++) {
1264     Double_t xx = x[i]*x[i];
1265     Double_t yy = y[i]*y[i];
1266     Double_t xy = x[i]*y[i];
1267     Double_t xxyy = xx + yy;
1268     //
1269     sxx += xx;
1270     sxy += xy;
1271     syy += yy;
1272     sx  += x[i];
1273     sy  += y[i];
1274     //
1275     rhs0 -= xxyy*x[i];
1276     rhs1 -= xxyy*y[i];
1277     rhs2 -= xxyy;
1278     // 
1279     // remember Id of the point closest to origin, to determine the charge  
1280     if (xxyy<minR) { minR   = xxyy; minRId = i; }
1281     //
1282     if (eloss) { // find layer id
1283       int lrid,volid = fkPoints->GetVolumeID()[i];
1284       if (volid>0) lrid = fgkActiveLrITS[AliGeomManager::VolUIDToLayer(fkPoints->GetVolumeID()[i])-1];
1285       else { // missing layer info, find from radius
1286         double r = TMath::Sqrt(xxyy);
1287         for (lrid=kMaxLrITS;lrid--;) if ( IsZero(r-fgkRLayITS[ lrid ],1.) ) break;
1288       }
1289       fElsDR[i] = (lrid>=0 && lrid<kMaxLrITS) ? fgRhoLITS[ lrid ] : 0;   // eloss for normal track
1290     }
1291     //
1292   }
1293   //
1294   Double_t mn00 = syy*np-sy*sy;
1295   Double_t mn01 = sxy*np-sy*sx;
1296   Double_t mn02 = sxy*sy-syy*sx;
1297   Double_t det  = sxx*mn00 - sxy*mn01 + sx*mn02; 
1298   if (TMath::Abs(det)<fgkAlmostZero) return kFALSE;
1299   //
1300   Double_t mn11 = sxx*np-sx*sx;
1301   Double_t mn12 = sxx*sy-sxy*sx;
1302   Double_t mn22 = sxx*syy-sxy*sxy;
1303   //
1304   Double_t mi00 =  mn00/det;
1305   Double_t mi01 = -mn01/det;
1306   Double_t mi02 =  mn02/det;
1307   Double_t mi11 =  mn11/det;
1308   Double_t mi12 = -mn12/det;
1309   Double_t mi22 =  mn22/det;
1310   //
1311   Double_t xc   = -(rhs0*mi00 + rhs1*mi01 + rhs2*mi02)/2;
1312   Double_t yc   = -(rhs0*mi01 + rhs1*mi11 + rhs2*mi12)/2;
1313   Double_t rho2 =  (rhs0*mi02 + rhs1*mi12 + rhs2*mi22);
1314
1315   //
1316   // check
1317   double bx = -2*xc;
1318   double by = -2*yc;
1319   double sm0=0,sm1=0;
1320   for (int i=fPntFirst;i<=fPntLast;i++) {
1321     double dst = bx*x[i]+by*y[i]+x[i]*x[i]+y[i]*y[i]+rho2;
1322     sm0 += dst;
1323     sm1 += dst*dst;
1324     printf("Point %d: %+e %+e %+e\n",i,dst,sm0,sm1);
1325   }
1326
1327   //
1328   Double_t rad = xc*xc + yc*yc - rho2;
1329   rad = (rad>fgkAlmostZero) ? (TMath::Sqrt(rad)):fgkAlmostZero;
1330   //
1331   //  printf("Rad: %+e xc: %+e yc: %+e\n",rad,xc,yc);
1332
1333   // linear circle fit --------------------------------------------------- <<<
1334   //
1335   // decide sign(Q*B) and fill cicrle parameters ------------------------- >>>
1336   int sqb;
1337   if (extQ) {
1338     SetCharge(extQ); 
1339     sqb = fBz<0 ? -GetCharge():GetCharge();
1340   }
1341   else { 
1342     // determine the charge from the collision type and field sign
1343     // the negative Q*B will have positive Vc x V0 product Z component
1344     // with Vc={-xc,-yc} : vector from circle center to the origin
1345     // and V0 - track direction vector (take {0,-1,1} for cosmics)
1346     // If Bz is not provided, assume positive Bz
1347     sqb = ( IsTypeCosmics() ? xc:(yc*x[minRId]-xc*y[minRId]) ) > 0 ? -1:1;
1348     SetCharge( fBz<0 ? -sqb : sqb);
1349   }
1350   //
1351   Double_t phi = TMath::ATan2(yc,xc);
1352   if (sqb<0) phi += TMath::Pi();
1353   if      (phi > TMath::Pi()) phi -= 2.*TMath::Pi();
1354   else if (phi <-TMath::Pi()) phi += 2.*TMath::Pi();
1355   fParams[kPhi0] = phi;  
1356   fParams[kR0]   = sqb<0 ? -rad:rad;  
1357   fParams[kD0] = xc*TMath::Cos(phi) + yc*TMath::Sin(phi) - fParams[kR0];
1358   //
1359   // decide sign(Q*B) and fill cicrle parameters ------------------------- <<<
1360   //
1361   // find z-offset and dip + the parameter t of closest approach to hits - >>>
1362   //
1363   UInt_t hitLrPos=0;  // pattern of hit layers at pos
1364   UInt_t hitLrNeg=0;  // and negative t's
1365
1366   Double_t ss=0,st=0,sz=0,stt=0,szt=0;
1367   for (int i=fPntFirst;i<=fPntLast;i++) {
1368     //
1369     Double_t ze2 = cov[i*6 + kZZ];
1370     Double_t t = TMath::ATan2(yc-y[i],xc-x[i]) - fParams[kPhi0]; // angle at measured z
1371     if (fParams[kR0]<0)  t += TMath::Pi();
1372     if      (t > TMath::Pi()) t -= TMath::Pi()*2;
1373     else if (t <-TMath::Pi()) t += TMath::Pi()*2;
1374     if (ze2<fgkAlmostZero) ze2 = 1E-8;
1375     ze2 = 1./ze2;
1376     ss += ze2;
1377     st += t*ze2;
1378     stt+= t*t*ze2;
1379     sz += z[i]*ze2;
1380     szt+= z[i]*t*ze2;
1381     //
1382     fCurT[i] = t; // parameter of the closest approach to the point
1383     //    printf("%d %+e %+e %+e %+e\n",i,x[i],y[i],z[i],t);
1384     if (eloss) {
1385       double r = TMath::Sqrt(x[i]*x[i]+y[i]*y[i]);
1386       int lr;
1387       for (lr=kMaxLrITS;lr--;) if ( IsZero(r-fgkRLayITS[ lr ],1.) ) break;
1388       if (lr<kMaxLrITS) {
1389         if (t>0) hitLrPos |= (1<<lr);  // set bit of the layer
1390         else     hitLrNeg |= (1<<lr);  // set bit of the layer
1391       }
1392     }
1393   }
1394   det = ss*stt - st*st;
1395   if (TMath::Abs(det)<fgkAlmostZero) { // no Z dependence
1396     fParams[kDZ]  = sz/ss;
1397     fParams[kDip] = 0;
1398   }
1399   else {
1400     fParams[kDZ]  =  (sz*stt-st*szt)/det;
1401     fParams[kDip] = -(ss*szt-st*sz)/det/fParams[kR0];
1402   }
1403   //
1404   // find z-offset and dip + the parameter t of closest approach to hits - <<<
1405   //
1406   // fill info needed to account for ELoss ------------------------------- >>>
1407   if (eloss) {
1408     fNElsPnt = fPntLast - fPntFirst + 1;
1409     //
1410     // to account for the energy loss in the passive volumes, calculate the relevant t-parameters 
1411     double* tcur = fCurT + fPntFirst;
1412     double* ecur = fElsDR+ fPntFirst;
1413     //
1414     for (int ilp=3;ilp--;) {
1415       int id = fgkPassivLrITS[ilp];
1416       double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1417       if (tp<0) continue; // does not hit this radius
1418       //
1419       tcur[fNElsPnt] = GetSignQB()>0 ? -tp : tp;
1420       ecur[fNElsPnt] = fgRhoLITS[ id ];
1421       fNElsPnt++;
1422       //      printf("Passive  on lr %d  %+e\n",ilp,tcur[fNElsPnt-1]);
1423       //
1424       if (IsTypeCosmics() && !IsZero(tp)) { // 2 crossings for cosmics
1425         tcur[fNElsPnt] = -tcur[fNElsPnt-1];
1426         ecur[fNElsPnt] =  ecur[fNElsPnt-1];
1427         fNElsPnt++;
1428         //printf("Passive* on lr %d  %+e\n",ilp,-tcur[fNElsPnt-1]);
1429       }
1430       //
1431     }
1432     // check if some active layers did not miss the hit, treat them as passive
1433     for (int ilp=6;ilp--;) {
1434       int id = fgkActiveLrITS[ilp];
1435       double tp = GetHelixParAtR( fgkRLayITS[ id ] );
1436       if (tp<0) continue; // does not hit this radius
1437       //
1438       if ( (GetSignQB()>0||IsTypeCosmics()) && !(hitLrNeg & (1<<id)) ) {
1439         tcur[fNElsPnt] = -tp;
1440         ecur[fNElsPnt] = fgRhoLITS[ id ];
1441         fNElsPnt++;
1442         //printf("Missed  on lr %d  %+e\n",ilp,-tp);
1443       }
1444       //
1445       if ( (GetSignQB()<0||IsTypeCosmics()) && !(hitLrPos & (1<<id)) ) {
1446         tcur[fNElsPnt] = tp;
1447         ecur[fNElsPnt] = fgRhoLITS[ id ];
1448         fNElsPnt++;
1449         //printf("Missed* on lr %d  %e\n",ilp,tp);
1450       }
1451     }
1452     //
1453     TMath::Sort(fNElsPnt,fCurT+fPntFirst,fElsId,kFALSE);    // index e-loss points in increasing order
1454     // find the position of smallest positive t-param
1455     for (fFirstPosT=0;fFirstPosT<fNElsPnt;fFirstPosT++) if (fCurT[ fElsId[ fFirstPosT ] ]>0) break;
1456     //
1457     Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]);
1458     Double_t ptot = TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum and energy
1459     Double_t etot = TMath::Sqrt(ptot*ptot + fMass*fMass);      // in the point of closest approach to beam
1460     Double_t normS[3];
1461     //
1462     // Positive t-params: along the track direction for negative track, against for positive
1463     Double_t pcur = ptot, ecurr = etot;
1464     for (int ip=fFirstPosT;ip<fNElsPnt;ip++) {
1465       int tID = fElsId[ip];
1466       Double_t t = fCurT[ tID ];
1467       //
1468       if (tID>fPntLast) { // this is not a hit layer but passive layer
1469         double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1470                                   xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1471         normS[0] = -TMath::Cos(php);  // normal to the cylinder at intersection point
1472         normS[1] = -TMath::Sin(php);
1473         normS[2] = 0;
1474       }
1475       else GetNormal(normS,fkPoints->GetCov()+tID*6);   // vector normal to hit module
1476       fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1477     }
1478     //
1479     // negaive t-params: against the track direction for negative track, along for positive
1480     pcur  = ptot;
1481     ecurr = etot;
1482     for (int ip=fFirstPosT;ip--;) {
1483       int tID = fElsId[ip];
1484       Double_t t = fCurT[ tID ];
1485       //
1486       if (tID>fPntLast) { // this is not a hit layer but passive layer
1487         double php = TMath::ATan2(yc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t),
1488                                   xc-fParams[kR0]*TMath::Cos(fParams[kPhi0]+t));
1489         normS[0] = -TMath::Cos(php);  // normal to the cylinder at intersection point
1490         normS[1] = -TMath::Sin(php);
1491         normS[2] = 0;   
1492       }
1493       else GetNormal(normS,fkPoints->GetCov()+tID*6);   // vector normal to hit module
1494       //
1495       fElsDR[tID] = GetDRofELoss(t,cdip,fElsDR[tID],normS,ptot,etot);
1496     }
1497   }
1498   // fill info needed to account for ELoss ------------------------------- <<<
1499   //
1500   return kTRUE;
1501 }
1502 */
1503 //____________________________________________________
1504 Double_t AliITSTPArrayFit::FitHelix(Int_t extQ, Double_t extPT,Double_t extPTerr)
1505 {
1506   // fit by helix accounting for the errors of all coordinates (and energy loss if requested)
1507   // 
1508   // If extQ is non-0, its sign is imposed as a charge of the track
1509   // If extPT>0 and extPTerr>=0, constrain to measured tr.momentum PT 
1510   // with corresponding error (err=0 -> rel.err=1e-6)
1511   //
1512   double chiprev=1e99;
1513   //const Double_t kMaxTEffect = 1E-6;
1514   Double_t dXYZdGlo[3*5],dXYZdLoc[3],xyzRes[3];
1515   //
1516   SetFitDone(kFALSE);
1517   fChi2NDF = -1;
1518   //
1519   if (!FitHelixCrude(extQ)) return -1; // get initial estimate, ignoring the errors
1520   //
1521   if (TMath::Abs(fParams[kR0])>fMaxRforHelix && extPT<=0) {
1522     fSwitch2Line = kTRUE;
1523     return FitLine();
1524   }
1525   //
1526   //
1527   if (!IsCovInv()) InvertPointsCovMat();    // prepare inverted errors
1528   if (!fParSol) fParSol = new AliParamSolver(5);
1529   fParSol->SetNGlobal(5);
1530   //
1531   //  printf("-1 | %+.2e %+.2e %+.2e %+.2e %+.2e | chi2: %+.4e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4],CalcChi2NDF());
1532   int iter = 0;
1533   fChi2NDF = 1e99;
1534   Bool_t converged = kFALSE;
1535   while(iter<fMaxIter) {
1536     chiprev = fChi2NDF;
1537     fParSol->Clear();
1538     for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1539       //
1540       GetResiduals(xyzRes, ip); // current residuals at point ip
1541       Double_t rrho = fParams[kR0]+fParams[kD0];
1542       Double_t cs0  = TMath::Cos(fParams[kPhi0]);
1543       Double_t sn0  = TMath::Sin(fParams[kPhi0]);
1544       Double_t cst  = TMath::Cos(fCurT[ip]+fParams[kPhi0]);
1545       Double_t snt  = TMath::Sin(fCurT[ip]+fParams[kPhi0]);
1546       //
1547       int offs = kD0;                  // dXYZ/dD0
1548       dXYZdGlo[offs + kX] = cs0;
1549       dXYZdGlo[offs + kY] = sn0;
1550       dXYZdGlo[offs + kZ] = 0;
1551       //
1552       offs = kPhi0*3;                  // dXYZ/dPhi0
1553       dXYZdGlo[offs + kX] = -rrho*sn0 + fParams[kR0]*snt;
1554       dXYZdGlo[offs + kY] =  rrho*cs0 - fParams[kR0]*cst;
1555       dXYZdGlo[offs + kZ] = 0;
1556       //
1557       offs = kR0*3;                   // dXYZ/dR0
1558       if (TMath::Abs(fParams[kR0])<0.9*fMaxRforHelix) {
1559         dXYZdGlo[offs + kX] =  cs0 - cst;
1560         dXYZdGlo[offs + kY] =  sn0 - snt;
1561         dXYZdGlo[offs + kZ] =  -fParams[kDip]*fCurT[ip];
1562       }
1563       else {
1564         dXYZdGlo[offs + kX] = dXYZdGlo[offs + kY] = dXYZdGlo[offs + kZ] = 0;
1565         fParSol->AddConstraint(kR0,0,1.e2);
1566       }
1567       //
1568       offs = kDZ*3;                   // dXYZ/dDZ
1569       dXYZdGlo[offs + kX] =  0;
1570       dXYZdGlo[offs + kY] =  0;
1571       dXYZdGlo[offs + kZ] =  1.;
1572       //
1573       offs = kDip*3;                  // dXYZ/dDip
1574       dXYZdGlo[offs + kX] =  0;
1575       dXYZdGlo[offs + kY] =  0;
1576       dXYZdGlo[offs + kZ] = -fParams[kR0]*fCurT[ip];
1577       //
1578       //      /*
1579       dXYZdLoc[kX] =  fParams[kR0]*snt;
1580       dXYZdLoc[kY] = -fParams[kR0]*cst;
1581       dXYZdLoc[kZ] = -fParams[kR0]*fParams[kDip];
1582       //      */
1583       //      dXYZdLoc[0] = dXYZdLoc[1] = dXYZdLoc[2] = 0;
1584       //
1585       fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip),GetCovIScale(ip));
1586     }
1587     //
1588     if (extPT>0) { // add constraint on pt
1589       if (extPTerr<fgkAlmostZero) extPTerr = 1e-6*extPT;
1590       Double_t cf = fBz*GetCharge()*fgkCQConv;
1591       Double_t err2i = extPTerr/cf;
1592       err2i = 1./err2i/err2i;
1593       //      printf("Constrain R to %+e\n",extPT/cf);
1594       fParSol->AddConstraint(kR0,-extPT/cf+fParams[kR0],err2i);
1595     }
1596     if (!fParSol->Solve()) { AliError("Failed to fit helix"); return -1; }
1597     Double_t *deltaG = fParSol->GetGlobals();
1598     //    Double_t *deltaT = fParSol->GetLocals();
1599     for (int ipar=5;ipar--;) fParams[ipar] -= deltaG[ipar];
1600     //
1601     if (TMath::Abs(fParams[kR0])>0.9*fMaxRforHelix) fParams[kR0] = TMath::Sign(0.9*fMaxRforHelix,fParams[kR0]);
1602     //
1603     for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1604       fCurT[ip] = CalcParPCA(ip);
1605       //      printf("iter%d : delta%2d : expl: %+e | %+e | R=%+e p0=%+e\n",iter,ip,deltaT[ip-fPntFirst], fCurT[ip],fParams[kR0],fParams[kPhi0]);
1606       //      fCurT[ip] -= deltaT[ip-fPntFirst];
1607     }
1608     iter++;
1609     //
1610     fChi2NDF = CalcChi2NDF();
1611     //        printf("%2d | %+.2e %+.2e %+.2e %+.2e %+.2e | chi2: %+.4e %+.4e\n",iter,deltaG[0],deltaG[1],deltaG[2],deltaG[3],deltaG[4],fChi2NDF,fChi2NDF-chiprev);
1612     //        printf("->>  %+.2e %+.2e %+.2e %+.2e %+.2e | Chi2: %+.6e %+.6e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4],fChi2NDF,fChi2NDF-chiprev);
1613     double difchi2 = chiprev - fChi2NDF;
1614     if ( difchi2<fEps && TMath::Abs(difchi2)<1e-4) {converged = kTRUE; break;} 
1615     //    if (errT*TMath::Abs(fParams[kR0])<kMaxTEffect && errP<fEps) {converged = kTRUE; break;} 
1616   }
1617   //
1618   if (!converged) {
1619     AliDebug(2,Form("Max number of %d iteration reached, Current chi2:%.3e, chi2 change %+.3e",iter,
1620                     fChi2NDF,chiprev-fChi2NDF));
1621     for (int ip=fPntFirst;ip<=fPntLast;ip++)
1622       AliDebug(2,Form("P%2d| %+.3e %+.3e %+.3e\n",ip,fkPoints->GetX()[ip],fkPoints->GetY()[ip],fkPoints->GetZ()[ip]));
1623
1624   }
1625   fIter = iter;
1626   SetCharge( fParams[kR0]>0 ? (fBz<0?-1:1):(fBz>0?-1:1) );
1627   SetFitDone();
1628   //  printf("F1>> %+.7e %+.7e %+.7e %+.7e %.7e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4]);
1629   //
1630   return fChi2NDF;
1631 }
1632
1633 //____________________________________________________
1634 Double_t AliITSTPArrayFit::FitLine()
1635 {
1636   // fit by helix accounting for the errors of all coordinates (and energy loss if requested)
1637   // 
1638   double chiprev=1e99;
1639   //  const Double_t kMaxTEffect = 1.e-6;
1640   Double_t dXYZdGlo[3*4],dXYZdLoc[3],xyzRes[3];
1641   SetFitDone(kFALSE);
1642   fChi2NDF = -1;
1643   //
1644   if (fParAxis<0) SetParAxis(ChoseParAxis());
1645   //
1646   const float *xyzp[3]={fkPoints->GetX(),fkPoints->GetY(),fkPoints->GetZ()};
1647   if (!IsCovInv()) InvertPointsCovMat();
1648   if (!FitLineCrude()) return -1; // get initial estimate, ignoring the errors
1649   //
1650   if (!fParSol) fParSol = new AliParamSolver(5);
1651   fParSol->SetNGlobal(4);
1652   // initial set of parameters
1653   for (int ip=fPntFirst;ip<=fPntLast;ip++) fCurT[ip] = xyzp[fParAxis][ip]; // use measured param-coordinate
1654   //
1655   int iter = 0;
1656   Bool_t converged = kFALSE;
1657   fChi2NDF = 1e99;
1658   while(iter<fMaxIter) {
1659     chiprev = fChi2NDF;
1660     fParSol->Clear();
1661     for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1662       //
1663       int offs;
1664       GetResiduals(xyzRes, ip); // current residuals at point ip
1665       //
1666       offs = kA0*3;                   // dXYZ/dA0
1667       dXYZdGlo[offs + fkAxID[kX]] = 1;
1668       dXYZdGlo[offs + fkAxID[kY]] = 0;
1669       dXYZdGlo[offs + fParAxis  ] = 0;
1670       //
1671       offs = kB0*3;                   // dXYZ/dB0
1672       dXYZdGlo[offs + fkAxID[kX]] = fCurT[ip];
1673       dXYZdGlo[offs + fkAxID[kY]] = 0;
1674       dXYZdGlo[offs + fParAxis  ] = 0;
1675       //
1676       offs = kA1*3;                   // dXYZ/dA1
1677       dXYZdGlo[offs + fkAxID[kX]] = 0;
1678       dXYZdGlo[offs + fkAxID[kY]] = 1;
1679       dXYZdGlo[offs + fParAxis  ] = 0;
1680       //
1681       offs = kB1*3;                   // dXYZ/dB1
1682       dXYZdGlo[offs + fkAxID[kX]] = 0;
1683       dXYZdGlo[offs + fkAxID[kY]] = fCurT[ip];
1684       dXYZdGlo[offs + fParAxis  ] = 0;
1685       //
1686       dXYZdLoc[ fkAxID[kX] ] =  fParams[kB0];  // dX/dt
1687       dXYZdLoc[ fkAxID[kY] ] =  fParams[kB1];  // dY/dt
1688       dXYZdLoc[ fParAxis   ] =  1;
1689       //
1690       fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip),GetCovIScale(ip));
1691     }
1692     //
1693     if (!fParSol->Solve()) { AliError("Failed to fit line"); return -1; }
1694     Double_t *deltaG = fParSol->GetGlobals();
1695     Double_t *deltaT = fParSol->GetLocals();
1696     for (int ipar=4;ipar--;) fParams[ipar] -= deltaG[ipar];
1697     for (int ip=fPntFirst;ip<=fPntLast;ip++) fCurT[ip] -= deltaT[ip-fPntFirst];
1698     iter++;
1699     fChi2NDF = CalcChi2NDF();
1700     //    printf("%d %+e %+e | %+.2e %+.2e %+.2e %+.2e | chi2: %+.4e %+.4e\n",iter,errP,errT, deltaG[0],deltaG[1],deltaG[2],deltaG[3],fChi2NDF,fChi2NDF-chiprev);
1701     //    printf("->> %+.2e %+.2e %+.2e %+.2e %+.2e | Chi2: %+.6e %+.6e\n",fParams[0],fParams[1],fParams[2],fParams[3],fParams[4],fChi2NDF,fChi2NDF-chiprev);
1702     double difchi2 = chiprev - fChi2NDF;
1703     if ( difchi2<fEps && TMath::Abs(difchi2)<1e-4) {converged = kTRUE; break;} 
1704     chiprev = fChi2NDF;
1705     //    if (errT<kMaxTEffect && errP<fEps) {converged = kTRUE; break;} 
1706   }
1707   //
1708   if (!converged) {
1709     AliDebug(2,Form("Max number of %d iteration reached, Current chi2:%.3e, chi2 change %+.3e",iter,
1710                     fChi2NDF,chiprev-fChi2NDF));
1711     for (int ip=fPntFirst;ip<=fPntLast;ip++)
1712       AliDebug(2,Form("P%2d| %+.3e %+.3e %+.3e\n",ip,fkPoints->GetX()[ip],fkPoints->GetY()[ip],fkPoints->GetZ()[ip]));
1713   }
1714   fIter = iter;
1715   SetFitDone();
1716   //printf("F1>> %+.2e %+.2e %+.2e %+.2e\n",fParams[0],fParams[1],fParams[2],fParams[3]);
1717   return fChi2NDF;
1718   //
1719 }
1720
1721 //____________________________________________________
1722 void AliITSTPArrayFit::GetNormal(Double_t *norm, const Float_t *covMat) 
1723 {
1724   // obtain the lab normal vector to the sensor from the covariance matrix
1725   // in such a way that when the local frame of the sensor coincides with 
1726   // the lab frame, the vector {0,1,0} is obtained
1727   Double_t tgxy = TMath::Tan(0.5*TMath::ATan2(2.*covMat[kXY],covMat[kYY]-covMat[kXX]));
1728   Double_t tgyz = TMath::Tan(0.5*TMath::ATan2(2.*covMat[kYZ],covMat[kZZ]-covMat[kYY]));
1729   norm[kY] = 1./TMath::Sqrt(1 + tgxy*tgxy + tgyz*tgyz);
1730   norm[kX] = norm[kY]*tgxy;
1731   norm[kZ] = norm[kY]*tgyz;
1732   //
1733 }
1734
1735 //____________________________________________________
1736 Double_t AliITSTPArrayFit::GetDRofELoss(Double_t t,Double_t cdip,Double_t rhoL,const Double_t *normS, 
1737                                         Double_t &p,Double_t &e) const
1738 {
1739   // Calculate energy loss of the particle at given t-param on the layer with rhoL (thickness*density) with
1740   // normal vector normS in the lab. The particle before eloss has energy "e" and momentum "p"
1741   // cdip = cosine of the dip angle = 1/sqrt(1+tgL^2)
1742   // Return the change DR of the radius due to the ELoss 
1743   //
1744   // NOTE: with B>0 the negative particles propagate along increasing t-param and positive 
1745   // particles - against.
1746   // t-param = 0 corresponds to the point of closest approach of the track to the beam.
1747   // Since the fitted helix parameters of the track are defined in this PCA point, when the correction
1748   // is applied upstream of the PCS, the energy must be increased (DR>0) rather than decreased (DR<0)
1749   //
1750   Double_t dirTr[3];
1751   dirTr[0] = -TMath::Sin(fParams[kPhi0]+t);
1752   dirTr[1] =  TMath::Cos(fParams[kPhi0]+t);
1753   dirTr[2] =  fParams[kDip];
1754   // cosine of the impact angle
1755   Double_t cosImp = cdip*TMath::Abs(dirTr[0]*normS[0]+dirTr[1]*normS[1]+dirTr[2]*normS[2]);
1756   //
1757   if (cosImp<0.3) cosImp = 0.3; //?
1758   Double_t dE = AliExternalTrackParam::BetheBlochSolid(p/fMass)*rhoL/cosImp;
1759   Double_t dP = e/p*dE;
1760   //
1761   if (t*GetSignQB() < 0) {
1762     dP =  -dP;
1763     dE =  -dE;
1764   }
1765   //
1766   if (p+dP<0) {
1767     AliInfo(Form("Estimated PLoss %.3f is larger than particle momentum %.3f. Skipping",dP,p));
1768     return 0;
1769   }
1770   //
1771   p += dP;
1772   e += dE;
1773   //
1774   return fCharge*dP*cdip/fBz/fgkCQConv;
1775 }
1776
1777 //_____________________________________________________________
1778 Double_t AliITSTPArrayFit::GetLineOffset(Int_t axis) const
1779 {
1780   // return intercept of the parameterization coord = intercept + slope*t for given axis
1781   if (fParAxis<0) return -1E6; // no line fit
1782   if (axis==fParAxis) return 0;
1783   if (fParAxis==kX) return fParams[axis==kY ? kA0 : kA1 ];
1784   if (fParAxis==kY) return fParams[axis==kZ ? kA0 : kA1 ];
1785   return fParams[axis==kX ? kA0 : kA1 ];
1786 }
1787
1788 //_____________________________________________________________
1789 Double_t AliITSTPArrayFit::GetLineSlope(Int_t axis) const
1790 {
1791   // return intercept of the parameterization coord = intercept + slope*t for given axis
1792   if (fParAxis<0) return -1E6; // no line fit
1793   if (axis==fParAxis) return 1.;
1794   if (fParAxis==kX) return fParams[axis==kY ? kB0 : kB1 ];
1795   if (fParAxis==kY) return fParams[axis==kZ ? kB0 : kB1 ];
1796   return fParams[axis==kX ? kB0 : kB1 ];
1797 }
1798
1799 //_____________________________________________________________
1800 void AliITSTPArrayFit::Print(Option_t *) const
1801 {
1802   // print results of the fit
1803   //
1804   const char kCxyz[] = "XYZ";
1805   if (!fkPoints) return;
1806   //
1807   printf("Track of %3d points in Bz=%+.1f |Fit ",fPntLast-fPntFirst+1,fBz);
1808   if ( IsFitDone() ) {
1809     if (IsHelix())
1810       printf("Helix: Chi2: %5.1f | %+.2e %+.2e %+.2e %+.2e %+.2e\n",
1811              fChi2NDF,fParams[kD0],fParams[kPhi0],fParams[kR0],fParams[kDZ],fParams[kDip]);
1812     else 
1813       printf("Line%c: Chi2: %5.1f | %+.2e %+.2e %+.2e %+.2e\n",
1814              kCxyz[fParAxis],fChi2NDF,fParams[kA0],fParams[kB0],fParams[kA1],fParams[kB1]);
1815   }
1816   else printf("N/A\n");
1817 }
1818
1819
1820
1821
1822 //____________________________________________________
1823 void AliITSTPArrayFit::BuildMaterialLUT(Int_t ntri) 
1824 {
1825   // Fill a look-up table with mean material a la AliITSTrackerMI
1826   //
1827   if (!AliGeomManager::GetGeometry()) AliFatal("Geometry is not loaded");
1828   //
1829   // detector layer to check: dX,dZ,Ymin,Ymax
1830   const double kLayr[9][4] =  {{0.  ,60. , 2.80,3.00},  // beam pipe
1831                                {1.28,7.07,-0.20,0.22},  // SPD1
1832                                {1.28,7.07,-0.20,0.22},  // SPD2
1833                                {0.  ,76.0, 10.4,11.8},  // Shield1
1834                                {7.02,7.53,-1.00,4.50},  // SDD1
1835                                {7.02,7.53,-1.00,4.50},  // SDD2
1836                                {0.  ,102., 29.0,30.0},  // Shield2
1837                                {7.50,4.20,-0.15,4.50},  // SSD1
1838                                {7.50,4.20,-0.15,4.50}}; // SSD2
1839   //
1840   //
1841   // build <dens*L> for detectors (track hitting the sensor in normal direction)
1842   double pg1[3],pg2[3],res[7];
1843   //
1844   int sID = 0;
1845   int actLrID = 0;
1846   for (int lr=0;lr<9;lr++) {
1847     //
1848     Bool_t active = kFALSE;
1849     const double* tpars = kLayr[lr];
1850     //
1851     if (IsZero(tpars[0])) { // passive layer
1852       active = kFALSE;
1853       AliInfo(Form("Probing passive layer (total layer #%d)",lr));
1854     }  
1855     else {
1856       active = kTRUE;
1857       sID += AliGeomManager::LayerSize(++actLrID);
1858       AliInfo(Form("Probing sensors of active layer #%d (total layers #%d)",actLrID,lr));
1859     }
1860     double shift = TMath::Abs(tpars[2]-tpars[3])*1E-4;
1861     double rhol = 0;
1862     for (int i=ntri;i--;) {
1863       //
1864       if (active) {
1865         int ssID = sID -1 - (int)(AliGeomManager::LayerSize(actLrID)*gRandom->Rndm());
1866         pg1[0] = pg2[0] = (gRandom->Rndm()-0.5)*tpars[0] + shift; // local X
1867         pg2[0] -= 2*shift;
1868         pg1[1] = tpars[2];
1869         pg2[1] = tpars[3];
1870         pg1[2] = pg2[2] = (gRandom->Rndm()-0.5)*tpars[1] + shift; // local Z
1871         pg2[2] -= 2*shift;
1872         AliITSgeomTGeo::LocalToGlobal(ssID,pg1,pg1);    
1873         AliITSgeomTGeo::LocalToGlobal(ssID,pg2,pg2);    
1874       }
1875       else {
1876         double ang = gRandom->Rndm()*TMath::Pi()*2;
1877         pg1[0] = tpars[2]*TMath::Cos(ang)+shift;
1878         pg2[0] = tpars[3]*TMath::Cos(ang)-shift;
1879         pg1[1] = tpars[2]*TMath::Sin(ang);
1880         pg2[1] = tpars[3]*TMath::Sin(ang);
1881         pg1[2] = pg2[2] = (gRandom->Rndm()-0.5)*tpars[1]+shift; // local Z
1882         pg2[2] -= 2*shift;
1883       }
1884
1885       //
1886       AliTracker::MeanMaterialBudget(pg1,pg2,res);
1887       rhol += res[0]*res[4];   // rho*L
1888     }
1889     fgRhoLITS[lr] = rhol/ntri;
1890     AliInfo(Form("Obtained <rho*L> = %e\n",fgRhoLITS[lr]));
1891   }
1892   //
1893   return;
1894 }
1895
1896 //____________________________________________________
1897 Double_t AliITSTPArrayFit::GetPCA2PlaneInfo(Double_t *xyz, Double_t *dir, Int_t axis, Double_t axval) const
1898 {
1899   // calculate the PCA to plane normal ti axis and crossing it at axval
1900   // fill the position and direction cosines at this point
1901   //
1902   double xyzp[3] = {0,0,0};                // create fake point
1903   xyzp[axis] = axval;
1904   double covI[6] = {1e-4,0,0,1e-4,0,1e-4}; // fake cov.matrix loose in all directions
1905   covI[4*axis - axis*(axis+1)/2] = 1e8;    // except axis
1906   //
1907   double t = GetPosition(xyz, xyzp, covI); // got pca
1908   //
1909   if (dir) GetDirCos(dir,t);
1910   return t;
1911 }
1912
1913 //____________________________________________________
1914 void AliITSTPArrayFit::GetT0Info(Double_t* xyz, Double_t *dir) const
1915 {
1916   // get direction cosines for t = 0;
1917   GetPosition(xyz, 0);
1918   if (dir) GetDirCos(dir,0);
1919 }
1920
1921 //____________________________________________________
1922 Bool_t AliITSTPArrayFit::CalcErrorMatrix()
1923 {
1924   // compute covariance matrix in lenear approximation of residuals on parameters
1925   static AliSymMatrix cv(5);
1926   static Double_t dres[5][3]; 
1927   cv.Reset();
1928   int npar = IsHelix() ? 5:4;
1929   //
1930   for (int ip=fPntFirst;ip<=fPntLast;ip++) {
1931     GetDResDParams(&dres[0][0],ip);
1932     Double_t* covI = GetCovI(ip);
1933     for (int i=npar;i--;) 
1934       for (int j=i+1;j--;) {
1935         double cvadd = dres[i][kX]*(dres[j][kX]*covI[ kXX ] + dres[j][kY]*covI[ kXY ] + dres[j][kZ]*covI[ kXZ ])
1936           +            dres[i][kY]*(dres[j][kX]*covI[ kXY ] + dres[j][kY]*covI[ kYY ] + dres[j][kZ]*covI[ kYZ ])
1937           +            dres[i][kZ]*(dres[j][kX]*covI[ kXZ ] + dres[j][kY]*covI[ kYZ ] + dres[j][kZ]*covI[ kZZ ]);
1938         if (covI[kScl]>0) cvadd *= covI[kScl];
1939         cv(i,j) += cvadd;
1940       }
1941   }
1942   cv.SetSizeUsed(npar);
1943   if (cv.InvertChol()) {
1944     //cv.Print("l");
1945     int cnt = 0;
1946     for (int i=npar;i--;) for (int j=i+1;j--;)fParamsCov[cnt++] = cv(i,j);
1947     return kTRUE;
1948   }
1949   //
1950   return kFALSE;
1951 }
1952
1953 //____________________________________________________
1954 Double_t AliITSTPArrayFit::CalcParPCA(Int_t ipnt) const
1955 {
1956   // get parameter for the point with least weighted distance to the point
1957   const double *xyz  = GetPoint(ipnt);
1958   const double *covI = GetCovI(ipnt);
1959   if (IsHelix()) return GetParPCAHelix(xyz,covI,covI[kScl]);
1960   else           return GetParPCALine(xyz,covI/*,covI[kScl]*/);
1961 }
1962
1963 //____________________________________________________
1964 Double_t AliITSTPArrayFit::GetPt() const 
1965 {
1966   return IsFieldON()&&IsHelix() ? TMath::Abs(fParams[kR0]*fBz*fgkCQConv) : -1;
1967 }
1968
1969 //____________________________________________________
1970 Double_t AliITSTPArrayFit::GetP() const 
1971 {
1972   if (!IsFieldON()) return -1;
1973   Double_t cdip = 1./TMath::Sqrt(1.+fParams[kDip]*fParams[kDip]);
1974   return TMath::Abs(fParams[kR0]*fgkCQConv*fBz/cdip); // momentum
1975 }
1976