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