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