First implementation of neural network PID
[u/mrichter/AliRoot.git] / TRD / AliTRDgtuTrack.cxx
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  *                                                                        *
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  **************************************************************************/
15
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //                                                                           //
21 //  TRD module global track (GTU)                                            //
22 //                                                                           //
23 //                                                                           //
24 ///////////////////////////////////////////////////////////////////////////////
25
26 #include <TMath.h>
27 #include <TMatrixD.h>
28 #include <TObjArray.h>
29
30 #include "AliLog.h"
31
32 #include "AliTRDgeometry.h"
33 #include "AliTRDcalibDB.h"
34 #include "AliTRDltuTracklet.h"
35 #include "AliTRDgtuTrack.h"
36 #include "Cal/AliTRDCalPID.h"
37
38 ClassImp(AliTRDgtuTrack)
39
40 //_____________________________________________________________________________
41 AliTRDgtuTrack::AliTRDgtuTrack()
42   :TObject()
43   ,fTracklets(new TObjArray(400))
44   ,fYproj(0)
45   ,fZproj(0)
46   ,fSlope(0)
47   ,fDetector(-1)
48   ,fNtracklets(0)
49   ,fNplanes(0)
50   ,fNclusters(0)
51   ,fPt(0)
52   ,fPhi(0)
53   ,fEta(0)
54   ,fLabel(-1)
55   ,fPID(0)
56   ,fIsElectron(kFALSE)
57 {
58   //
59   // Default constructor
60   //
61
62 }
63
64 //_____________________________________________________________________________
65 AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack &t)
66   :TObject(t)
67   ,fTracklets(NULL)
68   ,fYproj(t.fYproj)
69   ,fZproj(t.fZproj)
70   ,fSlope(t.fSlope)
71   ,fDetector(t.fDetector)
72   ,fNtracklets(t.fNtracklets)
73   ,fNplanes(t.fNplanes)
74   ,fNclusters(t.fNclusters)
75   ,fPt(t.fPt)
76   ,fPhi(t.fPhi)
77   ,fEta(t.fEta)
78   ,fLabel(t.fLabel)
79   ,fPID(t.fPID)
80   ,fIsElectron(t.fIsElectron)
81 {
82   //
83   // Copy contructor
84   //
85
86 }
87
88 //_____________________________________________________________________________
89 AliTRDgtuTrack &AliTRDgtuTrack::operator=(const AliTRDgtuTrack &t)
90 {
91   //
92   // Assignment operator
93   //
94
95   if (this != &t) ((AliTRDgtuTrack &) t).Copy(*this); 
96
97   return *this;
98
99
100
101 //_____________________________________________________________________________
102 void AliTRDgtuTrack::Copy(TObject &t) const
103 {
104   //
105   // Copy function
106   //
107
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;
122
123 }
124
125 //_____________________________________________________________________________
126 AliTRDgtuTrack::~AliTRDgtuTrack()
127 {
128   //
129   // Destructor
130   //
131
132   if (fTracklets) {
133     //fTracklets->Delete();
134     fTracklets->Clear();
135     delete fTracklets;
136     fTracklets = 0;
137   }
138
139 }
140
141 //_____________________________________________________________________________
142 void AliTRDgtuTrack::AddTracklet(AliTRDltuTracklet *trk)
143 {
144   //
145   // Add a LTU tracklet to this track
146   //
147
148   Tracklets()->Add(trk);
149
150 }
151
152 //_____________________________________________________________________________
153 AliTRDltuTracklet *AliTRDgtuTrack::GetTracklet(Int_t pos) const
154 {
155   //
156   // Return LTU tracklet at position "pos"
157   //
158
159   if (fTracklets == 0) {
160     return 0;
161   }
162   void *trk = fTracklets->UncheckedAt(pos);
163   if (trk == 0) {
164     return 0;
165   }
166
167   return (AliTRDltuTracklet *) trk;
168
169 }
170
171 //_____________________________________________________________________________
172 Int_t AliTRDgtuTrack::Compare(const TObject *o) const
173 {
174   //
175   // Compare function for sorting the tracks
176   //
177
178   AliTRDgtuTrack *gtutrack = (AliTRDgtuTrack *) o;
179
180   if (fYproj <  gtutrack->GetYproj()) {
181     return -1;
182   }
183   if (fYproj == gtutrack->GetYproj()) {
184     return  0;
185   }
186
187   return +1;
188
189 }
190
191 //_____________________________________________________________________________
192 void AliTRDgtuTrack::Reset()
193 {
194   //
195   // Reset the track information
196   //
197
198   fYproj      = 0.0;
199   fZproj      = 0.0;
200   fSlope      = 0.0;
201   fDetector   = -1;
202   fNtracklets = 0;
203   fNplanes    = 0;
204   fNclusters  = 0;
205   fPt         = 0.0;
206   fPhi        = 0.0;
207   fEta        = 0.0;
208   fLabel      = -1;
209   fPID        = 0.0;
210   fIsElectron = kFALSE;
211   
212 }
213
214 //_____________________________________________________________________________
215 void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
216 {
217   //
218   // Calculate the kinematics of the found track
219   //
220
221   Int_t  i    = 0;
222   Int_t  iDet = 0;
223   Int_t  nDet = 0;
224   Bool_t newDetector;
225
226   AliTRDltuTracklet *trk;
227   Int_t nTracklets = GetNtracklets();
228   Float_t fC[kNmaxTrk][3];            // X, Y, Z  coordinates of segments
229
230   fYproj      = 0.0;
231   fZproj      = 0.0;
232   fSlope      = 0.0;
233   fNclusters  = 0;
234   fNplanes    = 0;
235   fNtracklets = GetNtracklets();
236   Int_t inDetector[kNplan];
237   for (i = 0; i < kNplan; i++) {
238     inDetector[i] = -1;
239   }
240
241   for (i = 0; i < nTracklets; i++) {
242
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();
249
250     newDetector = kTRUE;
251     for (Int_t id = 0; id < nDet; id++) {
252       if (iDet == inDetector[id]) {
253         newDetector = kFALSE;
254         break;
255       }
256     }
257     if (newDetector) {
258       inDetector[nDet++] = iDet;
259       fNplanes++;
260     }
261
262     fC[i][0] = trk->GetTime0();
263     fC[i][1] = trk->GetOffset();
264     fC[i][2] = trk->GetRowz();
265
266   }
267   fYproj /= (Float_t) nTracklets;
268   fZproj /= (Float_t) nTracklets;
269   fSlope /= (Float_t) nTracklets;
270
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++) {
276     count[i] = kFALSE;
277   }
278
279   Int_t iXmin = -1;
280   Int_t j     =  0;
281   x[0] = y[0] = z[0] = 0.0;
282   while (j < nTracklets) {
283     iXmin = -1;
284     for (i = 0; i < nTracklets; i++) {
285       if (count[i]) continue;
286       if (iXmin == -1) {
287         iXmin = i;
288         continue;
289       }
290       if (fC[i][0] < fC[iXmin][0]) {
291         iXmin = i;
292       }
293     }
294     x[j+1] = fC[iXmin][0];
295     y[j+1] = fC[iXmin][1];
296     z[j+1] = fC[iXmin][2];
297     j++;
298     count[iXmin] = kTRUE;
299   }
300
301   TMatrixD smatrix(2,2);
302   TMatrixD sums(2,1);
303   TMatrixD res(2,1);
304   Double_t xv, yv;
305   Float_t a, b;
306
307   smatrix.Zero();
308   sums.Zero();
309   for (i = 0; i < nTracklets; i++) {
310     xv = (Double_t) x[i+1];
311     yv = (Double_t) y[i+1];
312     smatrix(0,0) += 1.0;
313     smatrix(1,1) += xv*xv;
314     smatrix(0,1) += xv;
315     smatrix(1,0) += xv;
316     sums(0,0)    += yv;
317     sums(1,0)    += xv*yv;
318   }
319   res = smatrix.Invert() * sums;
320   a   = res(0,0);
321   b   = res(1,0);
322
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);
331
332   fPt = 0.3 * field * 0.01 * r;
333
334   Float_t d1 = fx1*fx1 + fy1*fy1;
335   Float_t d2 = fx2*fx2 + fy2*fy2;
336   Float_t d  = fx1*fy2 - fx2*fy1;
337   
338   Float_t xc = (d1*fy2 - d2*fy1) / (2.0*d);
339   Float_t yc = (d2*fx1 - d1*fx2) / (2.0*d);
340
341   if (yc != 0.0) {
342     fPhi = TMath::ATan(xc/yc);
343   } 
344   else {
345     fPhi = TMath::PiOver2();
346   }
347
348   fPhi *= 180.0/TMath::Pi();
349
350   smatrix.Zero();
351   sums.Zero();
352   for (i = 0; i < nTracklets+1; i++) {
353     xv = (Double_t) z[i];
354     yv = (Double_t) x[i];
355     smatrix(0,0) += 1.0;
356     smatrix(1,1) += xv*xv;
357     smatrix(0,1) += xv;
358     smatrix(1,0) += xv;
359     sums(0,0)    += yv;
360     sums(1,0)    += xv*yv;
361   }
362   res = smatrix.Invert() * sums;
363   a   = res(0,0);
364   b   = res(1,0);
365   Float_t theta = TMath::ATan(b);
366   
367   if (theta < 0.0) {
368     theta = TMath::Pi() + theta;
369   }
370   if (theta == 0.0) {
371     fEta = 0.0;
372   } 
373   else {
374     fEta = -TMath::Log(TMath::Tan(theta/2.0));
375   }
376   
377 }
378
379 //_____________________________________________________________________________
380 void AliTRDgtuTrack::MakePID()
381 {
382   //
383   // Electron likelihood signal
384   //
385
386   Int_t i = 0;
387
388   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
389   if (!calibration) {
390     AliError("No instance of AliTRDcalibDB.");
391     return;  
392   }
393   const AliTRDCalPID *pd = calibration->GetPIDObject(1);
394   
395   AliTRDltuTracklet *trk;
396   Int_t   nTracklets = GetNtracklets();
397   Int_t   det;
398   Int_t   pla;
399   Float_t sl;
400   Float_t th;
401   Float_t q;
402   Float_t probPio = 1.0;
403   Float_t probEle = 1.0;
404
405   for (i = 0; i < nTracklets; i++) {
406
407     trk = GetTracklet(i);
408
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));
412
413     q   = trk->GetQ() 
414         * TMath::Cos(sl/180.0*TMath::Pi()) 
415         * TMath::Cos(th/180.0*TMath::Pi());
416
417     det = trk->GetDetector();
418     pla = trk->GetPlane(det);
419
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
423     //q = q * 5.8;
424
425     // HEAD28Mar06
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
431     // as folows:
432     /*
433       if (k == 0) {
434       // electrons
435       q = 4.3 * q + 95.0;
436       } else {
437       // others
438       q = 4.2 * q + 70.0;
439       }
440     */
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];
447     // version v0_s1
448
449     // HEAD28Mar06 new distributions: factor = 6.0
450
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
454                 Float_t dedx[3];
455                 dedx[0] = dedx[1] = q*3.; dedx[2] = 0.;
456                 Float_t length = 3.7;
457
458     probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length, 0);
459     probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length, 0);
460
461   }
462
463   if ((probEle + probPio) > 0.0) {
464     fPID = probEle/(probEle + probPio);
465   } 
466   else {
467     fPID = 0.0;
468   }
469
470   // Thresholds for LQ cut at 90% electron efficiency (from AliRoot v4-03-07)
471   // P [GeV/c]  fPIDcut (between 0 and 1)
472   // 2          0.925
473   // 3          0.915
474   // 4          0.875
475   // 5          0.855
476   // 6          0.845
477   // 8          0.785
478   // 10         0.735
479   //
480   // PIDcut = 0.978 - 0.024 * P[GeV/c]
481   //Float_t PIDcut = 0.978 - 0.024 * TMath::Abs(fPt);
482
483   // HEAD28Mar06 with modified distributions (A. Bercuci changes, P. Shukla distributions)
484   //Float_t fPIDcut = 0.829 - 0.032 * TMath::Abs(fPt);
485
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
490
491   fIsElectron = kFALSE;
492   if (fPID >= fPIDcut) {
493     fIsElectron = kTRUE;
494   }
495
496 }
497
498 //_____________________________________________________________________________
499 void AliTRDgtuTrack::CookLabel()
500 {
501   //
502   // Cook the track label from tracklets labels
503   //
504
505   AliTRDltuTracklet *trk;
506
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++) {
511     trackLabel[it] = -1;
512     trackCount[it] =  0;
513   }
514
515   Bool_t counted;
516   Int_t  label;
517   Int_t  nTracks = 0;
518   for (Int_t itrk = 0; itrk < fNtracklets; itrk++) {
519
520     trk = GetTracklet(itrk);
521
522     if (trk->GetLabel() == -1) {
523       continue;
524     }
525
526     label = trk->GetLabel();
527
528     counted = kFALSE;
529     for (Int_t it = 0; it < nTracks; it++) {
530       if (label == trackLabel[it]) {
531         trackCount[it]++;
532         counted = kTRUE;
533         break;
534       }
535     }
536     if (!counted) {
537       trackLabel[nTracks] = label;
538       trackCount[nTracks]++;
539       nTracks++;
540       if (nTracks == kMaxTracks) {
541         Warning("CookLabel","Too many tracks for this tracklet.");
542         nTracks--;
543         break;
544       }
545     }
546     
547   }
548
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];
553       break;
554     }
555   }
556
557 }
558
559 //_____________________________________________________________________________
560 void AliTRDgtuTrack::ResetTracklets() 
561
562   //
563   // Resets the list of tracklets
564   //
565
566   if (fTracklets) {
567     //fTracklets->Delete(); 
568     fTracklets->Clear();
569   }
570
571 }
572
573 //_____________________________________________________________________________
574 TObjArray *AliTRDgtuTrack::Tracklets() 
575
576   //
577   // Returns the list of tracklets
578   //
579
580   if (!fTracklets) {
581     fTracklets = new TObjArray(400); 
582   }
583
584   return fTracklets; 
585
586 }
587
588 //_____________________________________________________________________________
589 Int_t AliTRDgtuTrack::GetNtracklets() const 
590 {
591   //
592   // Returns the number of tracklets
593   //
594
595   if (fTracklets) {
596     return fTracklets->GetEntriesFast();
597   }
598
599   return 0;
600
601 }