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