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