Checking if the field map is set and propagating the tracks to the found vertex
[u/mrichter/AliRoot.git] / TRD / AliTRDtrack.cxx
CommitLineData
46d29e70 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
0fa7dfa7 16/* $Id$ */
46d29e70 17
a2cb5b3d 18#include <Riostream.h>
53c17fbf 19#include <TMath.h>
e1e6896f 20#include <TVector2.h>
46d29e70 21
6c94f330 22#include "AliTracker.h"
53c17fbf 23#include "AliESDtrack.h"
46d29e70 24#include "AliTRDgeometry.h"
25#include "AliTRDcluster.h"
26#include "AliTRDtrack.h"
53c17fbf 27#include "AliTRDtracklet.h"
b3a5a838 28
46d29e70 29ClassImp(AliTRDtrack)
30
53c17fbf 31///////////////////////////////////////////////////////////////////////////////
32// //
33// Represents a reconstructed TRD track //
34// Local TRD Kalman track //
35// //
36///////////////////////////////////////////////////////////////////////////////
7ad19338 37
6c94f330 38AliTRDtrack::AliTRDtrack():
39 AliKalmanTrack(),
40 fSeedLab(-1),
41 fdEdx(0),
42 fdEdxT(0),
43 fDE(0),
44 fStopped(kFALSE),
45 fLhElectron(0),
46 fNWrong(0),
47 fNRotate(0),
48 fNCross(0),
49 fNExpected(0),
50 fNLast(0),
51 fNExpectedLast(0),
52 fNdedx(0),
53 fChi2Last(1e10),
54 fBackupTrack(0x0)
f146da82 55{
6c94f330 56 for (Int_t i=0; i<kNplane; i++) {
57 for (Int_t j=0; j<kNslice; j++) {
6d45eaef 58 fdEdxPlane[i][j] = 0;
59 }
f146da82 60 fTimBinPlane[i] = -1;
61 }
6c94f330 62 for (UInt_t i=0; i<kMAXCLUSTERSPERTRACK; i++) {
63 fIndex[i] = 0;
64 fIndexBackup[i] = 0;
65 fdQdl[i] = 0;
f146da82 66 }
6c94f330 67 for (Int_t i=0; i<3; i++) fBudget[i] = 0;
f146da82 68}
6d45eaef 69
46d29e70 70//_____________________________________________________________________________
6c94f330 71AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, Int_t index,
72 const Double_t p[5], const Double_t cov[15],
73 Double_t x, Double_t alpha) :
74 AliKalmanTrack(),
75 fSeedLab(-1),
76 fdEdx(0),
77 fdEdxT(0),
78 fDE(0),
79 fStopped(kFALSE),
80 fLhElectron(0),
81 fNWrong(0),
82 fNRotate(0),
83 fNCross(0),
84 fNExpected(0),
85 fNLast(0),
86 fNExpectedLast(0),
87 fNdedx(0),
88 fChi2Last(1e10),
89 fBackupTrack(0x0)
15ed8ba1 90{
6c94f330 91 //-----------------------------------------------------------------
92 // This is the main track constructor.
93 //-----------------------------------------------------------------
94 Double_t cnv=1./(GetBz()*kB2C);
95
96 Double_t pp[5]={
97 p[0],
98 p[1],
99 x*p[4] - p[2],
100 p[3],
101 p[4]*cnv
102 };
103
104 Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[5];
105 Double_t c32 = x*cov[13] - cov[8];
106 Double_t c20 = x*cov[10] - cov[3],
107 c21 = x*cov[11] - cov[4], c42 = x*cov[14] - cov[12];
108
109 Double_t cc[15]={
110 cov[0 ],
111 cov[1 ], cov[2 ],
112 c20, c21, c22,
113 cov[6 ], cov[7 ], c32, cov[9 ],
114 cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv
115 };
116
117 Set(x,alpha,pp,cc);
d49e463b 118
5443e65e 119 SetNumberOfClusters(1);
6c94f330 120 fIndex[0]=index;
5443e65e 121
6c94f330 122 for (Int_t i=0;i<kNplane;i++){
123 for (Int_t j=0; j<kNslice; j++) {
6d45eaef 124 fdEdxPlane[i][j] = 0;
125 }
126 fTimBinPlane[i] = -1;
eab5961e 127 }
a2b90f83 128
5443e65e 129 Double_t q = TMath::Abs(c->GetQ());
6c94f330 130 Double_t s = GetSnp(), t=GetTgl();
131 if(s*s < 1) q *= TMath::Sqrt((1-s*s)/(1+t*t));
53c17fbf 132
6c94f330 133 fdQdl[0] = q;
134 // initialisation [SR, GSI 18.02.2003] (i startd for 1)
135 for(UInt_t i=1; i<kMAXCLUSTERSPERTRACK; i++) {
136 fdQdl[i] = 0;
137 fIndex[i] = 0;
138 fIndexBackup[i] = 0; //backup indexes MI
15ed8ba1 139 }
6c94f330 140 for (Int_t i=0;i<3;i++) fBudget[i]=0;
46d29e70 141}
142
143//_____________________________________________________________________________
6c94f330 144AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) :
145 AliKalmanTrack(t),
146 fSeedLab(t.GetSeedLabel()),
147 fdEdx(t.fdEdx),
148 fdEdxT(t.fdEdx),
149 fDE(t.fDE),
150 fStopped(t.fStopped),
151 fLhElectron(0),
152 fNWrong(t.fNWrong),
153 fNRotate(t.fNRotate),
154 fNCross(t.fNCross),
155 fNExpected(t.fNExpected),
156 fNLast(t.fNLast),
157 fNExpectedLast(t.fNExpectedLast),
158 fNdedx(t.fNdedx),
159 fChi2Last(t.fChi2Last),
160 fBackupTrack(0x0)
53c17fbf 161{
46d29e70 162 //
163 // Copy constructor.
164 //
6c94f330 165 for (Int_t i=0;i<kNplane;i++){
166 for (Int_t j=0; j<kNslice; j++) {
6d45eaef 167 fdEdxPlane[i][j] = t.fdEdxPlane[i][j];
168 }
169 fTimBinPlane[i] = t.fTimBinPlane[i];
170 fTracklets[i] = t.fTracklets[i];
eab5961e 171 }
46d29e70 172
6c94f330 173 Int_t n=t.GetNumberOfClusters();
174 SetNumberOfClusters(n);
175 for (Int_t i=0; i<n; i++) {
176 fIndex[i]=t.fIndex[i];
177 fIndexBackup[i]=t.fIndex[i]; // MI - backup indexes
178 fdQdl[i]=t.fdQdl[i];
03b0452e 179 }
15ed8ba1 180
6c94f330 181 // initialisation (i starts from n) [SR, GSI, 18.02.2003]
182 for(UInt_t i=n; i<kMAXCLUSTERSPERTRACK; i++) {
183 fdQdl[i] = 0;
184 fIndex[i] = 0;
185 fIndexBackup[i] = 0; //MI backup indexes
15ed8ba1 186 }
6c94f330 187 for (Int_t i=0;i<3;i++) fBudget[i]=t.fBudget[i];
5443e65e 188}
189
190//_____________________________________________________________________________
6c94f330 191AliTRDtrack::AliTRDtrack(const AliKalmanTrack& t, Double_t /*alpha*/):
192 AliKalmanTrack(t),
193 fSeedLab(-1),
194 fdEdx(t.GetPIDsignal()),
195 fdEdxT(0),
196 fDE(0),
197 fStopped(kFALSE),
198 fLhElectron(0.0),
199 fNWrong(0),
200 fNRotate(0),
201 fNCross(0),
202 fNExpected(0),
203 fNLast(0),
204 fNExpectedLast(0),
205 fNdedx(0),
206 fChi2Last(0.0),
207 fBackupTrack(0x0)
53c17fbf 208{
5443e65e 209 //
210 // Constructor from AliTPCtrack or AliITStrack .
211 //
212
6c94f330 213 SetLabel(t.GetLabel());
214 SetChi2(0.);
215 SetMass(t.GetMass());
5443e65e 216 SetNumberOfClusters(0);
217
6c94f330 218 for (Int_t i=0;i<kNplane;i++){
219 for (Int_t j=0;j<kNslice;j++){
6d45eaef 220 fdEdxPlane[i][j] = 0.0;
221 }
eab5961e 222 fTimBinPlane[i] = -1;
223 }
5443e65e 224
6c94f330 225 // Initialization [SR, GSI, 18.02.2003]
226 for(UInt_t i=0; i<kMAXCLUSTERSPERTRACK; i++) {
227 fdQdl[i] = 0;
228 fIndex[i] = 0;
229 fIndexBackup[i] = 0; // MI backup indexes
79e94bf8 230 }
4f1c04d3 231
6c94f330 232 for (Int_t i=0;i<3;i++) { fBudget[i]=0;};
79e94bf8 233}
53c17fbf 234
79e94bf8 235//_____________________________________________________________________________
6c94f330 236AliTRDtrack::AliTRDtrack(const AliESDtrack &t):
237 AliKalmanTrack(),
238 fSeedLab(-1),
239 fdEdx(t.GetTRDsignal()),
240 fdEdxT(0),
241 fDE(0),
242 fStopped(kFALSE),
243 fLhElectron(0),
244 fNWrong(0),
245 fNRotate(0),
246 fNCross(0),
247 fNExpected(0),
248 fNLast(0),
249 fNExpectedLast(0),
250 fNdedx(0),
251 fChi2Last(1e10),
252 fBackupTrack(0x0)
53c17fbf 253{
79e94bf8 254 //
255 // Constructor from AliESDtrack
256 //
79e94bf8 257 SetLabel(t.GetLabel());
258 SetChi2(0.);
259 SetMass(t.GetMass());
1e9bb598 260 SetNumberOfClusters(t.GetTRDclusters(fIndex));
15ed8ba1 261
46e2d86c 262 Int_t ncl = t.GetTRDclusters(fIndexBackup);
6c94f330 263 for (UInt_t i=ncl;i<kMAXCLUSTERSPERTRACK;i++) {
264 fIndexBackup[i]=0;
265 fIndex[i] = 0; //MI store indexes
15ed8ba1 266 }
6c94f330 267 for (Int_t i=0;i<kNplane;i++){
268 for (Int_t j=0;j<kNslice;j++){
6d45eaef 269 fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
270 }
eab5961e 271 fTimBinPlane[i] = t.GetTRDTimBin(i);
272 }
79e94bf8 273
79e94bf8 274
6c94f330 275 const AliExternalTrackParam *par=&t;
276 if (t.GetStatus()&AliESDtrack::kTRDbackup) {
277 par=t.GetOuterParam();
278 if (!par) {AliError("***** No backup info !!! ****"); par=&t;}
c5a8e3df 279 }
6c94f330 280 Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
79e94bf8 281
6c94f330 282
283 for (UInt_t i=0; i<kMAXCLUSTERSPERTRACK; i++) fdQdl[i] = 0;
5443e65e 284
6c94f330 285 for (Int_t i=0;i<3;i++) fBudget[i]=0;
c630aafd 286
6c94f330 287 if ((t.GetStatus()&AliESDtrack::kTIME) == 0) return;
c630aafd 288 StartTimeIntegral();
6c94f330 289 Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
c630aafd 290 SetIntegratedLength(t.GetIntegratedLength());
291
16d9fbba 292}
293
53c17fbf 294//____________________________________________________________________________
295AliTRDtrack::~AliTRDtrack()
3fad3d32 296{
297 //
53c17fbf 298 // Destructor
299 //
3fad3d32 300
6c94f330 301 if (fBackupTrack) delete fBackupTrack;
53c17fbf 302 fBackupTrack = 0;
3fad3d32 303
53c17fbf 304}
305
306//____________________________________________________________________________
53c17fbf 307Float_t AliTRDtrack::StatusForTOF()
7ad19338 308{
53c17fbf 309 //
310 // Defines the status of the TOF extrapolation
311 //
03b0452e 312
6c94f330 313 Float_t res = (0.2 + 0.8*(fN/(fNExpected+5.)))*(0.4+0.6*fTracklets[5].GetN()/20.);
314 res *= (0.25+0.8*40./(40.+fBudget[2]));
03b0452e 315 return res;
316
6c94f330 317 Int_t status=0;
318 if (GetNumberOfClusters()<20) return 0; //
319 if (fN>110&&fChi2/(Float_t(fN))<3) return 3; //gold
320 if (fNLast>30&&fChi2Last/(Float_t(fNLast))<3) return 3; //gold
321 if (fNLast>20&&fChi2Last/(Float_t(fNLast))<2) return 3; //gold
322 if (fNLast/(fNExpectedLast+3.)>0.8 && fChi2Last/Float_t(fNLast)<5&&fNLast>20) return 2; //silber
323 if (fNLast>5 &&((fNLast+1.)/(fNExpectedLast+1.))>0.8&&fChi2Last/(fNLast-5.)<6) return 1;
53c17fbf 324
4f1c04d3 325 return status;
53c17fbf 326
4f1c04d3 327}
46d29e70 328
329//_____________________________________________________________________________
53c17fbf 330Int_t AliTRDtrack::Compare(const TObject *o) const
331{
332 //
333 // Compares tracks according to their Y2 or curvature
334 //
46d29e70 335
d49e463b 336 AliTRDtrack *t = (AliTRDtrack *) o;
46d29e70 337 // Double_t co=t->GetSigmaY2();
338 // Double_t c =GetSigmaY2();
339
d49e463b 340 Double_t co = TMath::Abs(t->GetC());
341 Double_t c = TMath::Abs(GetC());
46d29e70 342
d49e463b 343 if (c > co) {
344 return 1;
345 }
346 else if (c < co) {
347 return -1;
348 }
46d29e70 349 return 0;
53c17fbf 350
46d29e70 351}
352
353//_____________________________________________________________________________
15ed8ba1 354void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
355{
356 //
d49e463b 357 // Calculates the truncated dE/dx within the "low" and "up" cuts.
15ed8ba1 358 //
359
d49e463b 360 Int_t i = 0;
361
362 // Array to sort the dEdx values according to amplitude
363 Float_t sorted[kMAXCLUSTERSPERTRACK];
15ed8ba1 364
365 // Number of clusters used for dedx
d49e463b 366 Int_t nc = fNdedx;
367
15ed8ba1 368 // Require at least 10 clusters for a dedx measurement
15ed8ba1 369 if (nc < 10) {
4f1c04d3 370 SetdEdx(0);
371 return;
372 }
a819a5f7 373
15ed8ba1 374 // Lower and upper bound
375 Int_t nl = Int_t(low * nc);
376 Int_t nu = Int_t( up * nc);
377
d49e463b 378 // Can fdQdl be negative ????
15ed8ba1 379 for (i = 0; i < nc; i++) {
380 sorted[i] = TMath::Abs(fdQdl[i]);
381 }
382
383 // Sort the dedx values by amplitude
384 Int_t *index = new Int_t[nc];
385 TMath::Sort(nc,sorted,index,kFALSE);
d49e463b 386
15ed8ba1 387 // Sum up the truncated charge between nl and nu
d49e463b 388 Float_t dedx = 0.0;
15ed8ba1 389 for (i = nl; i <= nu; i++) {
390 dedx += sorted[index[i]];
391 }
d49e463b 392 dedx /= (nu - nl + 1.0);
3551db50 393 SetdEdx(dedx);
394
a819a5f7 395}
396
a819a5f7 397//_____________________________________________________________________________
6c94f330 398Bool_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho)
46d29e70 399{
400 // Propagates a track of particle with mass=pm to a reference plane
401 // defined by x=xk through media of density=rho and radiationLength=x0
15ed8ba1 402
6c94f330 403 if (xk == GetX()) return kTRUE;
15ed8ba1 404
6c94f330 405 Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ();
46d29e70 406
6c94f330 407 Double_t bz=GetBz();
408 if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE;
9c9d2487 409
6c94f330 410 Double_t x=GetX(), y=GetY(), z=GetZ();
411 Double_t d=TMath::Sqrt((x-oldX)*(x-oldX)+(y-oldY)*(y-oldY)+(z-oldZ)*(z-oldZ));
412 if (oldX < xk)
413 if (IsStartedTimeIntegral()) {
414 Double_t l2=d;
415 Double_t crv=GetC();
416 if (TMath::Abs(l2*crv)>0.0001){
417 // make correction for curvature if neccesary
418 l2 = 0.5*TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY));
419 l2 = 2*TMath::ASin(l2*crv)/crv;
420 l2 = TMath::Sqrt(l2*l2+(z-oldZ)*(z-oldZ));
421 }
422 AddTimeStep(l2);
46d29e70 423 }
15ed8ba1 424
6c94f330 425 Double_t ll = (oldX < xk) ? -d : d;
426 if (!AliExternalTrackParam::CorrectForMaterial(ll*rho/x0,x0,GetMass()))
427 return kFALSE;
15ed8ba1 428
6c94f330 429 {//Energy losses************************
430 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
431 Double_t beta2=p2/(p2 + GetMass()*GetMass());
432 if ((5940*beta2/(1-beta2+1e-10) - beta2) < 0) return kFALSE;
433
434 Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2+1e-10)) - beta2)*d*rho;
435 Float_t budget = d*rho;
436 fBudget[0] += budget;
437 /*
438 // suspicious part - think about it ?
96c3a73c 439 Double_t kinE = TMath::Sqrt(p2);
6c94f330 440 if (dE>0.8*kinE) dE = 0.8*kinE; //
441 if (dE<0) dE = 0.0; // not valid region for Bethe bloch
442 */
443 //
444 fDE+=dE;
445 /*
446 // Suspicious ! I.B.
447 Double_t sigmade = 0.07*TMath::Sqrt(TMath::Abs(dE)); // energy loss fluctuation
448 Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
3fad3d32 449 fCcc += sigmac2;
6c94f330 450 fCee += fX*fX*sigmac2;
451 */
0d5b5c27 452 }
453
6c94f330 454 return kTRUE;
6d45eaef 455
46d29e70 456}
457
46d29e70 458//_____________________________________________________________________________
6c94f330 459Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index,
460 Double_t h01)
46d29e70 461{
6c94f330 462 // Assignes found cluster to the track and updates track information
46d29e70 463
b8dc2353 464 Bool_t fNoTilt = kTRUE;
6c94f330 465 if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
466 // add angular effect to the error contribution - MI
467 Float_t tangent2 = GetSnp()*GetSnp();
15ed8ba1 468 if (tangent2 < 0.90000){
6c94f330 469 tangent2 = tangent2/(1.-tangent2);
b8dc2353 470 }
6c94f330 471 //Float_t errang = tangent2*0.04; //
b8dc2353 472
6c94f330 473 Double_t p[2]={c->GetY(), c->GetZ()};
474 //Double_t cov[3]={c->GetSigmaY2()+errang, 0., c->GetSigmaZ2()*100.};
475 Double_t sy2=c->GetSigmaY2()*4;
476 Double_t sz2=c->GetSigmaZ2()*4;
477 Double_t cov[3]={sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2};
46e2d86c 478
6c94f330 479 if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
46d29e70 480
6c94f330 481 Int_t n=GetNumberOfClusters();
482 fIndex[n]=index;
b8dc2353 483 SetNumberOfClusters(n+1);
fd621f36 484
b8dc2353 485 SetChi2(GetChi2()+chisq);
5443e65e 486
6c94f330 487 return kTRUE;
46d29e70 488}
53c17fbf 489
46e2d86c 490//_____________________________________________________________________________
6c94f330 491Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, Int_t index, Double_t h01, Int_t /*plane*/) {
46e2d86c 492 // Assignes found cluster to the track and updates track information
46e2d86c 493 Bool_t fNoTilt = kTRUE;
6c94f330 494 if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
495 // add angular effect to the error contribution and make correction - MI
3c625a9b 496 //
6c94f330 497 Double_t tangent2 = GetSnp()*GetSnp();
498 if (tangent2 < 0.90000){
499 tangent2 = tangent2/(1.-tangent2);
15ed8ba1 500 }
6c94f330 501 Double_t tangent = TMath::Sqrt(tangent2);
502 if (GetSnp()<0) tangent*=-1;
503 // Double_t correction = 0*plane;
504 /*
505 Double_t errang = tangent2*0.04; //
506 Double_t errsys =0.025*0.025*20; //systematic error part
15ed8ba1 507
6c94f330 508 Float_t extend =1;
509 if (c->GetNPads()==4) extend=2;
510 */
511 //if (c->GetNPads()==5) extend=3;
512 //if (c->GetNPads()==6) extend=3;
513 //if (c->GetQ()<15) return 1;
15ed8ba1 514
c5a8e3df 515 /*
6c94f330 516 if (corrector!=0){
3c625a9b 517 //if (0){
518 correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent);
519 if (TMath::Abs(correction)>0){
520 //if we have info
521 errang = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent);
522 errang *= errang;
523 errang += tangent2*0.04;
524 }
525 }
c5a8e3df 526 */
6c94f330 527 //
528 //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
529 /*
530 {
531 Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ();
532 printf("%e %e %e %e\n",dy,dz,padlength/2,h01);
c84a5e9e 533 }
6c94f330 534 */
535 Double_t p[2]={c->GetY(), c->GetZ()};
536 /*
537 Double_t cov[3]={(c->GetSigmaY2()+errang+errsys)*extend, 0.,
538 c->GetSigmaZ2()*10000.};
539 */
540 Double_t sy2=c->GetSigmaY2()*4;
541 Double_t sz2=c->GetSigmaZ2()*4;
542 Double_t cov[3]={sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2};
15ed8ba1 543
6c94f330 544 if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
46e2d86c 545
6c94f330 546 Int_t n=GetNumberOfClusters();
547 fIndex[n]=index;
46e2d86c 548 SetNumberOfClusters(n+1);
6c94f330 549 SetChi2(GetChi2()+chisq);
46e2d86c 550
6c94f330 551 return kTRUE;
53c17fbf 552}
6c94f330 553/*
7ad19338 554//_____________________________________________________________________________
555Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
556{
557 //
558 // Assignes found tracklet to the track and updates track information
559 //
6c94f330 560 //
561 Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
562 r00+=fCyy; r01+=fCzy; r11+=fCzz;
563 //
564 Double_t det=r00*r11 - r01*r01;
565 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
566 //
7ad19338 567
6c94f330 568 Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
7ad19338 569
15ed8ba1 570
6c94f330 571 Double_t s00 = tracklet.GetTrackletSigma2(); // error pad
572 Double_t s11 = 100000; // error pad-row
573 Float_t h01 = tracklet.GetTilt();
574 //
575 // r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
576 r00 = fCyy + fCzz*h01*h01+s00;
577 // r01 = fCzy + fCzz*h01;
578 r01 = fCzy ;
7ad19338 579 r11 = fCzz + s11;
580 det = r00*r11 - r01*r01;
6c94f330 581 // inverse matrix
582 tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
583
584 Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
585 Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
586 Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
587 Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
588 Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
15ed8ba1 589
6c94f330 590 // K matrix
591// k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
592// k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
593// k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
594// k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
595// k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);
596 //
597 //Update measurement
598 Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;
599 // cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
600 if (TMath::Abs(cur*fX-eta) >= 0.90000) {
601 //Int_t n=GetNumberOfClusters();
602 // if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
7ad19338 603 return 0;
604 }
6c94f330 605// k01+=h01*k00;
606// k11+=h01*k10;
607// k21+=h01*k20;
608// k31+=h01*k30;
609// k41+=h01*k40;
610
7ad19338 611
612 fY += k00*dy + k01*dz;
613 fZ += k10*dy + k11*dz;
614 fE = eta;
615 fT += k30*dy + k31*dz;
616 fC = cur;
617
6c94f330 618
619 //Update covariance
620 //
621 //
622 Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
623 Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
624 Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
625 //Double_t oldte = fCte, oldce = fCce;
7ad19338 626 //Double_t oldct = fCct;
627
6c94f330 628 fCyy-=k00*oldyy+k01*oldzy;
629 fCzy-=k10*oldyy+k11*oldzy;
630 fCey-=k20*oldyy+k21*oldzy;
631 fCty-=k30*oldyy+k31*oldzy;
632 fCcy-=k40*oldyy+k41*oldzy;
633 //
634 fCzz-=k10*oldzy+k11*oldzz;
635 fCez-=k20*oldzy+k21*oldzz;
636 fCtz-=k30*oldzy+k31*oldzz;
637 fCcz-=k40*oldzy+k41*oldzz;
638 //
639 fCee-=k20*oldey+k21*oldez;
640 fCte-=k30*oldey+k31*oldez;
641 fCce-=k40*oldey+k41*oldez;
642 //
643 fCtt-=k30*oldty+k31*oldtz;
644 fCct-=k40*oldty+k41*oldtz;
645 //
646 fCcc-=k40*oldcy+k41*oldcz;
647 //
15ed8ba1 648
6c94f330 649 //Int_t n=GetNumberOfClusters();
650 //fIndex[n]=index;
651 //SetNumberOfClusters(n+1);
652
653 //SetChi2(GetChi2()+chisq);
654 // cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
7ad19338 655
53c17fbf 656 return 1;
7ad19338 657
53c17fbf 658}
6c94f330 659*/
c84a5e9e 660
46d29e70 661//_____________________________________________________________________________
6c94f330 662Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
46d29e70 663{
6c94f330 664 // Rotates track parameters in R*phi plane
665 // if absolute rotation alpha is in global system
666 // otherwise alpha rotation is relative to the current rotation angle
667
3fad3d32 668 if (absolute) {
6c94f330 669 alpha -= GetAlpha();
3fad3d32 670 }
671 else{
672 fNRotate++;
673 }
46d29e70 674
6c94f330 675 return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
53c17fbf 676}
7ad19338 677
46d29e70 678//_____________________________________________________________________________
fd621f36 679Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
46d29e70 680{
53c17fbf 681 //
682 // Returns the track chi2
683 //
684
6c94f330 685 Double_t p[2]={c->GetY(), c->GetZ()};
686 Double_t sy2=c->GetSigmaY2()*4;
687 Double_t sz2=c->GetSigmaZ2()*4;
688 Double_t cov[3]={sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2};
689
690 return AliExternalTrackParam::GetPredictedChi2(p,cov);
15ed8ba1 691
6c94f330 692 /*
693 Bool_t fNoTilt = kTRUE;
694 if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
15ed8ba1 695
6c94f330 696 return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2();
697 */
15ed8ba1 698
6c94f330 699 /*
700 Double_t chi2, dy, r00, r01, r11;
b8dc2353 701
6c94f330 702 if(fNoTilt) {
703 dy=c->GetY() - fY;
704 r00=c->GetSigmaY2();
705 chi2 = (dy*dy)/r00;
46d29e70 706 }
b8dc2353 707 else {
6c94f330 708 Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
709 //
710 r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2();
711 r00+=fCyy; r01+=fCzy; r11+=fCzz;
b8dc2353 712
6c94f330 713 Double_t det=r00*r11 - r01*r01;
b8dc2353 714 if (TMath::Abs(det) < 1.e-10) {
6c94f330 715 Int_t n=GetNumberOfClusters();
716 if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n";
b8dc2353 717 return 1e10;
718 }
6c94f330 719 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
720 Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
4f1c04d3 721 Double_t tiltdz = dz;
6c94f330 722 if (TMath::Abs(tiltdz)>padlength/2.) {
723 tiltdz = TMath::Sign(padlength/2,dz);
4f1c04d3 724 }
6c94f330 725 // dy=dy+h01*dz;
726 dy=dy+h01*tiltdz;
a819a5f7 727
6c94f330 728 chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
b8dc2353 729 }
53c17fbf 730
b8dc2353 731 return chi2;
6c94f330 732 */
53c17fbf 733}
7ad19338 734
53c17fbf 735//_____________________________________________________________________________
16d9fbba 736void AliTRDtrack::MakeBackupTrack()
737{
738 //
53c17fbf 739 // Creates a backup track
16d9fbba 740 //
53c17fbf 741
6c94f330 742 if (fBackupTrack) delete fBackupTrack;
16d9fbba 743 fBackupTrack = new AliTRDtrack(*this);
4f1c04d3 744
745}
746
53c17fbf 747//_____________________________________________________________________________
748Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
749{
4f1c04d3 750 //
751 // Find prolongation at given x
6c94f330 752 // return 0 if not exist
4f1c04d3 753
6c94f330 754 Double_t bz=GetBz();
15ed8ba1 755
6c94f330 756 if (!AliExternalTrackParam::GetYAt(xk,bz,y)) return 0;
757 if (!AliExternalTrackParam::GetZAt(xk,bz,z)) return 0;
4f1c04d3 758
6c94f330 759 return 1;
16d9fbba 760}
3fad3d32 761
53c17fbf 762//_____________________________________________________________________________
6c94f330 763Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
3fad3d32 764{
765 //
6c94f330 766 // Propagate track to given x position
767 // works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF)
3fad3d32 768 //
6c94f330 769 // material budget from geo manager
770 //
771 Double_t xyz0[3], xyz1[3],y,z;
772 const Double_t kAlphac = TMath::Pi()/9.;
15ed8ba1 773 const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
6c94f330 774 // critical alpha - cross sector indication
775 //
776 Double_t dir = (GetX()>xr) ? -1.:1.;
777 // direction +-
778 for (Double_t x=GetX()+dir*step;dir*x<dir*xr;x+=dir*step){
779 //
780 GetXYZ(xyz0);
3fad3d32 781 GetProlongation(x,y,z);
6c94f330 782 xyz1[0] = x*TMath::Cos(GetAlpha())+y*TMath::Sin(GetAlpha());
783 xyz1[1] = x*TMath::Sin(GetAlpha())-y*TMath::Cos(GetAlpha());
3fad3d32 784 xyz1[2] = z;
785 Double_t param[7];
786 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
6c94f330 787 //
788 if (param[0]>0&&param[1]>0) PropagateTo(x,param[1],param[0]);
789 if (GetY()>GetX()*kTalphac){
53c17fbf 790 Rotate(-kAlphac);
3fad3d32 791 }
6c94f330 792 if (GetY()<-GetX()*kTalphac){
53c17fbf 793 Rotate(kAlphac);
3fad3d32 794 }
795 }
6c94f330 796 //
3fad3d32 797 PropagateTo(xr);
53c17fbf 798
3fad3d32 799 return 0;
3fad3d32 800
53c17fbf 801}
3fad3d32 802
53c17fbf 803//_____________________________________________________________________________
6c94f330 804Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step)
3fad3d32 805{
806 //
6c94f330 807 // propagate track to the radial position
808 // rotation always connected to the last track position
3fad3d32 809 //
6c94f330 810 Double_t xyz0[3], xyz1[3],y,z;
811 Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
812 Double_t dir = (radius>r) ? -1.:1.; // direction +-
813 //
814 for (Double_t x=radius+dir*step;dir*x<dir*r;x+=dir*step){
815 GetXYZ(xyz0);
3fad3d32 816 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
817 Rotate(alpha,kTRUE);
6c94f330 818 GetXYZ(xyz0);
3fad3d32 819 GetProlongation(x,y,z);
6c94f330 820 xyz1[0] = x*TMath::Cos(alpha)+y*TMath::Sin(alpha);
821 xyz1[1] = x*TMath::Sin(alpha)-y*TMath::Cos(alpha);
3fad3d32 822 xyz1[2] = z;
823 Double_t param[7];
824 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
6c94f330 825 if (param[1]<=0) param[1] =100000000;
3fad3d32 826 PropagateTo(x,param[1],param[0]);
827 }
6c94f330 828 GetXYZ(xyz0);
3fad3d32 829 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
830 Rotate(alpha,kTRUE);
6c94f330 831 GetXYZ(xyz0);
3fad3d32 832 GetProlongation(r,y,z);
6c94f330 833 xyz1[0] = r*TMath::Cos(alpha)+y*TMath::Sin(alpha);
834 xyz1[1] = r*TMath::Sin(alpha)-y*TMath::Cos(alpha);
3fad3d32 835 xyz1[2] = z;
836 Double_t param[7];
837 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
6c94f330 838 //
839 if (param[1]<=0) param[1] =100000000;
3fad3d32 840 PropagateTo(r,param[1],param[0]);
53c17fbf 841
3fad3d32 842 return 0;
53c17fbf 843
844}
845
846//_____________________________________________________________________________
4e28c495 847Int_t AliTRDtrack::GetSector() const
53c17fbf 848{
849 //
850 // Return the current sector
851 //
852
6c94f330 853 return Int_t(TVector2::Phi_0_2pi(GetAlpha())
53c17fbf 854 / AliTRDgeometry::GetAlpha())
855 % AliTRDgeometry::kNsect;
856
3fad3d32 857}
858
53c17fbf 859//_____________________________________________________________________________
4e28c495 860void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
53c17fbf 861{
862 //
863 // The sampled energy loss
864 //
865
866 Double_t s = GetSnp();
867 Double_t t = GetTgl();
6c94f330 868 q *= TMath::Sqrt((1-s*s)/(1+t*t));
53c17fbf 869 fdQdl[i] = q;
870
871}
872
873 //_____________________________________________________________________________
4e28c495 874void AliTRDtrack::SetSampledEdx(Float_t q)
53c17fbf 875{
876 //
877 // The sampled energy loss
878 //
879
880 Double_t s = GetSnp();
881 Double_t t = GetTgl();
6c94f330 882 q*= TMath::Sqrt((1-s*s)/(1+t*t));
53c17fbf 883 fdQdl[fNdedx] = q;
884 fNdedx++;
885
886}
887
6c94f330 888Double_t AliTRDtrack::GetBz() const {
15ed8ba1 889 //
6c94f330 890 // returns Bz component of the magnetic field (kG)
15ed8ba1 891 //
6c94f330 892 if (AliTracker::UniformField()) return AliTracker::GetBz();
893 Double_t r[3]; GetXYZ(r);
894 return AliTracker::GetBz(r);
895}
53c17fbf 896
53c17fbf 897