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