Added ignores
[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 * *
69dee759 7 * Permission to use, copy, modify and distribute this software and its *
46d29e70 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
e1e6896f 18#include <TVector2.h>
46d29e70 19
6c94f330 20#include "AliTracker.h"
53c17fbf 21#include "AliESDtrack.h"
fe33d239 22
46d29e70 23#include "AliTRDgeometry.h"
24#include "AliTRDcluster.h"
25#include "AliTRDtrack.h"
43bfc8af 26#include "AliTRDcalibDB.h"
720a0a16 27#include "Cal/AliTRDCalPID.h"
43bfc8af 28
46d29e70 29ClassImp(AliTRDtrack)
30
53c17fbf 31///////////////////////////////////////////////////////////////////////////////
32// //
33// Represents a reconstructed TRD track //
34// Local TRD Kalman track //
35// //
36///////////////////////////////////////////////////////////////////////////////
7ad19338 37
fe33d239 38//_____________________________________________________________________________
39AliTRDtrack::AliTRDtrack()
40 :AliKalmanTrack()
41 ,fSeedLab(-1)
42 ,fdEdx(0)
43 ,fDE(0)
0906e73e 44 ,fPIDquality(0)
79143cbf 45 ,fClusterOwner(kFALSE)
69dee759 46 ,fPIDmethod(kLQ)
fe33d239 47 ,fStopped(kFALSE)
48 ,fLhElectron(0)
49 ,fNWrong(0)
50 ,fNRotate(0)
51 ,fNCross(0)
52 ,fNExpected(0)
53 ,fNLast(0)
54 ,fNExpectedLast(0)
55 ,fNdedx(0)
56 ,fChi2Last(1e10)
57 ,fBackupTrack(0x0)
f146da82 58{
fe33d239 59 //
60 // AliTRDtrack default constructor
61 //
62
63 for (Int_t i = 0; i < kNplane; i++) {
64 for (Int_t j = 0; j < kNslice; j++) {
65 fdEdxPlane[i][j] = 0.0;
6d45eaef 66 }
f146da82 67 fTimBinPlane[i] = -1;
af26ce80 68 fMom[i] = -1.;
69 fSnp[i] = 0.;
70 fTgl[i] = 0.;
71 }
fe33d239 72
73 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
74 fIndex[i] = 0;
6c94f330 75 fIndexBackup[i] = 0;
fe33d239 76 fdQdl[i] = 0;
af26ce80 77 fClusters[i] = 0x0;
78 }
fe33d239 79
80 for (Int_t i = 0; i < 3; i++) {
81 fBudget[i] = 0;
f146da82 82 }
fe33d239 83
f146da82 84}
6d45eaef 85
46d29e70 86//_____________________________________________________________________________
43bfc8af 87AliTRDtrack::AliTRDtrack(AliTRDcluster *c, Int_t index
fe33d239 88 , const Double_t p[5], const Double_t cov[15]
89 , Double_t x, Double_t alpha)
90 :AliKalmanTrack()
91 ,fSeedLab(-1)
92 ,fdEdx(0)
93 ,fDE(0)
0906e73e 94 ,fPIDquality(0)
79143cbf 95 ,fClusterOwner(kFALSE)
69dee759 96 ,fPIDmethod(kLQ)
fe33d239 97 ,fStopped(kFALSE)
98 ,fLhElectron(0)
99 ,fNWrong(0)
100 ,fNRotate(0)
101 ,fNCross(0)
102 ,fNExpected(0)
103 ,fNLast(0)
104 ,fNExpectedLast(0)
105 ,fNdedx(0)
106 ,fChi2Last(1e10)
43bfc8af 107 ,fBackupTrack(0x0)
15ed8ba1 108{
fe33d239 109 //
110 // The main AliTRDtrack constructor.
111 //
112
af26ce80 113 Double_t cnv = 1.0 / (GetBz() * kB2C);
fe33d239 114
115 Double_t pp[5] = { p[0]
116 , p[1]
117 , x*p[4] - p[2]
118 , p[3]
119 , p[4]*cnv };
6c94f330 120
af26ce80 121 Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
122 Double_t c32 = x*cov[13] - cov[ 8];
123 Double_t c20 = x*cov[10] - cov[ 3];
124 Double_t c21 = x*cov[11] - cov[ 4];
125 Double_t c42 = x*cov[14] - cov[12];
6c94f330 126
af26ce80 127 Double_t cc[15] = { cov[ 0]
128 , cov[ 1], cov[ 2]
fe33d239 129 , c20, c21, c22
af26ce80 130 , cov[ 6], cov[ 7], c32, cov[ 9]
fe33d239 131 , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
6c94f330 132
133 Set(x,alpha,pp,cc);
5443e65e 134 SetNumberOfClusters(1);
af26ce80 135 fIndex[0] = index;
136 fClusters[0] = c;
5443e65e 137
fe33d239 138 for (Int_t i = 0; i < kNplane; i++) {
139 for (Int_t j = 0; j < kNslice; j++) {
140 fdEdxPlane[i][j] = 0.0;
6d45eaef 141 }
142 fTimBinPlane[i] = -1;
af26ce80 143 fMom[i] = -1.;
144 fSnp[i] = 0.;
145 fTgl[i] = 0.;
eab5961e 146 }
a2b90f83 147
5443e65e 148 Double_t q = TMath::Abs(c->GetQ());
fe33d239 149 Double_t s = GetSnp();
150 Double_t t = GetTgl();
151 if (s*s < 1) {
152 q *= TMath::Sqrt((1-s*s)/(1+t*t));
153 }
53c17fbf 154
af26ce80 155 fdQdl[0] = q;
fe33d239 156 for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) {
157 fdQdl[i] = 0;
158 fIndex[i] = 0;
159 fIndexBackup[i] = 0;
af26ce80 160 fClusters[i] = 0x0;
15ed8ba1 161 }
fe33d239 162
163 for (Int_t i = 0; i < 3;i++) {
164 fBudget[i] = 0;
165 }
166
46d29e70 167}
168
169//_____________________________________________________________________________
43bfc8af 170AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/)
fe33d239 171 :AliKalmanTrack(t)
172 ,fSeedLab(t.GetSeedLabel())
173 ,fdEdx(t.fdEdx)
174 ,fDE(t.fDE)
0906e73e 175 ,fPIDquality(t.fPIDquality)
79143cbf 176 ,fClusterOwner(kTRUE)
69dee759 177 ,fPIDmethod(t.fPIDmethod)
fe33d239 178 ,fStopped(t.fStopped)
179 ,fLhElectron(0)
180 ,fNWrong(t.fNWrong)
181 ,fNRotate(t.fNRotate)
182 ,fNCross(t.fNCross)
183 ,fNExpected(t.fNExpected)
184 ,fNLast(t.fNLast)
185 ,fNExpectedLast(t.fNExpectedLast)
186 ,fNdedx(t.fNdedx)
187 ,fChi2Last(t.fChi2Last)
43bfc8af 188 ,fBackupTrack(0x0)
53c17fbf 189{
46d29e70 190 //
191 // Copy constructor.
192 //
fe33d239 193
194 for (Int_t i = 0; i < kNplane; i++) {
195 for (Int_t j = 0; j < kNslice; j++) {
6d45eaef 196 fdEdxPlane[i][j] = t.fdEdxPlane[i][j];
197 }
198 fTimBinPlane[i] = t.fTimBinPlane[i];
199 fTracklets[i] = t.fTracklets[i];
af26ce80 200 fMom[i] = t.fMom[i];
201 fSnp[i] = t.fSnp[i];
202 fTgl[i] = t.fTgl[i];
eab5961e 203 }
46d29e70 204
fe33d239 205 Int_t n = t.GetNumberOfClusters();
6c94f330 206 SetNumberOfClusters(n);
fe33d239 207
208 for (Int_t i = 0; i < n; i++) {
209 fIndex[i] = t.fIndex[i];
210 fIndexBackup[i] = t.fIndex[i];
211 fdQdl[i] = t.fdQdl[i];
af26ce80 212 if (fClusterOwner && t.fClusters[i]) {
213 fClusters[i] = new AliTRDcluster(*(t.fClusters[i]));
214 }
215 else {
216 fClusters[i] = t.fClusters[i];
217 }
218 }
fe33d239 219
220 for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) {
221 fdQdl[i] = 0;
222 fIndex[i] = 0;
223 fIndexBackup[i] = 0;
af26ce80 224 fClusters[i] = 0x0;
03b0452e 225 }
15ed8ba1 226
fe33d239 227 for (Int_t i = 0; i < 3;i++) {
228 fBudget[i] = t.fBudget[i];
15ed8ba1 229 }
fe33d239 230
0906e73e 231 for(Int_t ispec = 0; ispec<AliPID::kSPECIES; ispec++) fPID[ispec] = t.fPID[ispec];
232}
5443e65e 233
234//_____________________________________________________________________________
fe33d239 235AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
236 :AliKalmanTrack(t)
237 ,fSeedLab(-1)
238 ,fdEdx(t.GetPIDsignal())
239 ,fDE(0)
0906e73e 240 ,fPIDquality(0)
79143cbf 241 ,fClusterOwner(kFALSE)
69dee759 242 ,fPIDmethod(kLQ)
fe33d239 243 ,fStopped(kFALSE)
244 ,fLhElectron(0.0)
245 ,fNWrong(0)
246 ,fNRotate(0)
247 ,fNCross(0)
248 ,fNExpected(0)
249 ,fNLast(0)
250 ,fNExpectedLast(0)
251 ,fNdedx(0)
252 ,fChi2Last(0.0)
253 ,fBackupTrack(0x0)
53c17fbf 254{
5443e65e 255 //
fe33d239 256 // Constructor from AliTPCtrack or AliITStrack
5443e65e 257 //
258
6c94f330 259 SetLabel(t.GetLabel());
fe33d239 260 SetChi2(0.0);
6c94f330 261 SetMass(t.GetMass());
5443e65e 262 SetNumberOfClusters(0);
263
fe33d239 264 for (Int_t i = 0; i < kNplane; i++) {
265 for (Int_t j = 0; j < kNslice; j++) {
6d45eaef 266 fdEdxPlane[i][j] = 0.0;
267 }
eab5961e 268 fTimBinPlane[i] = -1;
af26ce80 269 fMom[i] = -1.;
270 fSnp[i] = 0.;
271 fTgl[i] = 0.;
eab5961e 272 }
5443e65e 273
fe33d239 274 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
275 fdQdl[i] = 0;
276 fIndex[i] = 0;
277 fIndexBackup[i] = 0;
af26ce80 278 fClusters[i] = 0x0;
79e94bf8 279 }
4f1c04d3 280
fe33d239 281 for (Int_t i = 0; i < 3; i++) {
282 fBudget[i] = 0;
283 }
284
0906e73e 285}
53c17fbf 286
79e94bf8 287//_____________________________________________________________________________
fe33d239 288AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
289 :AliKalmanTrack()
290 ,fSeedLab(-1)
291 ,fdEdx(t.GetTRDsignal())
292 ,fDE(0)
0906e73e 293 ,fPIDquality(0)
79143cbf 294 ,fClusterOwner(kFALSE)
69dee759 295 ,fPIDmethod(kLQ)
fe33d239 296 ,fStopped(kFALSE)
297 ,fLhElectron(0)
298 ,fNWrong(0)
299 ,fNRotate(0)
300 ,fNCross(0)
301 ,fNExpected(0)
302 ,fNLast(0)
303 ,fNExpectedLast(0)
304 ,fNdedx(0)
305 ,fChi2Last(1e10)
306 ,fBackupTrack(0x0)
53c17fbf 307{
79e94bf8 308 //
309 // Constructor from AliESDtrack
310 //
fe33d239 311
79e94bf8 312 SetLabel(t.GetLabel());
fe33d239 313 SetChi2(0.0);
79e94bf8 314 SetMass(t.GetMass());
1e9bb598 315 SetNumberOfClusters(t.GetTRDclusters(fIndex));
15ed8ba1 316
46e2d86c 317 Int_t ncl = t.GetTRDclusters(fIndexBackup);
fe33d239 318 for (UInt_t i = ncl; i < kMAXCLUSTERSPERTRACK; i++) {
319 fIndexBackup[i] = 0;
320 fIndex[i] = 0;
15ed8ba1 321 }
fe33d239 322
323 for (Int_t i = 0; i < kNplane; i++) {
324 for (Int_t j = 0; j < kNslice; j++) {
6d45eaef 325 fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
326 }
eab5961e 327 fTimBinPlane[i] = t.GetTRDTimBin(i);
af26ce80 328 fMom[i] = -1.;
329 fSnp[i] = 0.;
330 fTgl[i] = 0.;
eab5961e 331 }
79e94bf8 332
fe33d239 333 const AliExternalTrackParam *par = &t;
334 if (t.GetStatus() & AliESDtrack::kTRDbackup) {
335 par = t.GetOuterParam();
336 if (!par) {
337 AliError("No backup info!");
338 par = &t;
339 }
c5a8e3df 340 }
69dee759 341 Set(par->GetX()
342 ,par->GetAlpha()
343 ,par->GetParameter()
344 ,par->GetCovariance());
79e94bf8 345
fe33d239 346 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
43bfc8af 347 fdQdl[i] = 0;
af26ce80 348 fClusters[i] = 0x0;
349 }
fe33d239 350
351 for (Int_t i = 0; i < 3; i++) {
352 fBudget[i] = 0;
353 }
5443e65e 354
fe33d239 355 if ((t.GetStatus() & AliESDtrack::kTIME) == 0) {
356 return;
357 }
c630aafd 358
c630aafd 359 StartTimeIntegral();
fe33d239 360 Double_t times[10];
361 t.GetIntegratedTimes(times);
362 SetIntegratedTimes(times);
c630aafd 363 SetIntegratedLength(t.GetIntegratedLength());
364
16d9fbba 365}
366
53c17fbf 367//____________________________________________________________________________
368AliTRDtrack::~AliTRDtrack()
3fad3d32 369{
370 //
53c17fbf 371 // Destructor
372 //
3fad3d32 373
fe33d239 374 if (fBackupTrack) {
375 delete fBackupTrack;
376 }
43bfc8af 377 fBackupTrack = 0x0;
69dee759 378
af26ce80 379 if (fClusterOwner) {
380 UInt_t icluster = 0;
381 while ((icluster < kMAXCLUSTERSPERTRACK) && fClusters[icluster]) {
382 delete fClusters[icluster];
383 fClusters[icluster] = 0x0;
384 icluster++;
385 }
43bfc8af 386 }
3fad3d32 387
53c17fbf 388}
389
390//____________________________________________________________________________
53c17fbf 391Float_t AliTRDtrack::StatusForTOF()
7ad19338 392{
53c17fbf 393 //
394 // Defines the status of the TOF extrapolation
395 //
03b0452e 396
fe33d239 397 // Definition of res ????
398 Float_t res = (0.2 + 0.8 * (fN / (fNExpected + 5.0)))
399 * (0.4 + 0.6 * fTracklets[5].GetN() / 20.0);
400 res *= (0.25 + 0.8 * 40.0 / (40.0 + fBudget[2]));
03b0452e 401 return res;
402
4f1c04d3 403}
46d29e70 404
405//_____________________________________________________________________________
53c17fbf 406Int_t AliTRDtrack::Compare(const TObject *o) const
407{
408 //
409 // Compares tracks according to their Y2 or curvature
410 //
46d29e70 411
d49e463b 412 AliTRDtrack *t = (AliTRDtrack *) o;
46d29e70 413
d49e463b 414 Double_t co = TMath::Abs(t->GetC());
415 Double_t c = TMath::Abs(GetC());
46d29e70 416
d49e463b 417 if (c > co) {
418 return 1;
419 }
420 else if (c < co) {
421 return -1;
422 }
69dee759 423
46d29e70 424 return 0;
53c17fbf 425
46d29e70 426}
427
428//_____________________________________________________________________________
43bfc8af 429void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
15ed8ba1 430{
431 //
d49e463b 432 // Calculates the truncated dE/dx within the "low" and "up" cuts.
15ed8ba1 433 //
434
d49e463b 435 // Array to sort the dEdx values according to amplitude
af26ce80 436 Float_t sorted[kMAXCLUSTERSPERTRACK];
69dee759 437 fdEdx = 0.0;
43bfc8af 438
15ed8ba1 439 // Require at least 10 clusters for a dedx measurement
69dee759 440 if (fNdedx < 10) {
441 return;
442 }
15ed8ba1 443
d49e463b 444 // Can fdQdl be negative ????
af26ce80 445 for (Int_t i = 0; i < fNdedx; i++) {
446 sorted[i] = TMath::Abs(fdQdl[i]);
447 }
15ed8ba1 448 // Sort the dedx values by amplitude
c6011b06 449 Int_t *index = new Int_t[fNdedx];
450 TMath::Sort(fNdedx, sorted, index, kFALSE);
451
43bfc8af 452 // Sum up the truncated charge between lower and upper bounds
c6011b06 453 Int_t nl = Int_t(low * fNdedx);
454 Int_t nu = Int_t( up * fNdedx);
af26ce80 455 for (Int_t i = nl; i <= nu; i++) {
456 fdEdx += sorted[index[i]];
457 }
458 fdEdx /= (nu - nl + 1.0);
459
460 delete[] index;
43bfc8af 461
43bfc8af 462}
463
464//_____________________________________________________________________________
44dbae42 465void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/)
43bfc8af 466{
467 //
468 // Set fdEdxPlane and fTimBinPlane and also get the
469 // Time bin for Max. Cluster
af26ce80 470 //
471 // Authors:
43bfc8af 472 // Prashant Shukla (shukla@physi.uni-heidelberg.de)
473 // Alexandru Bercuci (A.Bercuci@gsi.de)
af26ce80 474 //
43bfc8af 475
69dee759 476 // Max charge in chamber
477 Double_t maxcharge[AliESDtrack::kNPlane];
af26ce80 478 // Number of clusters attached to track per chamber and slice
479 Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
480 // Number of time bins in chamber
481 Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
482 Int_t plane; // Plane of current cluster
483 Int_t tb; // Time bin of current cluster
484 Int_t slice; // Current slice
485 AliTRDcluster *cluster = 0x0; // Pointer to current cluster
486
487 // Reset class and local counters/variables
43bfc8af 488 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
af26ce80 489 fTimBinPlane[iPlane] = -1;
69dee759 490 maxcharge[iPlane] = 0.0;
af26ce80 491 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
69dee759 492 fdEdxPlane[iPlane][iSlice] = 0.0;
43bfc8af 493 nCluster[iPlane][iSlice] = 0;
494 }
15ed8ba1 495 }
c6011b06 496
af26ce80 497 // Start looping over clusters attached to this track
498 for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
499
43bfc8af 500 cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
69dee759 501 if (!cluster) continue;
43bfc8af 502
af26ce80 503 // Read info from current cluster
504 plane = AliTRDgeometry::GetPlane(cluster->GetDetector());
43bfc8af 505 if (plane < 0 || plane >= AliESDtrack::kNPlane) {
506 AliError(Form("Wrong plane %d", plane));
507 continue;
508 }
509
af26ce80 510 tb = cluster->GetLocalTimeBin();
e4f2f73d 511 if ((tb < 0) || (tb >= ntb)) {
512 AliWarning(Form("time bin < 0 or > %d in cluster %d", ntb, iClus));
af26ce80 513 AliInfo(Form("dQ/dl %f", fdQdl[iClus]));
514 continue;
515 }
43bfc8af 516
af26ce80 517 slice = tb * AliESDtrack::kNSlice / ntb;
518
519 fdEdxPlane[plane][slice] += fdQdl[iClus];
520 if (fdQdl[iClus] > maxcharge[plane]) {
521 maxcharge[plane] = fdQdl[iClus];
43bfc8af 522 fTimBinPlane[plane] = tb;
523 }
43bfc8af 524
af26ce80 525 nCluster[plane][slice]++;
43bfc8af 526
af26ce80 527 } // End of loop over cluster
43bfc8af 528
529 // Normalize fdEdxPlane to number of clusters and set track segments
530 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
531 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
af26ce80 532 if (nCluster[iPlane][iSlice]) {
533 fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice];
534 }
43bfc8af 535 }
af26ce80 536 }
43bfc8af 537
af26ce80 538}
43bfc8af 539
540//_____________________________________________________________________________
69dee759 541void AliTRDtrack::CookdEdxNN(Float_t *dedx)
542{
543 //
44dbae42 544 // This function calcuates the input for the neural networks
545 // which are used for the PID.
69dee759 546 //
547
548 //number of time bins in chamber
549 Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
550
551 Int_t plane; // plane of current cluster
552 Int_t tb; // time bin of current cluster
553 Int_t slice; // curent slice
554 AliTRDcluster *cluster = 0x0; // pointer to current cluster
555 const Int_t kMLPscale = 16000; // scaling of the MLP input to be smaller than 1
556
557 // Reset class and local contors/variables
558 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++){
559 for (Int_t iSlice = 0; iSlice < kNMLPslice; iSlice++) {
560 *(dedx + (kNMLPslice * iPlane) + iSlice) = 0.0;
561 }
562 }
563
564 // Start looping over clusters attached to this track
565 for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
44dbae42 566
69dee759 567 cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
568 if (!cluster) {
569 continue;
570 }
44dbae42 571
69dee759 572 // Read info from current cluster
573 plane = AliTRDgeometry::GetPlane(cluster->GetDetector());
574 if (plane < 0 || plane >= AliESDtrack::kNPlane) {
575 AliError(Form("Wrong plane %d",plane));
576 continue;
577 }
578
579 tb = cluster->GetLocalTimeBin();
580 if (tb == 0 || tb >= ntb) {
581 AliWarning(Form("time bin 0 or > %d in cluster %d",ntb,iClus));
582 AliInfo(Form("dQ/dl %f",fdQdl[iClus]));
583 continue;
584 }
585
586 slice = tb * kNMLPslice / ntb;
44dbae42 587
69dee759 588 *(dedx+(kNMLPslice * plane) + slice) += fdQdl[iClus]/kMLPscale;
589
590 } // End of loop over cluster
44dbae42 591
69dee759 592}
44dbae42 593
594//_____________________________________________________________________________
69dee759 595void AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane)
43bfc8af 596{
af26ce80 597 //
598 // Set the momenta for a track segment in a given plane
599 //
600
601 if ((plane < 0) ||
602 (plane>= kNplane)) {
603 AliError(Form("Trying to access out of range plane (%d)", plane));
604 return;
605 }
43bfc8af 606
af26ce80 607 fSnp[plane] = GetSnp();
608 fTgl[plane] = GetTgl();
609 Double_t p[3];
610 GetPxPyPz(p);
611 fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
612
43bfc8af 613}
614
615//_____________________________________________________________________________
616Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const
617{
af26ce80 618 //
619 // Calculate the track length of a track segment in a given plane
620 //
43bfc8af 621
69dee759 622 if ((plane < 0) ||
623 (plane >= kNplane)) {
624 return 0.0;
625 }
43bfc8af 626
af26ce80 627 return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())
628 / TMath::Sqrt((1.0 - fSnp[plane]*fSnp[plane])
629 / (1.0 + fTgl[plane]*fTgl[plane]));
630
631}
43bfc8af 632
633//_____________________________________________________________________________
44dbae42 634Bool_t AliTRDtrack::CookPID(Int_t &pidQuality)
43bfc8af 635{
636 //
637 // This function calculates the PID probabilities based on TRD signals
638 //
639 // The method produces probabilities based on the charge
640 // and the position of the maximum time bin in each layer.
641 // The dE/dx information can be used as global charge or 2 to 3
720a0a16 642 // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
43bfc8af 643 // implementation.
644 //
645 // Author
646 // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
af26ce80 647 //
43bfc8af 648
af26ce80 649 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
650 if (!calibration) {
651 AliError("No access to calibration data");
44dbae42 652 return kFALSE;
af26ce80 653 }
43bfc8af 654
44dbae42 655 // Retrieve the CDB container class with the probability distributions
656 const AliTRDCalPID *pd = calibration->GetPIDObject(fPIDmethod == kNN ? 0 : 1);
657 if (!pd) {
658 AliError("No access to AliTRDCalPID");
659 return kFALSE;
660 }
43bfc8af 661
44dbae42 662 // Calculate the input for the NN if fPIDmethod is kNN
663 Float_t ldEdxNN[AliTRDCalPID::kNPlane * kNMLPslice], *dedx = 0x0;
664 if(fPIDmethod == kNN) {
665 CookdEdxNN(&ldEdxNN[0]);
666 }
43bfc8af 667
44dbae42 668 // Sets the a priori probabilities
669 for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
69dee759 670 fPID[ispec] = 1.0 / AliPID::kSPECIES;
44dbae42 671 }
43bfc8af 672
44dbae42 673 if(AliPID::kSPECIES != 5){
674 AliError("Probabilities array defined only for 5 values. Please modify !!");
675 return kFALSE;
676 }
43bfc8af 677
44dbae42 678 pidQuality = 0;
679 Float_t length; // track segment length in chamber
43bfc8af 680
44dbae42 681 // Skip tracks which have no TRD signal at all
682 if (fdEdx == 0.) return kFALSE;
683
684 for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
685
69dee759 686 length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())
687 / TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane])
688 / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
44dbae42 689
690 // check data
691 if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <= 0.
692 || fTimBinPlane[iPlane] == -1.) continue;
693
694 // this track segment has fulfilled all requierments
695 pidQuality++;
696
697 // Get the probabilities for the different particle species
698 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
699 switch(fPIDmethod){
700 case kLQ:
701 dedx = fdEdxPlane[iPlane];
702 break;
703 case kNN:
704 dedx = &ldEdxNN[iPlane*kNMLPslice];
705 break;
706 }
707 fPID[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], dedx, length, iPlane);
708 }
69dee759 709
710 }
711
712 if (pidQuality == 0) {
713 return kTRUE;
44dbae42 714 }
43bfc8af 715
44dbae42 716 // normalize probabilities
69dee759 717 Double_t probTotal = 0.0;
718 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
719 probTotal += fPID[iSpecies];
720 }
43bfc8af 721
69dee759 722 if (probTotal <= 0.0) {
723 AliWarning("The total probability over all species <= 0.");
724 AliWarning("This may be caused by some error in the reference histograms.");
44dbae42 725 return kFALSE;
726 }
69dee759 727
728 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
729 fPID[iSpecies] /= probTotal;
730 }
731
44dbae42 732 return kTRUE;
43bfc8af 733
69dee759 734}
44dbae42 735
a819a5f7 736//_____________________________________________________________________________
826fe01c 737Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
46d29e70 738{
fe33d239 739 //
826fe01c 740 // Propagates this track to a reference plane defined by "xk" [cm]
741 // correcting for the mean crossed material.
fe33d239 742 //
826fe01c 743 // "xx0" - thickness/rad.length [units of the radiation length]
744 // "xrho" - thickness*density [g/cm^2]
745 //
fe33d239 746
747 if (xk == GetX()) {
748 return kTRUE;
749 }
15ed8ba1 750
fe33d239 751 Double_t oldX = GetX();
752 Double_t oldY = GetY();
753 Double_t oldZ = GetZ();
15ed8ba1 754
fe33d239 755 Double_t bz = GetBz();
46d29e70 756
fe33d239 757 if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
758 return kFALSE;
759 }
9c9d2487 760
fe33d239 761 Double_t x = GetX();
762 Double_t y = GetY();
763 Double_t z = GetZ();
826fe01c 764
fe33d239 765 if (oldX < xk) {
826fe01c 766 xrho = -xrho;
fe33d239 767 if (IsStartedTimeIntegral()) {
69dee759 768 Double_t l2 = TMath::Sqrt((x-oldX)*(x-oldX)
769 + (y-oldY)*(y-oldY)
770 + (z-oldZ)*(z-oldZ));
fe33d239 771 Double_t crv = GetC();
772 if (TMath::Abs(l2*crv) > 0.0001) {
773 // Make correction for curvature if neccesary
69dee759 774 l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX)
775 + (y-oldY)*(y-oldY));
fe33d239 776 l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
777 l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
778 }
779 AddTimeStep(l2);
6c94f330 780 }
46d29e70 781 }
15ed8ba1 782
826fe01c 783 if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) {
fe33d239 784 return kFALSE;
785 }
15ed8ba1 786
fe33d239 787 {
788
69dee759 789 // Energy losses
6c23ffed 790 Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
fe33d239 791 Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
69dee759 792 if ((beta2 < 1.0e-10) ||
793 ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) {
fe33d239 794 return kFALSE;
795 }
796
797 Double_t dE = 0.153e-3 / beta2
69dee759 798 * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
826fe01c 799 * xrho;
800 fBudget[0] += xrho;
fe33d239 801
802 /*
803 // Suspicious part - think about it ?
804 Double_t kinE = TMath::Sqrt(p2);
805 if (dE > 0.8*kinE) dE = 0.8 * kinE; //
806 if (dE < 0) dE = 0.0; // Not valid region for Bethe bloch
807 */
808
809 fDE += dE;
810
811 /*
812 // Suspicious ! I.B.
813 Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE)); // Energy loss fluctuation
814 Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
815 fCcc += sigmac2;
816 fCee += fX*fX * sigmac2;
817 */
6c94f330 818
0d5b5c27 819 }
820
6c94f330 821 return kTRUE;
6d45eaef 822
46d29e70 823}
824
46d29e70 825//_____________________________________________________________________________
69dee759 826Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq
827 , Int_t index, Double_t h01)
46d29e70 828{
fe33d239 829 //
69dee759 830 // Assignes the found cluster <c> to the track and updates
831 // track information.
832 // <chisq> : predicted chi2
833 // <index> : ??
834 // <h01> : Tilting factor
fe33d239 835 //
46d29e70 836
69dee759 837 Double_t p[2] = { c->GetY()
838 , c->GetZ() };
fe33d239 839 Double_t sy2 = c->GetSigmaY2() * 4.0;
840 Double_t sz2 = c->GetSigmaZ2() * 4.0;
69dee759 841 Double_t cov[3] = { sy2 + h01*h01*sz2
842 , h01 * (sy2-sz2)
843 , sz2 + h01*h01*sy2 };
46e2d86c 844
fe33d239 845 if (!AliExternalTrackParam::Update(p,cov)) {
846 return kFALSE;
847 }
46d29e70 848
af26ce80 849 Int_t n = GetNumberOfClusters();
fe33d239 850 fIndex[n] = index;
b8dc2353 851 SetNumberOfClusters(n+1);
af26ce80 852
b8dc2353 853 SetChi2(GetChi2()+chisq);
5443e65e 854
af26ce80 855 return kTRUE;
fe33d239 856
857}
53c17fbf 858
46e2d86c 859//_____________________________________________________________________________
af26ce80 860Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index
69dee759 861 , Double_t h01, Int_t /*plane*/, Int_t /*tid*/)
fe33d239 862{
863 //
69dee759 864 // Assignes the found cluster <c> to the track and
865 // updates track information
866 // <chisq> : predicted chi2
867 // <index> : ??
868 // <h01> : Tilting factor
fe33d239 869 //
69dee759 870 // Difference to Update(AliTRDcluster *c): cluster is added to fClusters
6c94f330 871 //
15ed8ba1 872
69dee759 873 Double_t p[2] = { c->GetY()
874 , c->GetZ() };
fe33d239 875 Double_t sy2 = c->GetSigmaY2() * 4.0;
876 Double_t sz2 = c->GetSigmaZ2() * 4.0;
69dee759 877 Double_t cov[3] = { sy2 + h01*h01*sz2
878 , h01 * (sy2-sz2)
879 , sz2 + h01*h01*sy2 };
46e2d86c 880
fe33d239 881 if (!AliExternalTrackParam::Update(p,cov)) {
882 return kFALSE;
883 }
99910b5a 884
885 AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
fe33d239 886
af26ce80 887 // Register cluster to track
888 Int_t n = GetNumberOfClusters();
889 fIndex[n] = index;
890 fClusters[n] = c;
46e2d86c 891 SetNumberOfClusters(n+1);
fe33d239 892 SetChi2(GetChi2() + chisq);
893
894 return kTRUE;
46e2d86c 895
53c17fbf 896}
7ad19338 897
69dee759 898//_____________________________________________________________________________
899Bool_t AliTRDtrack::Update(const AliTRDtracklet &t, Double_t chisq, Int_t index)
900{
901 //
902 // Assignes the found tracklet <t> to the track
903 // and updates track information
904 // <chisq> : predicted chi2
905 // <index> : ??
906 //
7ad19338 907
69dee759 908 Double_t h01 = t.GetTilt(); // Is this correct????
909 Double_t p[2] = { t.GetY(), t.GetZ() };
910 Double_t sy2 = t.GetTrackletSigma2(); // Error pad-column
911 Double_t sz2 = 100000.0; // Error pad-row (????)
912 Double_t cov[3] = { sy2 + h01*h01*sz2 // Does this have any sense????
913 , h01 * (sy2 - sz2)
914 , sz2 + h01*h01*sy2 };
915 if (!AliExternalTrackParam::Update(p,cov)) {
916 return kFALSE;
917 }
6c94f330 918
69dee759 919 Int_t n = GetNumberOfClusters();
920 fIndex[n] = index;
921 SetNumberOfClusters(n+1);
922 SetChi2(GetChi2()+chisq);
7ad19338 923
69dee759 924 return kTRUE;
7ad19338 925
69dee759 926}
c84a5e9e 927
46d29e70 928//_____________________________________________________________________________
6c94f330 929Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
46d29e70 930{
fe33d239 931 //
6c94f330 932 // Rotates track parameters in R*phi plane
933 // if absolute rotation alpha is in global system
934 // otherwise alpha rotation is relative to the current rotation angle
fe33d239 935 //
936
3fad3d32 937 if (absolute) {
6c94f330 938 alpha -= GetAlpha();
3fad3d32 939 }
940 else{
941 fNRotate++;
942 }
46d29e70 943
6c94f330 944 return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
fe33d239 945
69dee759 946}
7ad19338 947
46d29e70 948//_____________________________________________________________________________
fd621f36 949Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
46d29e70 950{
53c17fbf 951 //
952 // Returns the track chi2
953 //
954
69dee759 955 Double_t p[2] = { c->GetY()
956 , c->GetZ() };
fe33d239 957 Double_t sy2 = c->GetSigmaY2() * 4.0;
958 Double_t sz2 = c->GetSigmaZ2() * 4.0;
69dee759 959 Double_t cov[3] = { sy2 + h01*h01*sz2
960 , h01*(sy2-sz2)
961 , sz2 + h01*h01*sy2 };
6c94f330 962
963 return AliExternalTrackParam::GetPredictedChi2(p,cov);
15ed8ba1 964
fe33d239 965}
7ad19338 966
53c17fbf 967//_____________________________________________________________________________
16d9fbba 968void AliTRDtrack::MakeBackupTrack()
969{
970 //
53c17fbf 971 // Creates a backup track
16d9fbba 972 //
53c17fbf 973
fe33d239 974 if (fBackupTrack) {
975 delete fBackupTrack;
976 }
16d9fbba 977 fBackupTrack = new AliTRDtrack(*this);
4f1c04d3 978
979}
980
53c17fbf 981//_____________________________________________________________________________
982Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
983{
4f1c04d3 984 //
fe33d239 985 // Find a prolongation at given x
986 // Return 0 if it does not exist
987 //
988
989 Double_t bz = GetBz();
15ed8ba1 990
fe33d239 991 if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
992 return 0;
993 }
994 if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
995 return 0;
996 }
4f1c04d3 997
6c94f330 998 return 1;
fe33d239 999
16d9fbba 1000}
3fad3d32 1001
53c17fbf 1002//_____________________________________________________________________________
fe33d239 1003Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
3fad3d32 1004{
1005 //
6c94f330 1006 // Propagate track to given x position
af26ce80 1007 // Works inside of the 20 degree segmentation
1008 // (local cooordinate frame for TRD , TPC, TOF)
3fad3d32 1009 //
fe33d239 1010 // Material budget from geo manager
6c94f330 1011 //
fe33d239 1012
1013 Double_t xyz0[3];
1014 Double_t xyz1[3];
1015 Double_t y;
1016 Double_t z;
1017
1018 const Double_t kAlphac = TMath::Pi()/9.0;
15ed8ba1 1019 const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
fe33d239 1020
1021 // Critical alpha - cross sector indication
af26ce80 1022 Double_t dir = (GetX() > xr) ? -1.0 : 1.0;
fe33d239 1023
1024 // Direction +-
1025 for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
1026
6c94f330 1027 GetXYZ(xyz0);
3fad3d32 1028 GetProlongation(x,y,z);
fe33d239 1029 xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
1030 xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
3fad3d32 1031 xyz1[2] = z;
1032 Double_t param[7];
826fe01c 1033 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 1034
1035 if ((param[0] > 0) &&
1036 (param[1] > 0)) {
826fe01c 1037 PropagateTo(x,param[1],param[0]*param[4]);
fe33d239 1038 }
1039
1040 if (GetY() > GetX()*kTalphac) {
53c17fbf 1041 Rotate(-kAlphac);
3fad3d32 1042 }
fe33d239 1043 if (GetY() < -GetX()*kTalphac) {
1044 Rotate( kAlphac);
3fad3d32 1045 }
fe33d239 1046
3fad3d32 1047 }
fe33d239 1048
3fad3d32 1049 PropagateTo(xr);
53c17fbf 1050
3fad3d32 1051 return 0;
3fad3d32 1052
53c17fbf 1053}
3fad3d32 1054
53c17fbf 1055//_____________________________________________________________________________
6c94f330 1056Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step)
3fad3d32 1057{
1058 //
fe33d239 1059 // Propagate track to the radial position
1060 // Rotation always connected to the last track position
3fad3d32 1061 //
fe33d239 1062
1063 Double_t xyz0[3];
1064 Double_t xyz1[3];
1065 Double_t y;
1066 Double_t z;
1067
6c94f330 1068 Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
fe33d239 1069 // Direction +-
af26ce80 1070 Double_t dir = (radius > r) ? -1.0 : 1.0;
fe33d239 1071
1072 for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
1073
6c94f330 1074 GetXYZ(xyz0);
3fad3d32 1075 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1076 Rotate(alpha,kTRUE);
6c94f330 1077 GetXYZ(xyz0);
3fad3d32 1078 GetProlongation(x,y,z);
fe33d239 1079 xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1080 xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
3fad3d32 1081 xyz1[2] = z;
1082 Double_t param[7];
826fe01c 1083 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 1084 if (param[1] <= 0) {
1085 param[1] = 100000000;
1086 }
826fe01c 1087 PropagateTo(x,param[1],param[0]*param[4]);
fe33d239 1088
3fad3d32 1089 }
fe33d239 1090
6c94f330 1091 GetXYZ(xyz0);
3fad3d32 1092 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1093 Rotate(alpha,kTRUE);
6c94f330 1094 GetXYZ(xyz0);
3fad3d32 1095 GetProlongation(r,y,z);
fe33d239 1096 xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1097 xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
3fad3d32 1098 xyz1[2] = z;
1099 Double_t param[7];
826fe01c 1100 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 1101
1102 if (param[1] <= 0) {
1103 param[1] = 100000000;
1104 }
826fe01c 1105 PropagateTo(r,param[1],param[0]*param[4]);
53c17fbf 1106
3fad3d32 1107 return 0;
53c17fbf 1108
1109}
1110
1111//_____________________________________________________________________________
4e28c495 1112Int_t AliTRDtrack::GetSector() const
53c17fbf 1113{
1114 //
1115 // Return the current sector
1116 //
1117
fe33d239 1118 return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
af26ce80 1119 % AliTRDgeometry::kNsect;
53c17fbf 1120
3fad3d32 1121}
1122
53c17fbf 1123//_____________________________________________________________________________
4e28c495 1124void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
53c17fbf 1125{
1126 //
1127 // The sampled energy loss
1128 //
1129
1130 Double_t s = GetSnp();
1131 Double_t t = GetTgl();
fe33d239 1132 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
53c17fbf 1133 fdQdl[i] = q;
1134
1135}
1136
fe33d239 1137//_____________________________________________________________________________
4e28c495 1138void AliTRDtrack::SetSampledEdx(Float_t q)
53c17fbf 1139{
1140 //
1141 // The sampled energy loss
1142 //
1143
1144 Double_t s = GetSnp();
1145 Double_t t = GetTgl();
fe33d239 1146 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
53c17fbf 1147 fdQdl[fNdedx] = q;
1148 fNdedx++;
1149
1150}
1151
fe33d239 1152//_____________________________________________________________________________
1153Double_t AliTRDtrack::GetBz() const
1154{
15ed8ba1 1155 //
fe33d239 1156 // Returns Bz component of the magnetic field (kG)
15ed8ba1 1157 //
53c17fbf 1158
fe33d239 1159 if (AliTracker::UniformField()) {
1160 return AliTracker::GetBz();
1161 }
1162 Double_t r[3];
1163 GetXYZ(r);
1164
1165 return AliTracker::GetBz(r);
53c17fbf 1166
fe33d239 1167}