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