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 "AliTRDReconstructor.h"
33 //#include "AliTRDcalibDB.h"
34 #include "AliTRDgeometry.h"
35 #include "AliTRDltuTracklet.h"
36 #include "AliTRDgtuTrack.h"
37 #include "Cal/AliTRDCalPID.h"
39 ClassImp(AliTRDgtuTrack)
41 //_____________________________________________________________________________
42 AliTRDgtuTrack::AliTRDgtuTrack()
44 ,fTracklets(new TObjArray(400))
60 // Default constructor
65 //_____________________________________________________________________________
66 AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack &t)
72 ,fDetector(t.fDetector)
73 ,fNtracklets(t.fNtracklets)
75 ,fNclusters(t.fNclusters)
81 ,fIsElectron(t.fIsElectron)
89 //_____________________________________________________________________________
90 AliTRDgtuTrack &AliTRDgtuTrack::operator=(const AliTRDgtuTrack &t)
93 // Assignment operator
96 if (this != &t) ((AliTRDgtuTrack &) t).Copy(*this);
102 //_____________________________________________________________________________
103 void AliTRDgtuTrack::Copy(TObject &t) const
109 ((AliTRDgtuTrack &) t).fTracklets = NULL;
110 ((AliTRDgtuTrack &) t).fYproj = fYproj;
111 ((AliTRDgtuTrack &) t).fZproj = fZproj;
112 ((AliTRDgtuTrack &) t).fSlope = fSlope;
113 ((AliTRDgtuTrack &) t).fDetector = fDetector;
114 ((AliTRDgtuTrack &) t).fNtracklets = fNtracklets;
115 ((AliTRDgtuTrack &) t).fNplanes = fNplanes;
116 ((AliTRDgtuTrack &) t).fNclusters = fNclusters;
117 ((AliTRDgtuTrack &) t).fPt = fPt;
118 ((AliTRDgtuTrack &) t).fPhi = fPhi;
119 ((AliTRDgtuTrack &) t).fEta = fEta;
120 ((AliTRDgtuTrack &) t).fLabel = fLabel;
121 ((AliTRDgtuTrack &) t).fPID = fPID;
122 ((AliTRDgtuTrack &) t).fIsElectron = fIsElectron;
126 //_____________________________________________________________________________
127 AliTRDgtuTrack::~AliTRDgtuTrack()
134 //fTracklets->Delete();
142 //_____________________________________________________________________________
143 void AliTRDgtuTrack::AddTracklet(AliTRDltuTracklet *trk)
146 // Add a LTU tracklet to this track
149 Tracklets()->Add(trk);
153 //_____________________________________________________________________________
154 AliTRDltuTracklet *AliTRDgtuTrack::GetTracklet(Int_t pos) const
157 // Return LTU tracklet at position "pos"
160 if (fTracklets == 0) {
163 void *trk = fTracklets->UncheckedAt(pos);
168 return (AliTRDltuTracklet *) trk;
172 //_____________________________________________________________________________
173 Int_t AliTRDgtuTrack::Compare(const TObject *o) const
176 // Compare function for sorting the tracks
179 AliTRDgtuTrack *gtutrack = (AliTRDgtuTrack *) o;
181 if (fYproj < gtutrack->GetYproj()) {
184 if (fYproj == gtutrack->GetYproj()) {
192 //_____________________________________________________________________________
193 void AliTRDgtuTrack::Reset()
196 // Reset the track information
211 fIsElectron = kFALSE;
215 //_____________________________________________________________________________
216 void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
219 // Calculate the kinematics of the found track
227 AliTRDltuTracklet *trk;
228 Int_t nTracklets = GetNtracklets();
229 Float_t fC[kNmaxTrk][3]; // X, Y, Z coordinates of segments
236 fNtracklets = GetNtracklets();
237 Int_t inDetector[kNlayer];
238 for (i = 0; i < kNlayer; i++) {
242 for (i = 0; i < nTracklets; i++) {
244 trk = GetTracklet(i);
245 fYproj += trk->GetYproj(xpl);
246 fZproj += trk->GetZproj(xpl);
247 fSlope += trk->GetSlope();
248 fNclusters += trk->GetNclusters();
249 iDet = trk->GetDetector();
252 for (Int_t id = 0; id < nDet; id++) {
253 if (iDet == inDetector[id]) {
254 newDetector = kFALSE;
259 inDetector[nDet++] = iDet;
263 fC[i][0] = trk->GetTime0();
264 fC[i][1] = trk->GetOffset();
265 fC[i][2] = trk->GetRowz();
268 fYproj /= (Float_t) nTracklets;
269 fZproj /= (Float_t) nTracklets;
270 fSlope /= (Float_t) nTracklets;
272 Float_t x[kNmaxTrk+1];
273 Float_t y[kNmaxTrk+1];
274 Float_t z[kNmaxTrk+1];
275 Bool_t count[kNmaxTrk];
276 for (i = 0; i < kNmaxTrk; i++) {
282 x[0] = y[0] = z[0] = 0.0;
283 while (j < nTracklets) {
285 for (i = 0; i < nTracklets; i++) {
286 if (count[i]) continue;
291 if (fC[i][0] < fC[iXmin][0]) {
295 x[j+1] = fC[iXmin][0];
296 y[j+1] = fC[iXmin][1];
297 z[j+1] = fC[iXmin][2];
299 count[iXmin] = kTRUE;
302 TMatrixD smatrix(2,2);
310 for (i = 0; i < nTracklets; i++) {
311 xv = (Double_t) x[i+1];
312 yv = (Double_t) y[i+1];
314 smatrix(1,1) += xv*xv;
320 res = smatrix.Invert() * sums;
324 Float_t dist = AliTRDgeometry::GetTime0(1) - AliTRDgeometry::GetTime0(0);
325 Float_t fx1 = x[1] + dist * (Float_t) (nTracklets-1) / 6.0;
326 Float_t fy1 = a + b * fx1;
327 Float_t fx2 = x[nTracklets] - dist * (Float_t) (nTracklets-1) / 6.0;
328 Float_t fy2 = a + b * fx2;
329 Float_t d12 = TMath::Sqrt((fx2-fx1)*(fx2-fx1)+(fy2-fy1)*(fy2-fy1));
330 Float_t alpha = TMath::ATan(fy2/fx2) - TMath::ATan(fy1/fx1);
331 Float_t r = (d12 / 2.0) / TMath::Sin(alpha);
333 fPt = 0.3 * field * 0.01 * r;
335 Float_t d1 = fx1*fx1 + fy1*fy1;
336 Float_t d2 = fx2*fx2 + fy2*fy2;
337 Float_t d = fx1*fy2 - fx2*fy1;
339 Float_t xc = (d1*fy2 - d2*fy1) / (2.0*d);
340 Float_t yc = (d2*fx1 - d1*fx2) / (2.0*d);
343 fPhi = TMath::ATan(xc/yc);
346 fPhi = TMath::PiOver2();
349 fPhi *= 180.0/TMath::Pi();
353 for (i = 0; i < nTracklets+1; i++) {
354 xv = (Double_t) z[i];
355 yv = (Double_t) x[i];
357 smatrix(1,1) += xv*xv;
363 res = smatrix.Invert() * sums;
366 Float_t theta = TMath::ATan(b);
369 theta = TMath::Pi() + theta;
375 fEta = -TMath::Log(TMath::Tan(theta/2.0));
380 //_____________________________________________________________________________
381 void AliTRDgtuTrack::MakePID()
384 // Electron likelihood signal
389 // Avoid dependency on libTRDrec.pkg (CBL)
390 //AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
391 //if (!calibration) {
392 // AliError("No instance of AliTRDcalibDB.");
396 //AliTRDrecoParam *rec = AliTRDReconstructor::RecoParam();
398 // AliError("No TRD reco param.");
401 //const AliTRDCalPID *pd = calibration->GetPIDObject(rec->GetPIDMethod());
403 AliTRDltuTracklet *trk;
404 Int_t nTracklets = GetNtracklets();
410 Float_t probPio = 1.0;
411 Float_t probEle = 1.0;
413 for (i = 0; i < nTracklets; i++) {
415 trk = GetTracklet(i);
417 sl = TMath::Abs(trk->GetSlope()); // Tracklet inclination in X-y plane
418 th = trk->GetRowz()/trk->GetTime0(); // Tracklet inclination in X-z plane
419 th = TMath::ATan(TMath::Abs(th));
422 * TMath::Cos(sl/180.0*TMath::Pi())
423 * TMath::Cos(th/180.0*TMath::Pi());
425 det = trk->GetDetector();
426 pla = trk->GetPlane(det);
428 // Unclear yet factor to match the PS distributions = 5.8
429 // not explained only by the tail filter ...
430 // AliRoot v4-03-07 , v4-03-Release
434 // Temporary (B. Vulpescu):
435 // The charge distributions do not match the new changes in simulation (A. Bercuci),
436 // which are nevertheless now in agreement with the beam tests.
437 // Some tricks will be used to still have reasonable results
438 // To match the existing charge distributions, the charge per layer has to be modified
449 // Since at tracking time we have no information on the particle type, we will modify
450 // instead the charge distributions accordingly. This will slow down the sampling.
451 // The modified distributions are in TRDdEdxHistogramsV1_BV.root and the CDB has
452 // been regenerated with AliTRDCreateDummyCDB.C
453 // The new PIDLQ data base has the version :
454 // I-AliCDBLocal::Get: CDB object retrieved: path "TRD/Calib/PIDLQ"; run range [0,0];
457 // HEAD28Mar06 new distributions: factor = 6.0
459 // A.Bercuci on 2nd May 2007
460 // Temporary modification to deal with more energy slices. The values
461 // attached to charge slices and track length are dummy
463 dedx[0] = dedx[1] = q*3.; dedx[2] = 0.;
464 //Float_t length = 3.7;
466 //probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length, 0);
467 //probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length, 0);
471 if ((probEle + probPio) > 0.0) {
472 fPID = probEle/(probEle + probPio);
478 // Thresholds for LQ cut at 90% electron efficiency (from AliRoot v4-03-07)
479 // P [GeV/c] fPIDcut (between 0 and 1)
488 // PIDcut = 0.978 - 0.024 * P[GeV/c]
489 //Float_t PIDcut = 0.978 - 0.024 * TMath::Abs(fPt);
491 // HEAD28Mar06 with modified distributions (A. Bercuci changes, P. Shukla distributions)
492 //Float_t fPIDcut = 0.829 - 0.032 * TMath::Abs(fPt);
494 // HEAD28Mar06 with new generated distributions (pol2 fit)
495 Float_t absPt = TMath::Abs(fPt);
496 //Float_t fPIDcut = 0.9575 - 0.0832 * absPt + 0.0037 * absPt*absPt; // 800 bins in dEdx
497 Float_t fPIDcut = 0.9482 - 0.0795 * absPt + 0.0035 * absPt*absPt; // 200 bins in dEdx
499 fIsElectron = kFALSE;
500 if (fPID >= fPIDcut) {
506 //_____________________________________________________________________________
507 void AliTRDgtuTrack::CookLabel()
510 // Cook the track label from tracklets labels
513 AliTRDltuTracklet *trk;
515 const Int_t kMaxTracks = 10;
516 Int_t trackLabel[kMaxTracks];
517 Int_t trackCount[kMaxTracks];
518 for (Int_t it = 0; it < kMaxTracks; it++) {
526 for (Int_t itrk = 0; itrk < fNtracklets; itrk++) {
528 trk = GetTracklet(itrk);
530 if (trk->GetLabel() == -1) {
534 label = trk->GetLabel();
537 for (Int_t it = 0; it < nTracks; it++) {
538 if (label == trackLabel[it]) {
545 trackLabel[nTracks] = label;
546 trackCount[nTracks]++;
548 if (nTracks == kMaxTracks) {
549 Warning("CookLabel","Too many tracks for this tracklet.");
557 Float_t frac = 4.0 / 5.0;
558 for (Int_t it = 0; it < kMaxTracks; it++) {
559 if (trackCount[it] >= (Int_t) (frac*fNtracklets)) {
560 fLabel = trackLabel[it];
567 //_____________________________________________________________________________
568 void AliTRDgtuTrack::ResetTracklets()
571 // Resets the list of tracklets
575 //fTracklets->Delete();
581 //_____________________________________________________________________________
582 TObjArray *AliTRDgtuTrack::Tracklets()
585 // Returns the list of tracklets
589 fTracklets = new TObjArray(400);
596 //_____________________________________________________________________________
597 Int_t AliTRDgtuTrack::GetNtracklets() const
600 // Returns the number of tracklets
604 return fTracklets->GetEntriesFast();