1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
21 // TRD module global track (GTU) //
24 ///////////////////////////////////////////////////////////////////////////////
28 #include <TObjArray.h>
32 #include "AliTRDgeometry.h"
33 #include "AliTRDcalibDB.h"
34 #include "AliTRDltuTracklet.h"
35 #include "AliTRDgtuTrack.h"
36 #include "Cal/AliTRDCalPID.h"
38 ClassImp(AliTRDgtuTrack)
40 //_____________________________________________________________________________
41 AliTRDgtuTrack::AliTRDgtuTrack()
43 ,fTracklets(new TObjArray(400))
59 // Default constructor
64 //_____________________________________________________________________________
65 AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack &t)
71 ,fDetector(t.fDetector)
72 ,fNtracklets(t.fNtracklets)
74 ,fNclusters(t.fNclusters)
80 ,fIsElectron(t.fIsElectron)
88 //_____________________________________________________________________________
89 AliTRDgtuTrack &AliTRDgtuTrack::operator=(const AliTRDgtuTrack &t)
92 // Assignment operator
95 if (this != &t) ((AliTRDgtuTrack &) t).Copy(*this);
101 //_____________________________________________________________________________
102 void AliTRDgtuTrack::Copy(TObject &t) const
108 ((AliTRDgtuTrack &) t).fTracklets = NULL;
109 ((AliTRDgtuTrack &) t).fYproj = fYproj;
110 ((AliTRDgtuTrack &) t).fZproj = fZproj;
111 ((AliTRDgtuTrack &) t).fSlope = fSlope;
112 ((AliTRDgtuTrack &) t).fDetector = fDetector;
113 ((AliTRDgtuTrack &) t).fNtracklets = fNtracklets;
114 ((AliTRDgtuTrack &) t).fNplanes = fNplanes;
115 ((AliTRDgtuTrack &) t).fNclusters = fNclusters;
116 ((AliTRDgtuTrack &) t).fPt = fPt;
117 ((AliTRDgtuTrack &) t).fPhi = fPhi;
118 ((AliTRDgtuTrack &) t).fEta = fEta;
119 ((AliTRDgtuTrack &) t).fLabel = fLabel;
120 ((AliTRDgtuTrack &) t).fPID = fPID;
121 ((AliTRDgtuTrack &) t).fIsElectron = fIsElectron;
125 //_____________________________________________________________________________
126 AliTRDgtuTrack::~AliTRDgtuTrack()
133 //fTracklets->Delete();
141 //_____________________________________________________________________________
142 void AliTRDgtuTrack::AddTracklet(AliTRDltuTracklet *trk)
145 // Add a LTU tracklet to this track
148 Tracklets()->Add(trk);
152 //_____________________________________________________________________________
153 AliTRDltuTracklet *AliTRDgtuTrack::GetTracklet(Int_t pos) const
156 // Return LTU tracklet at position "pos"
159 if (fTracklets == 0) {
162 void *trk = fTracklets->UncheckedAt(pos);
167 return (AliTRDltuTracklet *) trk;
171 //_____________________________________________________________________________
172 Int_t AliTRDgtuTrack::Compare(const TObject *o) const
175 // Compare function for sorting the tracks
178 AliTRDgtuTrack *gtutrack = (AliTRDgtuTrack *) o;
180 if (fYproj < gtutrack->GetYproj()) {
183 if (fYproj == gtutrack->GetYproj()) {
191 //_____________________________________________________________________________
192 void AliTRDgtuTrack::Reset()
195 // Reset the track information
210 fIsElectron = kFALSE;
214 //_____________________________________________________________________________
215 void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
218 // Calculate the kinematics of the found track
226 AliTRDltuTracklet *trk;
227 Int_t nTracklets = GetNtracklets();
228 Float_t fC[kNmaxTrk][3]; // X, Y, Z coordinates of segments
235 fNtracklets = GetNtracklets();
236 Int_t inDetector[kNplan];
237 for (i = 0; i < kNplan; i++) {
241 for (i = 0; i < nTracklets; i++) {
243 trk = GetTracklet(i);
244 fYproj += trk->GetYproj(xpl);
245 fZproj += trk->GetZproj(xpl);
246 fSlope += trk->GetSlope();
247 fNclusters += trk->GetNclusters();
248 iDet = trk->GetDetector();
251 for (Int_t id = 0; id < nDet; id++) {
252 if (iDet == inDetector[id]) {
253 newDetector = kFALSE;
258 inDetector[nDet++] = iDet;
262 fC[i][0] = trk->GetTime0();
263 fC[i][1] = trk->GetOffset();
264 fC[i][2] = trk->GetRowz();
267 fYproj /= (Float_t) nTracklets;
268 fZproj /= (Float_t) nTracklets;
269 fSlope /= (Float_t) nTracklets;
271 Float_t x[kNmaxTrk+1];
272 Float_t y[kNmaxTrk+1];
273 Float_t z[kNmaxTrk+1];
274 Bool_t count[kNmaxTrk];
275 for (i = 0; i < kNmaxTrk; i++) {
281 x[0] = y[0] = z[0] = 0.0;
282 while (j < nTracklets) {
284 for (i = 0; i < nTracklets; i++) {
285 if (count[i]) continue;
290 if (fC[i][0] < fC[iXmin][0]) {
294 x[j+1] = fC[iXmin][0];
295 y[j+1] = fC[iXmin][1];
296 z[j+1] = fC[iXmin][2];
298 count[iXmin] = kTRUE;
301 TMatrixD smatrix(2,2);
309 for (i = 0; i < nTracklets; i++) {
310 xv = (Double_t) x[i+1];
311 yv = (Double_t) y[i+1];
313 smatrix(1,1) += xv*xv;
319 res = smatrix.Invert() * sums;
323 Float_t dist = AliTRDgeometry::GetTime0(1) - AliTRDgeometry::GetTime0(0);
324 Float_t fx1 = x[1] + dist * (Float_t) (nTracklets-1) / 6.0;
325 Float_t fy1 = a + b * fx1;
326 Float_t fx2 = x[nTracklets] - dist * (Float_t) (nTracklets-1) / 6.0;
327 Float_t fy2 = a + b * fx2;
328 Float_t d12 = TMath::Sqrt((fx2-fx1)*(fx2-fx1)+(fy2-fy1)*(fy2-fy1));
329 Float_t alpha = TMath::ATan(fy2/fx2) - TMath::ATan(fy1/fx1);
330 Float_t r = (d12 / 2.0) / TMath::Sin(alpha);
332 fPt = 0.3 * field * 0.01 * r;
334 Float_t d1 = fx1*fx1 + fy1*fy1;
335 Float_t d2 = fx2*fx2 + fy2*fy2;
336 Float_t d = fx1*fy2 - fx2*fy1;
338 Float_t xc = (d1*fy2 - d2*fy1) / (2.0*d);
339 Float_t yc = (d2*fx1 - d1*fx2) / (2.0*d);
342 fPhi = TMath::ATan(xc/yc);
345 fPhi = TMath::PiOver2();
348 fPhi *= 180.0/TMath::Pi();
352 for (i = 0; i < nTracklets+1; i++) {
353 xv = (Double_t) z[i];
354 yv = (Double_t) x[i];
356 smatrix(1,1) += xv*xv;
362 res = smatrix.Invert() * sums;
365 Float_t theta = TMath::ATan(b);
368 theta = TMath::Pi() + theta;
374 fEta = -TMath::Log(TMath::Tan(theta/2.0));
379 //_____________________________________________________________________________
380 void AliTRDgtuTrack::MakePID()
383 // Electron likelihood signal
388 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
390 AliError("No instance of AliTRDcalibDB.");
393 const AliTRDCalPID *pd = calibration->GetPIDLQObject();
395 AliTRDltuTracklet *trk;
396 Int_t nTracklets = GetNtracklets();
402 Float_t probPio = 1.0;
403 Float_t probEle = 1.0;
405 for (i = 0; i < nTracklets; i++) {
407 trk = GetTracklet(i);
409 sl = TMath::Abs(trk->GetSlope()); // Tracklet inclination in X-y plane
410 th = trk->GetRowz()/trk->GetTime0(); // Tracklet inclination in X-z plane
411 th = TMath::ATan(TMath::Abs(th));
414 * TMath::Cos(sl/180.0*TMath::Pi())
415 * TMath::Cos(th/180.0*TMath::Pi());
417 det = trk->GetDetector();
418 pla = trk->GetPlane(det);
420 // Unclear yet factor to match the PS distributions = 5.8
421 // not explained only by the tail filter ...
422 // AliRoot v4-03-07 , v4-03-Release
426 // Temporary (B. Vulpescu):
427 // The charge distributions do not match the new changes in simulation (A. Bercuci),
428 // which are nevertheless now in agreement with the beam tests.
429 // Some tricks will be used to still have reasonable results
430 // To match the existing charge distributions, the charge per layer has to be modified
441 // Since at tracking time we have no information on the particle type, we will modify
442 // instead the charge distributions accordingly. This will slow down the sampling.
443 // The modified distributions are in TRDdEdxHistogramsV1_BV.root and the CDB has
444 // been regenerated with AliTRDCreateDummyCDB.C
445 // The new PIDLQ data base has the version :
446 // I-AliCDBLocal::Get: CDB object retrieved: path "TRD/Calib/PIDLQ"; run range [0,0];
449 // HEAD28Mar06 new distributions: factor = 6.0
451 // A.Bercuci on 2nd May 2007
452 // Temporary modification to deal with more energy slices. The values
453 // attached to charge slices and track length are dummy
455 dedx[0] = dedx[1] = q*3.; dedx[2] = 0.;
456 Float_t length = 3.7;
458 probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length);
459 probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length);
463 if ((probEle + probPio) > 0.0) {
464 fPID = probEle/(probEle + probPio);
470 // Thresholds for LQ cut at 90% electron efficiency (from AliRoot v4-03-07)
471 // P [GeV/c] fPIDcut (between 0 and 1)
480 // PIDcut = 0.978 - 0.024 * P[GeV/c]
481 //Float_t PIDcut = 0.978 - 0.024 * TMath::Abs(fPt);
483 // HEAD28Mar06 with modified distributions (A. Bercuci changes, P. Shukla distributions)
484 //Float_t fPIDcut = 0.829 - 0.032 * TMath::Abs(fPt);
486 // HEAD28Mar06 with new generated distributions (pol2 fit)
487 Float_t absPt = TMath::Abs(fPt);
488 //Float_t fPIDcut = 0.9575 - 0.0832 * absPt + 0.0037 * absPt*absPt; // 800 bins in dEdx
489 Float_t fPIDcut = 0.9482 - 0.0795 * absPt + 0.0035 * absPt*absPt; // 200 bins in dEdx
491 fIsElectron = kFALSE;
492 if (fPID >= fPIDcut) {
498 //_____________________________________________________________________________
499 void AliTRDgtuTrack::CookLabel()
502 // Cook the track label from tracklets labels
505 AliTRDltuTracklet *trk;
507 const Int_t kMaxTracks = 10;
508 Int_t trackLabel[kMaxTracks];
509 Int_t trackCount[kMaxTracks];
510 for (Int_t it = 0; it < kMaxTracks; it++) {
518 for (Int_t itrk = 0; itrk < fNtracklets; itrk++) {
520 trk = GetTracklet(itrk);
522 if (trk->GetLabel() == -1) {
526 label = trk->GetLabel();
529 for (Int_t it = 0; it < nTracks; it++) {
530 if (label == trackLabel[it]) {
537 trackLabel[nTracks] = label;
538 trackCount[nTracks]++;
540 if (nTracks == kMaxTracks) {
541 Warning("CookLabel","Too many tracks for this tracklet.");
549 Float_t frac = 4.0 / 5.0;
550 for (Int_t it = 0; it < kMaxTracks; it++) {
551 if (trackCount[it] >= (Int_t) (frac*fNtracklets)) {
552 fLabel = trackLabel[it];
559 //_____________________________________________________________________________
560 void AliTRDgtuTrack::ResetTracklets()
563 // Resets the list of tracklets
567 //fTracklets->Delete();
573 //_____________________________________________________________________________
574 TObjArray *AliTRDgtuTrack::Tracklets()
577 // Returns the list of tracklets
581 fTracklets = new TObjArray(400);
588 //_____________________________________________________________________________
589 Int_t AliTRDgtuTrack::GetNtracklets() const
592 // Returns the number of tracklets
596 return fTracklets->GetEntriesFast();