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>
27 #include "AliTRDgeometry.h"
28 #include "AliTRDcalibDB.h"
29 #include "Cal/AliTRDCalPIDLQ.h"
31 #include "AliTRDltuTracklet.h"
33 #include "AliTRDgtuTrack.h"
35 ClassImp(AliTRDgtuTrack)
37 //_____________________________________________________________________________
38 AliTRDgtuTrack::AliTRDgtuTrack()
41 // Default constructor
58 fTracklets = new TObjArray(400);
62 //_____________________________________________________________________________
63 AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack& track):
69 fDetector(track.fDetector),
70 fNtracklets(track.fNtracklets),
71 fNplanes(track.fNplanes),
72 fNclusters(track.fNclusters),
78 fIsElectron(track.fIsElectron)
86 //_____________________________________________________________________________
87 AliTRDgtuTrack &AliTRDgtuTrack::operator=(const AliTRDgtuTrack &t)
90 // assignment operator
93 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()
131 //_____________________________________________________________________________
132 void AliTRDgtuTrack::AddTracklet(AliTRDltuTracklet *trk)
135 // Add a LTU tracklet to this track
138 Tracklets()->Add(trk);
142 //_____________________________________________________________________________
143 AliTRDltuTracklet* AliTRDgtuTrack::GetTracklet(Int_t pos) const
146 // Return LTU tracklet at position "pos"
149 if (fTracklets == 0) return 0;
150 void * trk = fTracklets->UncheckedAt(pos);
151 if (trk == 0) return 0;
153 return (AliTRDltuTracklet*)trk;
157 //_____________________________________________________________________________
158 Int_t AliTRDgtuTrack::Compare(const TObject * o) const
161 // Compare function for sorting the tracks
164 AliTRDgtuTrack *gtutrack = (AliTRDgtuTrack*)o;
166 if (fYproj < gtutrack->GetYproj()) return -1;
167 if (fYproj == gtutrack->GetYproj()) return 0;
173 //_____________________________________________________________________________
174 void AliTRDgtuTrack::Reset()
177 // Reset the track information
192 fIsElectron = kFALSE;
196 //_____________________________________________________________________________
197 void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
200 // Calculate the kinematics of the found track
203 AliTRDltuTracklet *trk;
204 Int_t nTracklets = GetNtracklets();
205 Float_t fC[kNmaxTrk][3]; // X, Y, Z coordinates of segments
212 fNtracklets = GetNtracklets();
213 Int_t inDetector[kNplan];
214 for (Int_t i = 0; i < kNplan; i++) inDetector[i] = -1;
215 Int_t iDet, nDet = 0;
217 for (Int_t i = 0; i < nTracklets; i++) {
219 trk = GetTracklet(i);
220 fYproj += trk->GetYproj(xpl);
221 fZproj += trk->GetZproj(xpl);
222 fSlope += trk->GetSlope();
223 fNclusters += trk->GetNclusters();
224 iDet = trk->GetDetector();
227 for (Int_t id = 0; id < nDet; id++) {
228 if (iDet == inDetector[id]) {
229 newDetector = kFALSE;
234 inDetector[nDet++] = iDet;
238 fC[i][0] = trk->GetTime0();
239 fC[i][1] = trk->GetOffset();
240 fC[i][2] = trk->GetRowz();
243 fYproj /= (Float_t)nTracklets;
244 fZproj /= (Float_t)nTracklets;
245 fSlope /= (Float_t)nTracklets;
247 Float_t x[kNmaxTrk+1], y[kNmaxTrk+1], z[kNmaxTrk+1];
248 Bool_t count[kNmaxTrk];
249 for (Int_t i = 0; i < kNmaxTrk; i++) count[i] = kFALSE;
253 x[0] = y[0] = z[0] = 0.0;
254 while (j < nTracklets) {
256 for (Int_t i = 0; i < nTracklets; i++) {
257 if (count[i]) continue;
262 if (fC[i][0] < fC[iXmin][0]) iXmin = i;
264 x[j+1] = fC[iXmin][0];
265 y[j+1] = fC[iXmin][1];
266 z[j+1] = fC[iXmin][2];
268 count[iXmin] = kTRUE;
271 TMatrixD smatrix(2,2);
279 for (Int_t i = 0; i < nTracklets; i++) {
280 xv = (Double_t)x[i+1];
281 yv = (Double_t)y[i+1];
283 smatrix(1,1) += xv*xv;
289 res = smatrix.Invert() * sums;
293 Float_t dist = AliTRDgeometry::GetTime0(1) - AliTRDgeometry::GetTime0(0);
295 Float_t fx1 = x[1] + dist * (Float_t)(nTracklets-1)/6.0;
296 Float_t fy1 = a + b * fx1;
297 Float_t fx2 = x[nTracklets] - dist * (Float_t)(nTracklets-1)/6.0;
298 Float_t fy2 = a + b * fx2;
299 Float_t d12 = TMath::Sqrt((fx2-fx1)*(fx2-fx1)+(fy2-fy1)*(fy2-fy1));
300 Float_t alpha = TMath::ATan(fy2/fx2) - TMath::ATan(fy1/fx1);
301 Float_t r = (d12/2.0)/TMath::Sin(alpha);
303 fPt = 0.3 * field * 0.01 * r;
305 Float_t d1 = fx1*fx1+fy1*fy1;
306 Float_t d2 = fx2*fx2+fy2*fy2;
307 Float_t d = fx1*fy2-fx2*fy1;
309 Float_t xc = (d1*fy2-d2*fy1)/(2*d);
310 Float_t yc = (d2*fx1-d1*fx2)/(2*d);
313 fPhi = TMath::ATan(xc/yc);
315 fPhi = TMath::PiOver2();
318 fPhi *= 180.0/TMath::Pi();
322 for (Int_t i = 0; i < nTracklets+1; i++) {
326 smatrix(1,1) += xv*xv;
332 res = smatrix.Invert() * sums;
335 Float_t theta = TMath::ATan(b);
337 if (theta < 0.0) theta = TMath::Pi() + theta;
342 fEta = -TMath::Log(TMath::Tan(theta/2.0));
347 //_____________________________________________________________________________
348 void AliTRDgtuTrack::MakePID()
351 // Electron likelihood signal
354 AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
357 Error("MakePID","No instance of AliTRDcalibDB.");
360 const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
362 AliTRDltuTracklet *trk;
363 Int_t nTracklets = GetNtracklets();
365 Float_t sl, th, q, probPio = 1.0, probEle = 1.0;
366 for (Int_t i = 0; i < nTracklets; i++) {
368 trk = GetTracklet(i);
370 sl = TMath::Abs(trk->GetSlope()); // tracklet inclination in X-y plane
371 th = trk->GetRowz()/trk->GetTime0(); // tracklet inclination in X-z plane
372 th = TMath::ATan(TMath::Abs(th));
375 * TMath::Cos(sl/180.0*TMath::Pi())
376 * TMath::Cos(th/180.0*TMath::Pi());
378 det = trk->GetDetector();
379 pla = trk->GetPlane(det);
381 // unclear yet factor to match the PS distributions = 5.8
382 // not explained only by the tail filter ...
384 // AliRoot v4-03-07 , v4-03-Release
388 // Temporary (B. Vulpescu):
389 // The charge distributions do not match the new changes in simulation (A. Bercuci),
390 // which are nevertheless now in agreement with the beam tests.
391 // Some tricks will be used to still have reasonable results
392 // To match the existing charge distributions, the charge per layer has to be modified
403 // Since at tracking time we have no information on the particle type, we will modify
404 // instead the charge distributions accordingly. This will slow down the sampling.
405 // The modified distributions are in TRDdEdxHistogramsV1_BV.root and the CDB has
406 // been regenerated with AliTRDCreateDummyCDB.C
407 // The new PIDLQ data base has the version :
408 // I-AliCDBLocal::Get: CDB object retrieved: path "TRD/Calib/PIDLQ"; run range [0,0];
411 // HEAD28Mar06 new distributions: factor = 6.0
413 probEle *= pd->GetProbability(0,TMath::Abs(fPt),q*6.0);
414 probPio *= pd->GetProbability(2,TMath::Abs(fPt),q*6.0);
418 if ((probEle+probPio) > 0.0) {
419 fPID = probEle/(probEle+probPio);
424 // Thresholds for LQ cut at 90% electron efficiency (from AliRoot v4-03-07)
425 // P [GeV/c] fPIDcut (between 0 and 1)
434 // PIDcut = 0.978 - 0.024 * P[GeV/c]
435 //Float_t PIDcut = 0.978 - 0.024 * TMath::Abs(fPt);
437 // HEAD28Mar06 with modified distributions (A. Bercuci changes, P. Shukla distributions)
438 //Float_t fPIDcut = 0.829 - 0.032 * TMath::Abs(fPt);
440 // HEAD28Mar06 with new generated distributions (pol2 fit)
441 Float_t absPt = TMath::Abs(fPt);
442 //Float_t fPIDcut = 0.9575 - 0.0832 * absPt + 0.0037 * absPt*absPt; // 800 bins in dEdx
443 Float_t fPIDcut = 0.9482 - 0.0795 * absPt + 0.0035 * absPt*absPt; // 200 bins in dEdx
445 fIsElectron = kFALSE;
446 if (fPID >= fPIDcut) {
452 //_____________________________________________________________________________
453 void AliTRDgtuTrack::CookLabel()
456 // Cook the track label from tracklets labels
459 AliTRDltuTracklet *trk;
461 const Int_t kMaxTracks = 10;
462 Int_t trackLabel[kMaxTracks];
463 Int_t trackCount[kMaxTracks];
464 for (Int_t it = 0; it < kMaxTracks; it++) {
470 Int_t label, nTracks = 0;
471 for (Int_t itrk = 0; itrk < fNtracklets; itrk++) {
473 trk = GetTracklet(itrk);
475 if (trk->GetLabel() == -1) continue;
477 label = trk->GetLabel();
480 for (Int_t it = 0; it < nTracks; it++) {
481 if (label == trackLabel[it]) {
488 trackLabel[nTracks] = label;
489 trackCount[nTracks]++;
491 if (nTracks == kMaxTracks) {
492 Warning("CookLabel","Too many tracks for this tracklet.");
500 Float_t frac = 4.0/5.0;
501 for (Int_t it = 0; it < kMaxTracks; it++) {
502 if (trackCount[it] >= (Int_t)(frac*fNtracklets)) {
503 fLabel = trackLabel[it];