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