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