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