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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////////
19 // TRD module global track (GTU) //
22 ///////////////////////////////////////////////////////////////////////////////
25 #include <TObjArray.h>
29 #include "AliTRDgeometry.h"
30 #include "AliTRDcalibDB.h"
31 #include "AliTRDltuTracklet.h"
32 #include "AliTRDgtuTrack.h"
33 #include "Cal/AliTRDCalPIDLQ.h"
35 ClassImp(AliTRDgtuTrack)
37 //_____________________________________________________________________________
38 AliTRDgtuTrack::AliTRDgtuTrack()
40 ,fTracklets(new TObjArray(400))
56 // Default constructor
61 //_____________________________________________________________________________
62 AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack &t)
68 ,fDetector(t.fDetector)
69 ,fNtracklets(t.fNtracklets)
71 ,fNclusters(t.fNclusters)
77 ,fIsElectron(t.fIsElectron)
85 //_____________________________________________________________________________
86 AliTRDgtuTrack &AliTRDgtuTrack::operator=(const AliTRDgtuTrack &t)
89 // Assignment operator
92 if (this != &t) ((AliTRDgtuTrack &) t).Copy(*this);
98 //_____________________________________________________________________________
99 void AliTRDgtuTrack::Copy(TObject &t) const
105 ((AliTRDgtuTrack &) t).fTracklets = NULL;
106 ((AliTRDgtuTrack &) t).fYproj = fYproj;
107 ((AliTRDgtuTrack &) t).fZproj = fZproj;
108 ((AliTRDgtuTrack &) t).fSlope = fSlope;
109 ((AliTRDgtuTrack &) t).fDetector = fDetector;
110 ((AliTRDgtuTrack &) t).fNtracklets = fNtracklets;
111 ((AliTRDgtuTrack &) t).fNplanes = fNplanes;
112 ((AliTRDgtuTrack &) t).fNclusters = fNclusters;
113 ((AliTRDgtuTrack &) t).fPt = fPt;
114 ((AliTRDgtuTrack &) t).fPhi = fPhi;
115 ((AliTRDgtuTrack &) t).fEta = fEta;
116 ((AliTRDgtuTrack &) t).fLabel = fLabel;
117 ((AliTRDgtuTrack &) t).fPID = fPID;
118 ((AliTRDgtuTrack &) t).fIsElectron = fIsElectron;
122 //_____________________________________________________________________________
123 AliTRDgtuTrack::~AliTRDgtuTrack()
130 //fTracklets->Delete();
138 //_____________________________________________________________________________
139 void AliTRDgtuTrack::AddTracklet(AliTRDltuTracklet *trk)
142 // Add a LTU tracklet to this track
145 Tracklets()->Add(trk);
149 //_____________________________________________________________________________
150 AliTRDltuTracklet *AliTRDgtuTrack::GetTracklet(Int_t pos) const
153 // Return LTU tracklet at position "pos"
156 if (fTracklets == 0) {
159 void *trk = fTracklets->UncheckedAt(pos);
164 return (AliTRDltuTracklet *) trk;
168 //_____________________________________________________________________________
169 Int_t AliTRDgtuTrack::Compare(const TObject *o) const
172 // Compare function for sorting the tracks
175 AliTRDgtuTrack *gtutrack = (AliTRDgtuTrack *) o;
177 if (fYproj < gtutrack->GetYproj()) {
180 if (fYproj == gtutrack->GetYproj()) {
188 //_____________________________________________________________________________
189 void AliTRDgtuTrack::Reset()
192 // Reset the track information
207 fIsElectron = kFALSE;
211 //_____________________________________________________________________________
212 void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
215 // Calculate the kinematics of the found track
223 AliTRDltuTracklet *trk;
224 Int_t nTracklets = GetNtracklets();
225 Float_t fC[kNmaxTrk][3]; // X, Y, Z coordinates of segments
232 fNtracklets = GetNtracklets();
233 Int_t inDetector[kNplan];
234 for (i = 0; i < kNplan; i++) {
238 for (i = 0; i < nTracklets; i++) {
240 trk = GetTracklet(i);
241 fYproj += trk->GetYproj(xpl);
242 fZproj += trk->GetZproj(xpl);
243 fSlope += trk->GetSlope();
244 fNclusters += trk->GetNclusters();
245 iDet = trk->GetDetector();
248 for (Int_t id = 0; id < nDet; id++) {
249 if (iDet == inDetector[id]) {
250 newDetector = kFALSE;
255 inDetector[nDet++] = iDet;
259 fC[i][0] = trk->GetTime0();
260 fC[i][1] = trk->GetOffset();
261 fC[i][2] = trk->GetRowz();
264 fYproj /= (Float_t) nTracklets;
265 fZproj /= (Float_t) nTracklets;
266 fSlope /= (Float_t) nTracklets;
268 Float_t x[kNmaxTrk+1];
269 Float_t y[kNmaxTrk+1];
270 Float_t z[kNmaxTrk+1];
271 Bool_t count[kNmaxTrk];
272 for (Int_t i = 0; i < kNmaxTrk; i++) {
278 x[0] = y[0] = z[0] = 0.0;
279 while (j < nTracklets) {
281 for (Int_t i = 0; i < nTracklets; i++) {
282 if (count[i]) continue;
287 if (fC[i][0] < fC[iXmin][0]) {
291 x[j+1] = fC[iXmin][0];
292 y[j+1] = fC[iXmin][1];
293 z[j+1] = fC[iXmin][2];
295 count[iXmin] = kTRUE;
298 TMatrixD smatrix(2,2);
306 for (i = 0; i < nTracklets; i++) {
307 xv = (Double_t) x[i+1];
308 yv = (Double_t) y[i+1];
310 smatrix(1,1) += xv*xv;
316 res = smatrix.Invert() * sums;
320 Float_t dist = AliTRDgeometry::GetTime0(1) - AliTRDgeometry::GetTime0(0);
321 Float_t fx1 = x[1] + dist * (Float_t) (nTracklets-1) / 6.0;
322 Float_t fy1 = a + b * fx1;
323 Float_t fx2 = x[nTracklets] - dist * (Float_t) (nTracklets-1) / 6.0;
324 Float_t fy2 = a + b * fx2;
325 Float_t d12 = TMath::Sqrt((fx2-fx1)*(fx2-fx1)+(fy2-fy1)*(fy2-fy1));
326 Float_t alpha = TMath::ATan(fy2/fx2) - TMath::ATan(fy1/fx1);
327 Float_t r = (d12 / 2.0) / TMath::Sin(alpha);
329 fPt = 0.3 * field * 0.01 * r;
331 Float_t d1 = fx1*fx1 + fy1*fy1;
332 Float_t d2 = fx2*fx2 + fy2*fy2;
333 Float_t d = fx1*fy2 - fx2*fy1;
335 Float_t xc = (d1*fy2 - d2*fy1) / (2.0*d);
336 Float_t yc = (d2*fx1 - d1*fx2) / (2.0*d);
339 fPhi = TMath::ATan(xc/yc);
342 fPhi = TMath::PiOver2();
345 fPhi *= 180.0/TMath::Pi();
349 for (i = 0; i < nTracklets+1; i++) {
350 xv = (Double_t) z[i];
351 yv = (Double_t) x[i];
353 smatrix(1,1) += xv*xv;
359 res = smatrix.Invert() * sums;
362 Float_t theta = TMath::ATan(b);
365 theta = TMath::Pi() + theta;
371 fEta = -TMath::Log(TMath::Tan(theta/2.0));
376 //_____________________________________________________________________________
377 void AliTRDgtuTrack::MakePID()
380 // Electron likelihood signal
385 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
387 AliError("No instance of AliTRDcalibDB.");
390 const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
392 AliTRDltuTracklet *trk;
393 Int_t nTracklets = GetNtracklets();
399 Float_t probPio = 1.0;
400 Float_t probEle = 1.0;
402 for (i = 0; i < nTracklets; i++) {
404 trk = GetTracklet(i);
406 sl = TMath::Abs(trk->GetSlope()); // Tracklet inclination in X-y plane
407 th = trk->GetRowz()/trk->GetTime0(); // Tracklet inclination in X-z plane
408 th = TMath::ATan(TMath::Abs(th));
411 * TMath::Cos(sl/180.0*TMath::Pi())
412 * TMath::Cos(th/180.0*TMath::Pi());
414 det = trk->GetDetector();
415 pla = trk->GetPlane(det);
417 // Unclear yet factor to match the PS distributions = 5.8
418 // not explained only by the tail filter ...
419 // AliRoot v4-03-07 , v4-03-Release
423 // Temporary (B. Vulpescu):
424 // The charge distributions do not match the new changes in simulation (A. Bercuci),
425 // which are nevertheless now in agreement with the beam tests.
426 // Some tricks will be used to still have reasonable results
427 // To match the existing charge distributions, the charge per layer has to be modified
438 // Since at tracking time we have no information on the particle type, we will modify
439 // instead the charge distributions accordingly. This will slow down the sampling.
440 // The modified distributions are in TRDdEdxHistogramsV1_BV.root and the CDB has
441 // been regenerated with AliTRDCreateDummyCDB.C
442 // The new PIDLQ data base has the version :
443 // I-AliCDBLocal::Get: CDB object retrieved: path "TRD/Calib/PIDLQ"; run range [0,0];
446 // HEAD28Mar06 new distributions: factor = 6.0
448 probEle *= pd->GetProbability(0,TMath::Abs(fPt),q*6.0);
449 probPio *= pd->GetProbability(2,TMath::Abs(fPt),q*6.0);
453 if ((probEle + probPio) > 0.0) {
454 fPID = probEle/(probEle + probPio);
460 // Thresholds for LQ cut at 90% electron efficiency (from AliRoot v4-03-07)
461 // P [GeV/c] fPIDcut (between 0 and 1)
470 // PIDcut = 0.978 - 0.024 * P[GeV/c]
471 //Float_t PIDcut = 0.978 - 0.024 * TMath::Abs(fPt);
473 // HEAD28Mar06 with modified distributions (A. Bercuci changes, P. Shukla distributions)
474 //Float_t fPIDcut = 0.829 - 0.032 * TMath::Abs(fPt);
476 // HEAD28Mar06 with new generated distributions (pol2 fit)
477 Float_t absPt = TMath::Abs(fPt);
478 //Float_t fPIDcut = 0.9575 - 0.0832 * absPt + 0.0037 * absPt*absPt; // 800 bins in dEdx
479 Float_t fPIDcut = 0.9482 - 0.0795 * absPt + 0.0035 * absPt*absPt; // 200 bins in dEdx
481 fIsElectron = kFALSE;
482 if (fPID >= fPIDcut) {
488 //_____________________________________________________________________________
489 void AliTRDgtuTrack::CookLabel()
492 // Cook the track label from tracklets labels
495 AliTRDltuTracklet *trk;
497 const Int_t kMaxTracks = 10;
498 Int_t trackLabel[kMaxTracks];
499 Int_t trackCount[kMaxTracks];
500 for (Int_t it = 0; it < kMaxTracks; it++) {
508 for (Int_t itrk = 0; itrk < fNtracklets; itrk++) {
510 trk = GetTracklet(itrk);
512 if (trk->GetLabel() == -1) {
516 label = trk->GetLabel();
519 for (Int_t it = 0; it < nTracks; it++) {
520 if (label == trackLabel[it]) {
527 trackLabel[nTracks] = label;
528 trackCount[nTracks]++;
530 if (nTracks == kMaxTracks) {
531 Warning("CookLabel","Too many tracks for this tracklet.");
539 Float_t frac = 4.0 / 5.0;
540 for (Int_t it = 0; it < kMaxTracks; it++) {
541 if (trackCount[it] >= (Int_t) (frac*fNtracklets)) {
542 fLabel = trackLabel[it];
549 //_____________________________________________________________________________
550 void AliTRDgtuTrack::ResetTracklets()
553 // Resets the list of tracklets
557 //fTracklets->Delete();
563 //_____________________________________________________________________________
564 TObjArray *AliTRDgtuTrack::Tracklets()
567 // Returns the list of tracklets
571 fTracklets = new TObjArray(400);
578 //_____________________________________________________________________________
579 Int_t AliTRDgtuTrack::GetNtracklets() const
582 // Returns the number of tracklets
586 return fTracklets->GetEntriesFast();