Another histos for lumi
[u/mrichter/AliRoot.git] / ITS / AliITStrackV2.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////
19 //                Implementation of the ITS track class
20 //
21 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
22 //     dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
23 ///////////////////////////////////////////////////////////////////////////
24 #include <TMath.h>
25
26 #include "AliCluster.h"
27 #include "AliESDVertex.h"
28 #include "AliITSReconstructor.h"
29 #include "AliITStrackV2.h"
30 #include "AliTracker.h"
31 #include "AliLog.h"
32
33 const Int_t AliITStrackV2::fgkWARN = 5;
34
35 ClassImp(AliITStrackV2)
36
37
38 //____________________________________________________________________________
39 AliITStrackV2::AliITStrackV2() : AliKalmanTrack(),
40   fCheckInvariant(kTRUE),
41   fdEdx(0),
42   fESDtrack(0)
43 {
44   for(Int_t i=0; i<2*AliITSgeomTGeo::kNLayers; i++) {fIndex[i]=-1; fModule[i]=-1;}
45   for(Int_t i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=0;}
46   for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
47 }
48
49
50 //____________________________________________________________________________
51 AliITStrackV2::AliITStrackV2(AliESDtrack& t,Bool_t c):
52   AliKalmanTrack(),
53   fCheckInvariant(kTRUE),
54   fdEdx(t.GetITSsignal()),
55   fESDtrack(&t)
56 {
57   //------------------------------------------------------------------
58   // Conversion ESD track -> ITS track.
59   // If c==kTRUE, create the ITS track out of the constrained params.
60   //------------------------------------------------------------------
61   const AliExternalTrackParam *par=&t;
62   if (c) {
63     par=t.GetConstrainedParam();
64     if (!par) AliError("AliITStrackV2: conversion failed !\n");
65   }
66   Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
67
68   SetLabel(t.GetITSLabel());
69   SetMass(t.GetMass());
70   SetNumberOfClusters(t.GetITSclusters(fIndex));
71
72   if (t.GetStatus()&AliESDtrack::kTIME) {
73     StartTimeIntegral();
74     Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
75     SetIntegratedLength(t.GetIntegratedLength());
76   }
77
78   for(Int_t i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=0;}
79   for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
80 }
81
82 //____________________________________________________________________________
83 void AliITStrackV2::ResetClusters() {
84   //------------------------------------------------------------------
85   // Reset the array of attached clusters.
86   //------------------------------------------------------------------
87   for (Int_t i=0; i<2*AliITSgeomTGeo::kNLayers; i++) fIndex[i]=-1;
88   for (Int_t i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=0;}
89   SetChi2(0.); 
90   SetNumberOfClusters(0);
91
92
93 //____________________________________________________________________________
94 void AliITStrackV2::UpdateESDtrack(ULong_t flags) const {
95   // Update track params
96   fESDtrack->UpdateTrackParams(this,flags);
97   //
98   // set correctly the global label
99   if (fESDtrack->IsOn(AliESDtrack::kTPCin)) { 
100     // for global track the GetLabel should be negative if
101     // 1) GetTPCLabel<0
102     // 2) this->GetLabel()<0
103     // 3) GetTPCLabel() != this->GetLabel()
104     int label = fESDtrack->GetTPCLabel();
105     int itsLabel = GetLabel();
106     if (label<0 || itsLabel<0 || itsLabel!=label) label = -TMath::Abs(label);
107     fESDtrack->SetLabel(label);
108   }
109   //
110   // copy the module indices
111   Int_t i;
112   for(i=0;i<2*AliITSgeomTGeo::kNLayers;i++) {
113     //   printf("     %d\n",GetModuleIndex(i));
114     fESDtrack->SetITSModuleIndex(i,GetModuleIndex(i));
115   }
116   // copy the map of shared clusters
117   if(flags==AliESDtrack::kITSin) {
118     UChar_t itsSharedMap=0;
119     for(i=0;i<AliITSgeomTGeo::kNLayers;i++) {
120       if(fSharedWeight[i]>0) SETBIT(itsSharedMap,i);
121       
122     }
123     fESDtrack->SetITSSharedMap(itsSharedMap);
124   }
125
126   // copy the 4 dedx samples
127   Double_t sdedx[4]={0.,0.,0.,0.};
128   for(i=0; i<4; i++) sdedx[i]=fdEdxSample[i];
129   fESDtrack->SetITSdEdxSamples(sdedx);
130 }
131
132 //____________________________________________________________________________
133 AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : 
134   AliKalmanTrack(t),
135   fCheckInvariant(t.fCheckInvariant),
136   fdEdx(t.fdEdx),
137   fESDtrack(t.fESDtrack) 
138 {
139   //------------------------------------------------------------------
140   //Copy constructor
141   //------------------------------------------------------------------
142   Int_t i;
143   for (i=0; i<4; i++) fdEdxSample[i]=t.fdEdxSample[i];
144   for (i=0; i<2*AliITSgeomTGeo::GetNLayers(); i++) {
145     fIndex[i]=t.fIndex[i];
146     fModule[i]=t.fModule[i];
147   }
148   for (i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=t.fSharedWeight[i];}
149 }
150
151 //_____________________________________________________________________________
152 Int_t AliITStrackV2::Compare(const TObject *o) const {
153   //-----------------------------------------------------------------
154   // This function compares tracks according to the their curvature
155   //-----------------------------------------------------------------
156   AliITStrackV2 *t=(AliITStrackV2*)o;
157   //Double_t co=OneOverPt();
158   //Double_t c =OneOverPt();
159   Double_t co=t->GetSigmaY2()*t->GetSigmaZ2();
160   Double_t c =GetSigmaY2()*GetSigmaZ2();
161   if (c>co) return 1;
162   else if (c<co) return -1;
163   return 0;
164 }
165
166 //____________________________________________________________________________
167 Bool_t 
168 AliITStrackV2::PropagateToVertex(const AliESDVertex *v,Double_t d,Double_t x0) 
169 {
170   //------------------------------------------------------------------
171   //This function propagates a track to the minimal distance from the origin
172   //------------------------------------------------------------------  
173   Double_t bz=GetBz();
174   if (PropagateToDCA(v,bz,kVeryBig)) {
175     Double_t xOverX0,xTimesRho; 
176     xOverX0 = d; xTimesRho = d*x0;
177     if (CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kTRUE;
178   }
179   return kFALSE;
180 }
181
182 //____________________________________________________________________________
183 Bool_t AliITStrackV2::
184 GetGlobalXYZat(Double_t xloc, Double_t &x, Double_t &y, Double_t &z) const {
185   //------------------------------------------------------------------
186   //This function returns a track position in the global system
187   //------------------------------------------------------------------
188   Double_t r[3];
189   Bool_t rc=GetXYZAt(xloc, GetBz(), r);
190   x=r[0]; y=r[1]; z=r[2]; 
191   return rc;
192 }
193
194 //_____________________________________________________________________________
195 Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const {
196   //-----------------------------------------------------------------
197   // This function calculates a predicted chi2 increment.
198   //-----------------------------------------------------------------
199   Double_t p[2]={c->GetY(), c->GetZ()};
200   Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()};
201   return AliExternalTrackParam::GetPredictedChi2(p,cov);
202 }
203
204 //____________________________________________________________________________
205 Bool_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) {
206   //------------------------------------------------------------------
207   //This function propagates a track
208   //------------------------------------------------------------------
209
210   Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ();
211   
212   //Double_t bz=GetBz();
213   //if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE;
214   Double_t b[3]; GetBxByBz(b);
215   if (!AliExternalTrackParam::PropagateToBxByBz(xk,b)) return kFALSE;
216   Double_t xOverX0,xTimesRho; 
217   xOverX0 = d; xTimesRho = d*x0;
218   if (!CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kFALSE;
219
220   Double_t x=GetX(), y=GetY(), z=GetZ();
221   if (IsStartedTimeIntegral() && x>oldX) {
222     Double_t l2 = (x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ);
223     AddTimeStep(TMath::Sqrt(l2));
224   }
225
226   return kTRUE;
227 }
228
229 //____________________________________________________________________________
230 Bool_t AliITStrackV2::PropagateToTGeo(Double_t xToGo, Int_t nstep, Double_t &xOverX0, Double_t &xTimesRho, Bool_t addTime) {
231   //-------------------------------------------------------------------
232   //  Propagates the track to a reference plane x=xToGo in n steps.
233   //  These n steps are only used to take into account the curvature.
234   //  The material is calculated with TGeo. (L.Gaudichet)
235   //-------------------------------------------------------------------
236   
237   Double_t startx = GetX(), starty = GetY(), startz = GetZ();
238   Double_t sign = (startx<xToGo) ? -1.:1.;
239   Double_t step = (xToGo-startx)/TMath::Abs(nstep);
240
241   Double_t start[3], end[3], mparam[7];
242   //Double_t bz = GetBz();
243   Double_t b[3]; GetBxByBz(b);
244   Double_t bz = b[2];
245
246   Double_t x = startx;
247   
248   for (Int_t i=0; i<nstep; i++) {
249     
250     GetXYZ(start);   //starting global position
251     x += step;
252     if (!GetXYZAt(x, bz, end)) return kFALSE;
253     //if (!AliExternalTrackParam::PropagateTo(x, bz)) return kFALSE;
254     if (!AliExternalTrackParam::PropagateToBxByBz(x, b)) return kFALSE;
255     AliTracker::MeanMaterialBudget(start, end, mparam);
256     xTimesRho = sign*mparam[4]*mparam[0];
257     xOverX0   = mparam[1];
258     if (mparam[1]<900000) {
259       if (!AliExternalTrackParam::CorrectForMeanMaterial(xOverX0,
260                            xTimesRho,GetMass())) return kFALSE;
261     } else { // this happens when MeanMaterialBudget cannot cross a boundary
262       return kFALSE;
263     }
264   }
265
266   if (addTime && IsStartedTimeIntegral() && GetX()>startx) {
267     Double_t l2 = ( (GetX()-startx)*(GetX()-startx) +
268                     (GetY()-starty)*(GetY()-starty) +
269                     (GetZ()-startz)*(GetZ()-startz) );
270     AddTimeStep(TMath::Sqrt(l2));
271   }
272
273   return kTRUE;
274 }
275
276 //____________________________________________________________________________
277 Bool_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, Int_t index) 
278 {
279   //------------------------------------------------------------------
280   //This function updates track parameters
281   //------------------------------------------------------------------
282   Double_t p[2]={c->GetY(), c->GetZ()};
283   Double_t cov[3]={c->GetSigmaY2(), c->GetSigmaYZ(), c->GetSigmaZ2()};
284
285   if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
286
287   Int_t n=GetNumberOfClusters();
288   if (!Invariant()) {
289     if (n>fgkWARN) AliDebug(1,"Wrong invariant !");
290      return kFALSE;
291   }
292
293   if (chi2<0) return kTRUE;
294
295   // fill residuals for ITS+TPC tracks 
296   if (fESDtrack) {
297     if (fESDtrack->GetStatus()&AliESDtrack::kTPCin) {
298       AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
299     }
300   }
301
302   fIndex[n]=index;
303   SetNumberOfClusters(n+1);
304   SetChi2(GetChi2()+chi2);
305
306   return kTRUE;
307 }
308
309 Bool_t AliITStrackV2::Invariant() const {
310   //------------------------------------------------------------------
311   // This function is for debugging purpose only
312   //------------------------------------------------------------------
313   if(!fCheckInvariant) return kTRUE;
314
315   Int_t n=GetNumberOfClusters();
316   static Float_t bz = GetBz();
317   // take into account the misalignment error
318   Float_t maxMisalErrY2=0,maxMisalErrZ2=0;
319   //RS
320   const AliITSRecoParam* recopar = AliITSReconstructor::GetRecoParam();
321   if (!recopar) recopar = AliITSRecoParam::GetHighFluxParam();
322
323   for (Int_t lay=0; lay<AliITSgeomTGeo::kNLayers; lay++) {
324     maxMisalErrY2 = TMath::Max(maxMisalErrY2,recopar->GetClusterMisalErrorY(lay,bz));
325     maxMisalErrZ2 = TMath::Max(maxMisalErrZ2,recopar->GetClusterMisalErrorZ(lay,bz));
326   }
327   maxMisalErrY2 *= maxMisalErrY2;
328   maxMisalErrZ2 *= maxMisalErrZ2;
329   // this is because when we reset before refitting, we multiply the
330   // matrix by 10
331   maxMisalErrY2 *= 10.; 
332   maxMisalErrZ2 *= 10.;
333
334   Double_t sP2=GetParameter()[2];
335   if (TMath::Abs(sP2) >= kAlmost1){
336     if (n>fgkWARN) AliDebug(1,Form("fP2=%f\n",sP2));
337      return kFALSE;
338   }
339   Double_t sC00=GetCovariance()[0];
340   if (sC00<=0 || sC00>(9.+maxMisalErrY2)) {
341     if (n>fgkWARN) AliDebug(1,Form("fC00=%f\n",sC00)); 
342      return kFALSE;
343   }
344   Double_t sC11=GetCovariance()[2];
345   if (sC11<=0 || sC11>(9.+maxMisalErrZ2)) {
346     if (n>fgkWARN) AliDebug(1,Form("fC11=%f\n",sC11)); 
347      return kFALSE;
348   }
349   Double_t sC22=GetCovariance()[5];
350   if (sC22<=0 || sC22>1.) {
351     if (n>fgkWARN) AliDebug(1,Form("fC22=%f\n",sC22)); 
352      return kFALSE;
353   }
354   Double_t sC33=GetCovariance()[9];
355   if (sC33<=0 || sC33>1.) {
356     if (n>fgkWARN) AliDebug(1,Form("fC33=%f\n",sC33)); 
357      return kFALSE;
358   }
359   Double_t sC44=GetCovariance()[14];
360   if (sC44<=0 /*|| sC44>6e-5*/) {
361     if (n>fgkWARN) AliDebug(1,Form("fC44=%f\n",sC44));
362      return kFALSE;
363   }
364
365   return kTRUE;
366 }
367
368 //____________________________________________________________________________
369 Bool_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
370   //------------------------------------------------------------------
371   //This function propagates a track
372   //------------------------------------------------------------------
373   //Double_t bz=GetBz();
374   //if (!AliExternalTrackParam::Propagate(alp,xk,bz)) return kFALSE;
375   Double_t b[3]; GetBxByBz(b);
376   if (!AliExternalTrackParam::PropagateBxByBz(alp,xk,b)) return kFALSE;
377
378   if (!Invariant()) {
379     Int_t n=GetNumberOfClusters();
380     if (n>fgkWARN) AliDebug(1,"Wrong invariant !");
381     return kFALSE;
382   }
383
384   return kTRUE;
385 }
386
387 Bool_t AliITStrackV2::MeanBudgetToPrimVertex(Double_t xyz[3], Double_t step, Double_t &d) const {
388
389   //-------------------------------------------------------------------
390   //  Get the mean material budget between the actual point and the
391   //  primary vertex. (L.Gaudichet)
392   //-------------------------------------------------------------------
393
394   Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
395   Double_t vertexX = xyz[0]*cs + xyz[1]*sn;
396
397   Int_t nstep = Int_t((GetX()-vertexX)/step);
398   if (nstep<1) nstep = 1;
399   step = (GetX()-vertexX)/nstep;
400
401   //  Double_t mparam[7], densMean=0, radLength=0, length=0;
402   Double_t mparam[7];
403   Double_t p1[3], p2[3], x = GetX(), bz = GetBz();
404   GetXYZ(p1);
405
406   d=0.;
407
408   for (Int_t i=0; i<nstep; i++) {
409     x  += step;
410     if (!GetXYZAt(x, bz, p2)) return kFALSE;
411     AliTracker::MeanMaterialBudget(p1, p2, mparam);
412     if (mparam[1]>900000) return kFALSE;
413     d  += mparam[1];
414
415     p1[0] = p2[0];
416     p1[1] = p2[1];
417     p1[2] = p2[2];
418   }
419
420   return kTRUE;
421 }
422
423 Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) {
424   //------------------------------------------------------------------
425   //This function improves angular track parameters
426   //------------------------------------------------------------------
427   //Store the initail track parameters
428
429   Double_t x = GetX();
430   Double_t alpha = GetAlpha();
431   Double_t par[] = {GetY(),GetZ(),GetSnp(),GetTgl(),GetSigned1Pt()};
432   Double_t cov[] = {
433     GetSigmaY2(),
434     GetSigmaZY(),
435     GetSigmaZ2(),
436     GetSigmaSnpY(),
437     GetSigmaSnpZ(),
438     GetSigmaSnp2(),
439     GetSigmaTglY(),
440     GetSigmaTglZ(),
441     GetSigmaTglSnp(),
442     GetSigmaTgl2(),
443     GetSigma1PtY(),
444     GetSigma1PtZ(),
445     GetSigma1PtSnp(),
446     GetSigma1PtTgl(),
447     GetSigma1Pt2()
448   }; 
449
450
451   Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
452   Double_t xv = xyz[0]*cs + xyz[1]*sn; // vertex
453   Double_t yv =-xyz[0]*sn + xyz[1]*cs; // in the
454   Double_t zv = xyz[2];                // local frame
455
456   Double_t dx = x - xv, dy = par[0] - yv, dz = par[1] - zv;
457   Double_t r2=dx*dx + dy*dy;
458   Double_t p2=(1.+ GetTgl()*GetTgl())/(GetSigned1Pt()*GetSigned1Pt());
459   Double_t beta2=p2/(p2 + GetMass()*GetMass());
460   x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
461   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
462   //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
463
464   Double_t bz=GetBz();
465   Double_t cnv=bz*kB2C;
466   Double_t curv=GetC(bz);
467   {
468     Double_t dummy = 4/r2 - curv*curv;
469     if (dummy < 0) return kFALSE;
470     Double_t parp = 0.5*(curv*dx + dy*TMath::Sqrt(dummy));
471     Double_t sigma2p = theta2*(1.-GetSnp())*(1.+GetSnp())*(1. + GetTgl()*GetTgl());
472     Double_t ovSqr2 = 1./TMath::Sqrt(r2);
473     Double_t tfact = ovSqr2*(1.-dy*ovSqr2)*(1.+dy*ovSqr2);
474     sigma2p += cov[0]*tfact*tfact;
475     sigma2p += ers[1]*ers[1]/r2;
476     sigma2p += 0.25*cov[14]*cnv*cnv*dx*dx;
477     Double_t eps2p=sigma2p/(cov[5] + sigma2p);
478     par[0] += cov[3]/(cov[5] + sigma2p)*(parp - GetSnp());
479     par[2]  = eps2p*GetSnp() + (1 - eps2p)*parp;
480     cov[5] *= eps2p;
481     cov[3] *= eps2p;
482   }
483   {
484     Double_t parl=0.5*curv*dz/TMath::ASin(0.5*curv*TMath::Sqrt(r2));
485     Double_t sigma2l=theta2;
486     sigma2l += cov[2]/r2 + cov[0]*dy*dy*dz*dz/(r2*r2*r2);
487     sigma2l += ers[2]*ers[2]/r2;
488     Double_t eps2l = sigma2l/(cov[9] + sigma2l);
489     par[1] += cov[7 ]/(cov[9] + sigma2l)*(parl - par[3]);
490     par[4] += cov[13]/(cov[9] + sigma2l)*(parl - par[3]);
491     par[3]  = eps2l*par[3] + (1-eps2l)*parl;
492     cov[9] *= eps2l; 
493     cov[13]*= eps2l; 
494     cov[7] *= eps2l; 
495   }
496
497   Set(x,alpha,par,cov);
498
499   if (!Invariant()) return kFALSE;
500
501   return kTRUE;
502 }
503
504 void AliITStrackV2::CookdEdx(Double_t /*low*/, Double_t /*up*/) {
505   //-----------------------------------------------------------------
506   // This function calculates dE/dX within the "low" and "up" cuts.
507   // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch 
508   // Updated: F. Prino 8-June-2009
509   //-----------------------------------------------------------------
510   // The cluster order is: SDD-1, SDD-2, SSD-1, SSD-2
511
512   Int_t nc=0;
513   Float_t dedx[4];
514   for (Int_t il=0; il<4; il++) { // count good (>0) dE/dx values
515     if(fdEdxSample[il]>0.){
516       dedx[nc]= fdEdxSample[il];
517       nc++;
518     }
519   }
520   if(nc<1){
521     SetdEdx(0.);
522     return;
523   }
524
525   Int_t swap; // sort in ascending order
526   do {
527     swap=0;
528     for (Int_t i=0; i<nc-1; i++) {
529       if (dedx[i]<=dedx[i+1]) continue;
530       Float_t tmp=dedx[i];
531       dedx[i]=dedx[i+1]; 
532       dedx[i+1]=tmp;
533       swap++;
534     }
535   } while (swap);
536
537
538   Double_t sumamp=0,sumweight=0;
539   Double_t weight[4]={1.,1.,0.,0.};
540   if(nc==3) weight[1]=0.5;
541   else if(nc<3) weight[1]=0.;
542   for (Int_t i=0; i<nc; i++) {
543     sumamp+= dedx[i]*weight[i];
544     sumweight+=weight[i];
545   }
546   SetdEdx(sumamp/sumweight);
547 }
548
549 //____________________________________________________________________________
550 Bool_t AliITStrackV2::
551 GetPhiZat(Double_t r, Double_t &phi, Double_t &z) const {
552   //------------------------------------------------------------------
553   // This function returns the global cylindrical (phi,z) of the track 
554   // position estimated at the radius r. 
555   // The track curvature is neglected.
556   //------------------------------------------------------------------
557   Double_t d=GetD(0.,0.);
558   if (TMath::Abs(d) > r) {
559     if (r>1e-1) return kFALSE;
560     r = TMath::Abs(d);
561   }
562
563   Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
564   if (TMath::Abs(d) > rcurr) return kFALSE;
565   Double_t globXYZcurr[3]; GetXYZ(globXYZcurr); 
566   Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]);
567
568   if (GetX()>=0.) {
569     phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
570   } else {
571     phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi();
572   }
573
574   // return a phi in [0,2pi[ 
575   if (phi<0.) phi+=2.*TMath::Pi();
576   else if (phi>=2.*TMath::Pi()) phi-=2.*TMath::Pi();
577   z=GetZ()+GetTgl()*(TMath::Sqrt((r-d)*(r+d))-TMath::Sqrt((rcurr-d)*(rcurr+d)));
578   return kTRUE;
579 }
580 //____________________________________________________________________________
581 Bool_t AliITStrackV2::
582 GetLocalXat(Double_t r,Double_t &xloc) const {
583   //------------------------------------------------------------------
584   // This function returns the local x of the track 
585   // position estimated at the radius r. 
586   // The track curvature is neglected.
587   //------------------------------------------------------------------
588   Double_t d=GetD(0.,0.);
589   if (TMath::Abs(d) > r) { 
590     if (r>1e-1) return kFALSE; 
591     r = TMath::Abs(d); 
592   } 
593
594   Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
595   Double_t globXYZcurr[3]; GetXYZ(globXYZcurr); 
596   Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]);
597   Double_t phi;
598   if (GetX()>=0.) {
599     phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
600   } else {
601     phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi();
602   }
603
604   xloc=r*(TMath::Cos(phi)*TMath::Cos(GetAlpha())
605          +TMath::Sin(phi)*TMath::Sin(GetAlpha())); 
606
607   return kTRUE;
608 }
609
610 //____________________________________________________________________________
611 Bool_t AliITStrackV2::ImproveKalman(Double_t xyz[3],Double_t ers[3], const Double_t* xlMS, const Double_t* x2X0MS, Int_t nMS)
612 {
613   // Substitute the state of the track (p_{k|k},C_{k|k}) at the k-th measumerent by its
614   // smoothed value from the k-th measurement + measurement at the vertex.
615   // Account for the MS on nMS layers at x-postions xlMS with x/x0 = x2X0MS
616   // p_{k|kv} = p_{k|k} + C_{k|k}*D^Tr_{k+1} B^{-1}_{k+1} ( vtx - D_{k+1}*p_{k|k})
617   // C_{k|kv} = C_{k|k}*( I - D^Tr_{k+1} B^{-1}_{k+1} D_{k+1} C^Tr_{k|k})
618   // 
619   // where D_{k} = H_{k} F_{k} with H being the matrix converting the tracks parameters
620   // to measurements m_{k} = H_{k} p_{k} and F_{k} the matrix propagating the track between the
621   // the point k-1 and k:  p_{k|k-1} = F_{k} p_{k-1|k-1}
622   //
623   // B_{k+1} = V_{k+1} + H_{k+1} C_{k+1|k} H^Tr_{k+1} with V_{k+1} being the error of the measurment
624   // at point k+1 (i.e. vertex), and C_{k+1|k} - error matrix extrapolated from k-th measurement to
625   // k+1 (vtx) and accounting for the MS inbetween
626   //
627   // H = {{1,0,0,0,0},{0,1,0,0,0}}
628   //
629   double covc[15], *cori = (double*) GetCovariance(),par[5] = {GetY(),GetZ(),GetSnp(),GetTgl(),GetSigned1Pt()},
630     &c00=cori[0],
631     &c01=cori[1],&c11=cori[2],
632     &c02=cori[3],&c12=cori[4],&c22=cori[5],
633     &c03=cori[6],&c13=cori[7],&c23=cori[8],&c33=cori[9],
634     &c04=cori[10],&c14=cori[11],&c24=cori[12],&c34=cori[13],&c44=cori[14],
635     // for smoothed cov matrix 
636     &cov00=covc[0],
637     &cov01=covc[1],&cov11=covc[2],
638     &cov02=covc[3],&cov12=covc[4],&cov22=covc[5],
639     &cov03=covc[6],&cov13=covc[7],&cov23=covc[8],&cov33=covc[9],
640     &cov04=covc[10],&cov14=covc[11],&cov24=covc[12],&cov34=covc[13],&cov44=covc[14];
641   //
642   double x = GetX(), alpha = GetAlpha();
643   // vertex in the track frame
644   double cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
645   double xv = xyz[0]*cs + xyz[1]*sn, yv =-xyz[0]*sn + xyz[1]*cs, zv = xyz[2];
646   double dx = xv - GetX();
647   if (TMath::Abs(dx)<=kAlmost0)  return kTRUE;
648   //
649   double cnv=GetBz()*kB2C, x2r=cnv*par[4]*dx, f1=par[2], f2=f1+x2r;
650   if (TMath::Abs(f1) >= kAlmost1 || TMath::Abs(f2) >= kAlmost1) {
651     AliInfo(Form("Fail: %+e %+e",f1,f2));
652     return kFALSE;
653   }
654   double r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2)), dx2r=dx/(r1+r2);
655   // elements of matrix F_{k+1} (1s on diagonal)
656   double f02 = 2*dx2r, f04 = cnv*dx*dx2r, f13/*, f24 = cnv*dx*/;
657   if (TMath::Abs(x2r)<0.05) f13 = dx*r2+f2*(f1+f2)*dx2r; // see AliExternalTrackParam::PropagateTo
658   else {
659     double dy2dx = (f1+f2)/(r1+r2);
660     f13 = 2*TMath::ASin(0.5*TMath::Sqrt(1+dy2dx*dy2dx)*x2r)/(cnv*par[4]);
661   }  
662   // elements of matrix D_{k+1} = H_{k+1} * F_{k+1}
663   // double d00 = 1., d11 = 1.;
664   double &d02 = f02,  &d04 = f04, &d13 = f13;
665   //
666   // elements of matrix DC = D_{k+1}*C_{kk}^T
667   double dc00 = c00+c02*d02+c04*d04,   dc10 = c01+c03*d13;
668   double dc01 = c01+c12*d02+c14*d04,   dc11 = c11+c13*d13;
669   double dc02 = c02+c22*d02+c24*d04,   dc12 = c12+c23*d13;
670   double dc03 = c03+c23*d02+c34*d04,   dc13 = c13+c33*d13;
671   double dc04 = c04+c24*d02+c44*d04,   dc14 = c14+c34*d13;
672   //
673   // difference between the vertex and the the track extrapolated to vertex
674   yv -= par[0] + par[2]*d02 + par[4]*d04;
675   zv -= par[1] + par[3]*d13;
676   //
677   // y,z part of the cov.matrix extrapolated to vtx (w/o MS contribution)
678   // C_{k+1,k} = H F_{k+1} C_{k,k} F^Tr_{k+1} H^Tr = D C D^Tr
679   double cv00 = dc00+dc02*d02+dc04*d04, cv01 = dc01+dc03*d13, cv11 = dc11+dc13*d13;
680   //
681   // add MS contribution layer by layer
682   double xCurr = x;
683   double p2Curr = par[2];
684   //
685   // precalculated factors of MS contribution matrix:
686   double ms22t = (1. + par[3]*par[3]);
687   double ms33t = ms22t*ms22t;
688   double p34 = par[3]*par[4];
689   double ms34t = p34*ms22t;
690   double ms44t = p34*p34;
691   //
692   double p2=(1.+ par[3]*par[3])/(par[4]*par[4]);
693   double beta2 = p2/(p2+GetMass()*GetMass());
694   double theta2t = 14.1*14.1/(beta2*p2*1e6) * (1. + par[3]*par[3]);
695   //
696   // account for the MS in the layers between the last measurement and the vertex
697   for (int il=0;il<nMS;il++) {
698     double dfx = xlMS[il] - xCurr;
699     xCurr = xlMS[il];
700     p2Curr += dfx*cnv*par[4];   // p2 at the scattering layer
701     double dxL=xv - xCurr;    // distance from scatering layer to vtx
702     double x2rL=cnv*par[4]*dxL, f1L=p2Curr, f2L=f1L+x2rL;
703     if (TMath::Abs(f1L) >= kAlmost1 || TMath::Abs(f2L) >= kAlmost1) {
704       AliInfo(Form("FailMS at step %d of %d: dfx:%e dxL:%e %e %e",il,nMS,dfx,dxL,f1L,f2L));
705       return kFALSE;
706     }
707     double r1L=TMath::Sqrt((1.-f1L)*(1.+f1L)), r2L=TMath::Sqrt((1.-f2L)*(1.+f2L)), dx2rL=dxL/(r1L+r2L);
708     // elements of matrix for propagation from scatering layer to vertex
709     double f02L = 2*dx2rL, f04L = cnv*dxL*dx2rL, f13L/*, f24L = cnv*dxL*/;
710     if (TMath::Abs(x2rL)<0.05) f13L = dxL*r2L+f2L*(f1L+f2L)*dx2rL; // see AliExternalTrackParam::PropagateTo
711     else {
712       double dy2dxL = (f1L+f2L)/(r1L+r2L);
713       f13L = 2*TMath::ASin(0.5*TMath::Sqrt(1+dy2dxL*dy2dxL)*x2rL)/(cnv*par[4]);
714     }
715     // MS contribution matrix:
716     double theta2 = theta2t*TMath::Abs(x2X0MS[il]);
717     double ms22 = theta2*(1.-p2Curr)*(1.+p2Curr)*ms22t;
718     double ms33 = theta2*ms33t;
719     double ms34 = theta2*ms34t;
720     double ms44 = theta2*ms44t;
721     //
722     // add  H F MS F^Tr H^Tr to cv
723     cv00 += f02L*f02L*ms22 + f04L*f04L*ms44;
724     cv01 += f04L*f13L*ms34;
725     cv11 += f13L*f13L*ms33;
726   }
727   //
728   // inverse of matrix B
729   double b11 = ers[1]*ers[1] + cv00;
730   double b00 = ers[2]*ers[2] + cv11;
731   double det = b11*b00 - cv01*cv01;
732   if (TMath::Abs(det)<kAlmost0) {
733     AliInfo(Form("Fail on det %e: %e %e %e",det,cv00,cv11,cv01));
734     return kFALSE;
735   }
736   det = 1./det;
737   b00 *= det; b11 *= det; 
738   double b01 = -cv01*det;
739   //
740   // elements of matrix DC^Tr * B^-1
741   double dcb00 = b00*dc00+b01*dc10, dcb01 = b01*dc00+b11*dc10;
742   double dcb10 = b00*dc01+b01*dc11, dcb11 = b01*dc01+b11*dc11;
743   double dcb20 = b00*dc02+b01*dc12, dcb21 = b01*dc02+b11*dc12;
744   double dcb30 = b00*dc03+b01*dc13, dcb31 = b01*dc03+b11*dc13;
745   double dcb40 = b00*dc04+b01*dc14, dcb41 = b01*dc04+b11*dc14;
746   //
747   // p_{k|k+1} = p_{k|k} +  C_{k|k}*D^Tr_{k+1} B^{-1}_{k+1} ( vtx - D_{k+1}*p_{k|k})
748   par[0] += dcb00*yv + dcb01*zv;
749   par[1] += dcb10*yv + dcb11*zv;
750   par[2] += dcb20*yv + dcb21*zv;
751   par[3] += dcb30*yv + dcb31*zv;
752   par[4] += dcb40*yv + dcb41*zv;
753   //
754   // C_{k|kv} = C_{k|k} - C_{k|k} D^Tr_{k+1} B^{-1}_{k+1} D_{k+1} C^Tr_{k|k})
755   //
756   cov00 = c00 - (dc00*dcb00 + dc10*dcb01);
757   cov01 = c01 - (dc01*dcb00 + dc11*dcb01);
758   cov02 = c02 - (dc02*dcb00 + dc12*dcb01);
759   cov03 = c03 - (dc03*dcb00 + dc13*dcb01);
760   cov04 = c04 - (dc04*dcb00 + dc14*dcb01);
761   //
762   cov11 = c11 - (dc01*dcb10 + dc11*dcb11);
763   cov12 = c12 - (dc02*dcb10 + dc12*dcb11);
764   cov13 = c13 - (dc03*dcb10 + dc13*dcb11);
765   cov14 = c14 - (dc04*dcb10 + dc14*dcb11);
766   //
767   cov22 = c22 - (dc02*dcb20 + dc12*dcb21);
768   cov23 = c23 - (dc03*dcb20 + dc13*dcb21);
769   cov24 = c24 - (dc04*dcb20 + dc14*dcb21);
770   //
771   cov33 = c33 - (dc03*dcb30 + dc13*dcb31);
772   cov34 = c34 - (dc04*dcb30 + dc14*dcb31);
773   //
774   cov44 = c44 - (dc04*dcb40 + dc14*dcb41);
775   //
776   Set(x,alpha,par,covc);
777   if (!Invariant()) {
778     AliInfo(Form("Fail on Invariant, X=%e",GetX()));
779     return kFALSE;
780   }
781   return kTRUE;
782   //
783 }