6f06f7138ac4f31cef69247b7e04141691571a63
[u/mrichter/AliRoot.git] / TRD / AliTRDpid.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 /*
17 $Log$
18 Revision 1.3  2001/11/14 10:50:46  cblume
19 Changes in digits IO. Add merging of summable digits
20
21 Revision 1.2  2001/11/06 17:19:41  cblume
22 Add detailed geometry and simple simulator
23
24 Revision 1.1  2001/05/07 08:08:05  cblume
25 Update of TRD code
26
27 */
28
29 ///////////////////////////////////////////////////////////////////////////////
30 //                                                                           //
31 //   The TRD particle identification base class                              //
32 //                                                                           //
33 //   Its main purposes are:                                                  //
34 //      - Provide I/O framework for all neccessary files                     //
35 //      - Assignment of a e/pi propability to a given track                  //
36 //                                                                           //
37 ///////////////////////////////////////////////////////////////////////////////
38
39 #include <stdlib.h>
40 #include <math.h>
41
42 #include <TROOT.h>
43 #include <TH1.h>
44 #include <TObjArray.h>
45 #include <TTree.h>
46 #include <TFile.h>
47 #include <TParticle.h>
48
49 #include "AliRun.h"
50 #include "AliTRD.h"
51 #include "AliTRDpid.h"
52 #include "AliTRDcluster.h"
53 #include "AliTRDtrack.h"
54 #include "AliTRDtracker.h"
55 #include "AliTRDgeometry.h"
56
57 ClassImp(AliTRDpid)
58
59 //_____________________________________________________________________________
60 AliTRDpid::AliTRDpid():TNamed()
61 {
62   //
63   // AliTRDpid default constructor
64   // 
65
66   fTrackArray    = NULL;
67   fClusterArray  = NULL;
68   fGeometry      = NULL;
69   fFileKine      = NULL;
70
71   fPIDratioMin   = 0.0;
72   fPIDpurePoints = kFALSE;
73   fPIDindexMin   = 0;
74   fPIDindexMax   = 0;
75
76   fEvent         = 0;
77
78   fThreePadOnly  = kFALSE;
79
80 }
81
82 //_____________________________________________________________________________
83 AliTRDpid::AliTRDpid(const char* name, const char* title):TNamed(name,title)
84 {
85   //
86   // AliTRDpid constructor
87   //
88
89   fTrackArray   = NULL;
90   fClusterArray = NULL;
91   fGeometry     = NULL;
92   fFileKine     = NULL;
93
94   fEvent        = 0;
95
96   Init();
97
98 }
99
100 //_____________________________________________________________________________
101 AliTRDpid::AliTRDpid(const AliTRDpid &p)
102 {
103   //
104   // AliTRDpid copy constructor
105   //
106
107   ((AliTRDpid &) p).Copy(*this);
108
109 }
110
111 //_____________________________________________________________________________
112 AliTRDpid::~AliTRDpid()
113 {
114   //
115   // AliTRDpid destructor
116   //
117
118   if (fClusterArray) {
119     fClusterArray->Delete();
120     delete fClusterArray;
121     fClusterArray = NULL;
122   }
123
124   if (fTrackArray) {
125     fTrackArray->Delete();
126     delete fTrackArray;
127     fTrackArray = NULL;
128   }
129
130   if (fFileKine) {
131     delete fFileKine;
132     fFileKine = NULL;
133   }
134
135 }
136
137 //_____________________________________________________________________________
138 AliTRDpid &AliTRDpid::operator=(const AliTRDpid &p)
139 {
140   //
141   // Assignment operator
142   //
143
144   if (this != &p) ((AliTRDpid &) p).Copy(*this);
145   return *this;
146
147 }
148
149 //_____________________________________________________________________________
150 void AliTRDpid::Copy(TObject &p)
151 {
152   //
153   // Copy function
154   //
155
156   ((AliTRDpid &) p).fTrackArray    = NULL;    
157   ((AliTRDpid &) p).fClusterArray  = NULL;    
158   ((AliTRDpid &) p).fGeometry      = NULL;    
159   ((AliTRDpid &) p).fFileKine      = NULL;
160   ((AliTRDpid &) p).fPIDratioMin   = fPIDratioMin;
161   ((AliTRDpid &) p).fPIDpurePoints = fPIDpurePoints;
162   ((AliTRDpid &) p).fPIDindexMin   = fPIDindexMin;
163   ((AliTRDpid &) p).fPIDindexMax   = fPIDindexMax;
164   ((AliTRDpid &) p).fThreePadOnly  = fThreePadOnly;
165   ((AliTRDpid &) p).fEvent         = fEvent;
166
167 }
168
169 //_____________________________________________________________________________
170 Bool_t AliTRDpid::Init()
171 {
172   //
173   // Initializes the PID object 
174   //
175
176   fClusterArray  = new TObjArray();
177   fTrackArray    = new TObjArray();
178
179   fPIDratioMin   = 0.7;
180   fPIDpurePoints = kTRUE;
181   fPIDindexMin   = -1;
182   fPIDindexMax   = -1;
183
184   fThreePadOnly  = kFALSE;
185
186   return kTRUE;
187
188 }
189
190 //_____________________________________________________________________________
191 Bool_t AliTRDpid::AssignLikelihood()
192 {
193   //
194   // Assigns the e / pi likelihood to all tracks.
195   //
196
197   return AssignLikelihood(fTrackArray);
198
199 }
200
201 //_____________________________________________________________________________
202 Bool_t AliTRDpid::AssignLikelihood(TObjArray *tarray)
203 {
204   //
205   // Assigns the e / pi likelihood to all tracks in the array <tarray>.
206   //
207
208   Bool_t status = kTRUE;
209
210   AliTRDtrack *track;
211
212   TIter nextTrack(tarray);  
213   while ((track = (AliTRDtrack *) nextTrack())) {
214     if (!AssignLikelihood(track)) {
215       status = kFALSE;
216       continue;
217     }
218   }
219
220   return status;
221
222 }
223
224 //_____________________________________________________________________________
225 Bool_t AliTRDpid::FillSpectra()
226 {
227   //
228   // Fills the energy loss distributions with all tracks.
229   //
230
231   return FillSpectra(fTrackArray);
232
233 }
234   
235 //_____________________________________________________________________________
236 Bool_t AliTRDpid::FillSpectra(TObjArray *tarray)
237 {
238   //
239   // Fills the energy loss distributions with all tracks in <tarray>.
240   //
241   
242   Bool_t status = kTRUE;
243
244   AliTRDtrack *track;
245
246   TIter nextTrack(tarray);
247   while ((track = (AliTRDtrack *) nextTrack())) {
248     if (!FillSpectra(track)) {
249       status = kFALSE;
250       continue;
251     }
252   }  
253
254   return kTRUE;
255
256 }
257
258 //_____________________________________________________________________________
259 Bool_t AliTRDpid::Open(const Char_t *name, Int_t event)
260 {
261   //
262   // Opens and reads the kine tree, the geometry, the cluster array
263   // and the track array from the file <name>.
264   //
265
266   Bool_t status = kTRUE;
267
268   status = ReadCluster(name);
269   status = ReadTracks(name);
270   status = ReadKine(name,event);
271
272   return status;
273
274 }
275
276 //_____________________________________________________________________________
277 Bool_t AliTRDpid::Open(const Char_t *namekine
278                      , const Char_t *namecluster
279                      , const Char_t *nametracks
280                      , Int_t event)
281 {
282   //
283   // Opens and reads the kine tree and the geometry from file <namekine>,
284   // the cluster array from file <namecluster>,
285   // and the track array from the file <nametracks>.
286   //
287
288   Bool_t status = kTRUE;
289
290   status = ReadCluster(namecluster);
291   status = ReadTracks(nametracks);
292   status = ReadKine(namekine,event);
293
294   return status;
295
296 }
297
298 //_____________________________________________________________________________
299 Bool_t AliTRDpid::ReadKine(const Char_t *name, Int_t event)
300 {
301   //
302   // Opens and reads the kine tree and the geometry from the file <name>.
303   //
304
305   TFile *fFileKine = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
306   if (!fFileKine) {
307     printf("AliTRDpid::ReadKine -- ");
308     printf("Open file %s\n",name);
309     fFileKine = new TFile(name);
310     if (!fFileKine) {
311       printf("AliTRDpid::ReadKine -- ");
312       printf("Cannot open file %s\n",name);
313       return kFALSE;
314     }
315   }
316
317   gAlice = (AliRun *) fFileKine->Get("gAlice");
318   if (!gAlice) {
319     printf("AliTRDpid::ReadKine -- ");
320     printf("No AliRun object found\n");    
321     return kFALSE;
322   }
323   gAlice->GetEvent(event);
324
325   AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
326   if (!trd) {
327     printf("AliTRDpid::ReadKine -- ");
328     printf("No TRD object found\n");    
329     return kFALSE;
330   }
331
332   fGeometry = trd->GetGeometry();
333   if (!fGeometry) {
334     printf("AliTRDpid::ReadKine -- ");
335     printf("No TRD geometry found\n");
336     return kFALSE;
337   }
338
339   fEvent = event;
340
341   return kTRUE;
342
343 }
344
345 //_____________________________________________________________________________
346 Bool_t AliTRDpid::ReadCluster(const Char_t *name)
347 {
348   //
349   // Opens and reads the cluster array from the file <name>.
350   //
351  
352   TDirectory *savedir = gDirectory;                                                   
353
354   if (fClusterArray) {
355     fClusterArray->Delete();
356   }
357   else {
358     fClusterArray = new TObjArray();
359   }
360
361   printf("AliTRDpid::ReadCluster -- ");
362   printf("Open file %s\n",name);
363
364   AliTRDtracker *tracker = new AliTRDtracker("dummy","dummy");
365   tracker->ReadClusters(fClusterArray,name);
366
367   if (!fClusterArray) {
368     printf("AliTRDpid::ReadCluster -- ");
369     printf("Error reading the cluster array from file %s\n",name);
370     return kFALSE;
371   }
372
373   delete tracker;
374
375   savedir->cd();
376
377   return kTRUE;
378
379 }
380
381 //_____________________________________________________________________________
382 Bool_t AliTRDpid::ReadTracks(const Char_t *name)
383 {
384   //
385   // Opens and reads the track array from the file <name>.
386   //
387  
388   TDirectory *savedir = gDirectory;                                                   
389
390   if (fTrackArray) {
391     fTrackArray->Delete();
392   }
393   else {
394     fTrackArray = new TObjArray();
395   }
396
397   TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
398   if (!file) {
399     printf("AliTRDpid::ReadTracks -- ");
400     printf("Open file %s\n",name);
401     file = new TFile(name);
402     if (!file) {
403       printf("AliTRDpid::ReadTracks -- ");
404       printf("Cannot open file %s\n",name);
405       return kFALSE;
406     }
407   }
408
409   Char_t treeName[12];
410   sprintf(treeName,"TreeT%d_TRD",fEvent);
411   TTree   *trackTree   = (TTree *) file->Get(treeName);
412   TBranch *trackBranch = trackTree->GetBranch("tracks");
413
414   Int_t nEntry = ((Int_t) trackTree->GetEntries());
415   for (Int_t iEntry = 0; iEntry < nEntry; iEntry++) {
416     AliTRDtrack *track = new AliTRDtrack();
417     trackBranch->SetAddress(&track);
418     trackTree->GetEvent(iEntry);
419     fTrackArray->AddLast(track);
420   }
421
422   file->Close();
423
424   savedir->cd();
425
426   return kTRUE;
427
428 }
429
430 //_____________________________________________________________________________
431 Int_t AliTRDpid::MCpid(const AliTRDtrack *t)
432 {
433   //
434   // Determines the pid of the MC track <t>.
435   //
436
437   // PDG codes
438   const Int_t kPdgEl =  11;
439   const Int_t kPdgPi = 211;
440
441   AliTRDcluster *cluster;
442   TParticle     *particle;
443
444   Int_t   nClusterEl   = 0;
445   Int_t   nClusterPi   = 0;
446   Int_t   nClusterPure = 0;
447   Int_t   nClusterMix  = 0;
448
449   Float_t ratioEl      = 0;
450   Float_t ratioPi      = 0;
451
452   Int_t   ipid         = -1;
453
454   if (!fClusterArray) {
455     printf("AliTRDpid::MCpid -- ");
456     printf("ClusterArray not defined. Initialize the PID object first\n");
457     return -1;  
458   }
459   
460   // Loop through all clusters associated to this track
461   Int_t nCluster = t->GetNclusters();
462   for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
463
464     // Get a cluster
465     Int_t index = t->GetClusterIndex(iCluster);
466     if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
467       break;
468     } 
469
470     // Get the first two MC track indices
471     Int_t track0 = cluster->GetTrackIndex(0);
472     Int_t track1 = cluster->GetTrackIndex(1);
473
474     // Check on the track index to find the right primaries
475     if ((track0 >  fPIDindexMin) && 
476         (track0 <= fPIDindexMax)) {
477
478       // Check whether it is a pure cluster, i.e. not overlapping
479       Bool_t accept = kTRUE;
480       if ((fPIDpurePoints) && (track1 != -1)) {
481         accept = kFALSE;
482       }
483       if (accept) {
484
485         particle = gAlice->Particle(track0);
486         if (particle->GetFirstMother() == -1) {
487           switch (TMath::Abs(particle->GetPdgCode())) {
488           case kPdgEl:
489             nClusterEl++;
490             break;
491           case kPdgPi:
492             nClusterPi++;
493             break;
494           };
495           nClusterPure++;
496         }
497
498       }
499       else {
500  
501         nClusterMix++;
502
503       }
504
505     }
506
507   }
508   
509   if (nCluster) {
510     ratioEl = ((Float_t) nClusterEl) / ((Float_t) nCluster);
511     ratioPi = ((Float_t) nClusterPi) / ((Float_t) nCluster);
512     if      (ratioEl > fPIDratioMin) {
513       ipid = kElectron;
514     }
515     else if (ratioPi > fPIDratioMin) {
516       ipid = kPion;
517     }
518   }
519 //   printf("AliTRDpid::MCpid -- ");
520 //   printf("nCluster = %d, nClusterEl = %d, nClusterPi = %d\n"
521 //       ,nCluster,nClusterEl,nClusterPi);
522 //   printf("AliTRDpid::MCpid -- ");
523 //   printf("nClusterPure = %d, nClusterMix = %d\n"
524 //       ,nClusterPure,nClusterMix);
525 //   printf("AliTRDpid::MCpid -- ");
526 //   printf("ratioEl = %f, ratioPi = %f, ipid = %d\n",ratioEl,ratioPi,ipid);
527
528   return ipid;
529
530 }
531
532 //_____________________________________________________________________________
533 Int_t AliTRDpid::MCpid(const AliTRDtrack *t, Int_t *pdg
534                                            , Int_t *nFound
535                                            , Int_t *indices)
536 {
537   //
538   // Determines the pid of the MC track <t>.
539   // Returns the number of different MC particles that contribute to this
540   // track. <pdg> contains their PDG code, ordered by their frequency.
541   // <nFound> contains the corresponding number of cluster that contain
542   // contributions to this cluster.
543   //
544
545   const Int_t kNtrack = 3;
546   const Int_t kNpart  = 200;
547
548   AliTRDcluster *cluster;
549   TParticle     *particle;
550
551   Int_t nPart    = 0;
552   Int_t iPart    = 0;
553   Int_t iTrack   = 0;
554   Int_t iCluster = 0;
555
556   if (!fClusterArray) {
557     printf("AliTRDpid::MCpid -- ");
558     printf("ClusterArray not defined. Initialize the PID object first\n");
559     return -1;  
560   }
561   
562   Int_t sort[kNpart];
563   Int_t pdgTmp[kNpart];
564   Int_t nFoundTmp[kNpart];
565   Int_t indicesTmp[kNpart];
566   for (iPart = 0; iPart < kNpart; iPart++) {
567     pdg[iPart]        = 0;
568     nFound[iPart]     = 0;
569     indices[iPart]    = 0;
570     pdgTmp[iPart]     = 0;
571     nFoundTmp[iPart]  = 0;
572     indicesTmp[iPart] = 0;
573     sort[iPart]       = 0;
574   }
575
576   // Loop through all clusters associated to this track
577   Int_t nCluster = t->GetNclusters();
578   for (iCluster = 0; iCluster < nCluster; iCluster++) {
579
580     // Get a cluster
581     Int_t index = t->GetClusterIndex(iCluster);
582     if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
583       break;
584     } 
585
586     // Get the MC track indices
587     for (iTrack = 0; iTrack < kNtrack; iTrack++) {
588
589       Int_t trackIndex = cluster->GetTrackIndex(iTrack);
590       if (trackIndex >= 0) {
591         particle = gAlice->Particle(trackIndex);
592         Int_t  pdgCode = particle->GetPdgCode(); 
593         Bool_t newPart = kTRUE;
594         for (iPart = 0; iPart < nPart; iPart++) {
595          if (trackIndex == indicesTmp[iPart]) {
596             nFoundTmp[iPart]++;
597             newPart = kFALSE;
598             break;
599           }
600         }
601         if ((newPart) && (nPart < kNpart)) {
602           indicesTmp[nPart] = trackIndex;
603           pdgTmp[nPart]     = pdgCode;
604           nFoundTmp[nPart]++;
605           nPart++;
606         }
607       }
608
609     }
610
611   }
612
613   // Sort the arrays by the frequency of the appearance of a MC track
614   TMath::Sort(nPart,nFoundTmp,sort,kTRUE);
615   for (iPart = 0; iPart < nPart; iPart++) {
616     pdg[iPart]     = pdgTmp[sort[iPart]];
617     nFound[iPart]  = nFoundTmp[sort[iPart]];
618     indices[iPart] = indicesTmp[sort[iPart]];
619 //     printf(" iPart = %d, pdg = %d, nFound = %d, indices = %d\n"
620 //        ,iPart,pdg[iPart],nFound[iPart],indices[iPart]);
621   }
622
623   return nPart;
624
625 }
626
627 //_____________________________________________________________________________
628 Bool_t AliTRDpid::SumCharge(const AliTRDtrack *t
629                            , Float_t *charge
630                            , Int_t   *nCluster)
631 {
632   //
633   // Sums up the charge in each plane for track <t>.
634   //
635
636   Bool_t status = kTRUE;
637
638   AliTRDcluster *cluster;
639
640   const Int_t kNplane = AliTRDgeometry::Nplan();
641   for (Int_t iPlane = 0; iPlane < kNplane; iPlane++) {
642     charge[iPlane]   = 0.0;
643     nCluster[iPlane] = 0;
644   }
645
646   if (!fClusterArray) {
647     printf("AliTRDpid::SumCharge -- ");
648     printf("ClusterArray not defined. Initialize the PID object first\n");
649     return kFALSE;  
650   }
651   if (!fGeometry) {
652     printf("AliTRDpid::SumCharge -- ");
653     printf("TRD geometry not defined. Initialize the PID object first\n");
654     return kFALSE;  
655   }
656   
657   // Loop through all clusters associated to this track
658   Int_t nClus = t->GetNclusters();
659   for (Int_t iClus = 0; iClus < nClus; iClus++) {
660
661     // Get a cluster
662     Int_t index = t->GetClusterIndex(iClus);
663     if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
664       status = kFALSE;
665       break;
666     } 
667
668     if (fThreePadOnly) {
669       if ((!cluster->From2pad()) &&
670           (!cluster->From3pad())) continue;
671     }
672
673     // Sum up the charge
674     Int_t plane    = fGeometry->GetPlane(cluster->GetDetector());
675     //charge[plane] += cluster->GetQ();
676     // Corrected charge
677     charge[plane] += t->GetClusterdQdl(iClus);
678     nCluster[plane]++;
679
680   }
681 //   printf("AliTRDpid::SumCharge -- ");
682 //   printf("charge = %d, %d, %d, %d, %d, %d\n"
683 //       ,(Int_t) charge[0]
684 //          ,(Int_t) charge[1]
685 //          ,(Int_t) charge[2]
686 //          ,(Int_t) charge[3]
687 //          ,(Int_t) charge[4]
688 //          ,(Int_t) charge[5]);
689 //   printf("AliTRDpid::SumCharge -- ");
690 //   printf("nCluster = %d, %d, %d, %d, %d, %d\n"
691 //          ,nCluster[0]
692 //          ,nCluster[1]
693 //          ,nCluster[2]
694 //          ,nCluster[3]
695 //          ,nCluster[4]
696 //          ,nCluster[5]);
697
698   return status;
699
700 }