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