ITS/AliITSTPArrayFit.cxx

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