ITS/AliITSTPArrayFit.cxx

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