1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // The TRD particle identification base class //
22 // Its main purposes are: //
23 // - Provide I/O framework for all neccessary files //
24 // - Assignment of a e/pi propability to a given track //
26 ///////////////////////////////////////////////////////////////////////////////
33 #include <TObjArray.h>
36 #include <TParticle.h>
39 #include "AliTRDpid.h"
40 #include "AliTRDcluster.h"
41 #include "AliTRDtrack.h"
42 #include "AliTRDtracker.h"
43 #include "AliTRDgeometry.h"
48 //_____________________________________________________________________________
49 AliTRDpid::AliTRDpid():TNamed()
52 // AliTRDpid default constructor
61 fPIDpurePoints = kFALSE;
67 fThreePadOnly = kFALSE;
71 //_____________________________________________________________________________
72 AliTRDpid::AliTRDpid(const char* name, const char* title):TNamed(name,title)
75 // AliTRDpid constructor
89 //_____________________________________________________________________________
90 AliTRDpid::AliTRDpid(const AliTRDpid &p):TNamed(p)
93 // AliTRDpid copy constructor
96 ((AliTRDpid &) p).Copy(*this);
100 //_____________________________________________________________________________
101 AliTRDpid::~AliTRDpid()
104 // AliTRDpid destructor
108 fClusterArray->Delete();
109 delete fClusterArray;
110 fClusterArray = NULL;
114 fTrackArray->Delete();
126 //_____________________________________________________________________________
127 AliTRDpid &AliTRDpid::operator=(const AliTRDpid &p)
130 // Assignment operator
133 if (this != &p) ((AliTRDpid &) p).Copy(*this);
138 //_____________________________________________________________________________
139 void AliTRDpid::Copy(TObject &p) const
145 ((AliTRDpid &) p).fTrackArray = NULL;
146 ((AliTRDpid &) p).fClusterArray = NULL;
147 ((AliTRDpid &) p).fGeometry = NULL;
148 ((AliTRDpid &) p).fFileKine = NULL;
149 ((AliTRDpid &) p).fPIDratioMin = fPIDratioMin;
150 ((AliTRDpid &) p).fPIDpurePoints = fPIDpurePoints;
151 ((AliTRDpid &) p).fPIDindexMin = fPIDindexMin;
152 ((AliTRDpid &) p).fPIDindexMax = fPIDindexMax;
153 ((AliTRDpid &) p).fThreePadOnly = fThreePadOnly;
154 ((AliTRDpid &) p).fEvent = fEvent;
158 //_____________________________________________________________________________
159 Bool_t AliTRDpid::Init()
162 // Initializes the PID object
165 fClusterArray = new TObjArray();
166 fTrackArray = new TObjArray();
169 fPIDpurePoints = kTRUE;
173 fThreePadOnly = kFALSE;
179 //_____________________________________________________________________________
180 Bool_t AliTRDpid::AssignLikelihood()
183 // Assigns the e / pi likelihood to all tracks.
186 return AssignLikelihood(fTrackArray);
190 //_____________________________________________________________________________
191 Bool_t AliTRDpid::AssignLikelihood(TObjArray *tarray)
194 // Assigns the e / pi likelihood to all tracks in the array <tarray>.
197 Bool_t status = kTRUE;
201 TIter nextTrack(tarray);
202 while ((track = (AliTRDtrack *) nextTrack())) {
203 if (!AssignLikelihood(track)) {
213 //_____________________________________________________________________________
214 Bool_t AliTRDpid::FillSpectra()
217 // Fills the energy loss distributions with all tracks.
220 return FillSpectra(fTrackArray);
224 //_____________________________________________________________________________
225 Bool_t AliTRDpid::FillSpectra(TObjArray *tarray)
228 // Fills the energy loss distributions with all tracks in <tarray>.
231 Bool_t status = kTRUE;
235 TIter nextTrack(tarray);
236 while ((track = (AliTRDtrack *) nextTrack())) {
237 if (!FillSpectra(track)) {
247 //_____________________________________________________________________________
248 Bool_t AliTRDpid::Open(const Char_t *name, Int_t event)
251 // Opens and reads the kine tree, the geometry, the cluster array
252 // and the track array from the file <name>.
255 Bool_t status = kTRUE;
257 status = ReadCluster(name);
258 status = ReadTracks(name);
259 status = ReadKine(name,event);
265 //_____________________________________________________________________________
266 Bool_t AliTRDpid::Open(const Char_t *namekine
267 , const Char_t *namecluster
268 , const Char_t *nametracks
272 // Opens and reads the kine tree and the geometry from file <namekine>,
273 // the cluster array from file <namecluster>,
274 // and the track array from the file <nametracks>.
277 Bool_t status = kTRUE;
279 status = ReadCluster(namecluster);
280 status = ReadTracks(nametracks);
281 status = ReadKine(namekine,event);
287 //_____________________________________________________________________________
288 Bool_t AliTRDpid::ReadKine(const Char_t *name, Int_t event)
291 // Opens and reads the kine tree and the geometry from the file <name>.
294 TFile *fFileKine = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
296 printf("AliTRDpid::ReadKine -- ");
297 printf("Open file %s\n",name);
298 fFileKine = new TFile(name);
300 printf("AliTRDpid::ReadKine -- ");
301 printf("Cannot open file %s\n",name);
306 gAlice = (AliRun *) fFileKine->Get("gAlice");
308 printf("AliTRDpid::ReadKine -- ");
309 printf("No AliRun object found\n");
312 gAlice->GetEvent(event);
314 fGeometry = AliTRDgeometry::GetGeometry(gAlice->GetRunLoader());
316 printf("AliTRDpid::ReadKine -- ");
317 printf("No TRD geometry found\n");
327 //_____________________________________________________________________________
328 Bool_t AliTRDpid::ReadCluster(const Char_t *name)
331 // Opens and reads the cluster array from the file <name>.
334 TDirectory *savedir = gDirectory;
337 fClusterArray->Delete();
340 fClusterArray = new TObjArray();
343 printf("AliTRDpid::ReadCluster -- ");
344 printf("Open file %s\n",name);
346 AliTRDtracker *tracker = new AliTRDtracker();
347 TFile* file = TFile::Open(name);
349 TTree* tree = (TTree*) file->Get("TreeD");
350 tracker->ReadClusters(fClusterArray,tree);
354 if (!fClusterArray) {
355 printf("AliTRDpid::ReadCluster -- ");
356 printf("Error reading the cluster array from file %s\n",name);
368 //_____________________________________________________________________________
369 Bool_t AliTRDpid::ReadTracks(const Char_t *name)
372 // Opens and reads the track array from the file <name>.
375 TDirectory *savedir = gDirectory;
378 fTrackArray->Delete();
381 fTrackArray = new TObjArray();
384 TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
386 printf("AliTRDpid::ReadTracks -- ");
387 printf("Open file %s\n",name);
388 file = new TFile(name);
390 printf("AliTRDpid::ReadTracks -- ");
391 printf("Cannot open file %s\n",name);
397 sprintf(treeName,"TreeT%d_TRD",fEvent);
398 TTree *trackTree = (TTree *) file->Get(treeName);
399 TBranch *trackBranch = trackTree->GetBranch("tracks");
401 Int_t nEntry = ((Int_t) trackTree->GetEntries());
402 for (Int_t iEntry = 0; iEntry < nEntry; iEntry++) {
403 AliTRDtrack *track = new AliTRDtrack();
404 trackBranch->SetAddress(&track);
405 trackTree->GetEvent(iEntry);
406 fTrackArray->AddLast(track);
417 //_____________________________________________________________________________
418 Int_t AliTRDpid::MCpid(const AliTRDtrack *t)
421 // Determines the pid of the MC track <t>.
425 const Int_t kPdgEl = 11;
426 const Int_t kPdgPi = 211;
428 AliTRDcluster *cluster;
431 Int_t nClusterEl = 0;
432 Int_t nClusterPi = 0;
433 Int_t nClusterPure = 0;
434 Int_t nClusterMix = 0;
441 if (!fClusterArray) {
442 printf("AliTRDpid::MCpid -- ");
443 printf("ClusterArray not defined. Initialize the PID object first\n");
447 // Loop through all clusters associated to this track
448 Int_t nCluster = t->GetNumberOfClusters();
449 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
452 Int_t index = t->GetClusterIndex(iCluster);
453 if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
457 // Get the first two MC track indices
458 Int_t track0 = cluster->GetLabel(0);
459 Int_t track1 = cluster->GetLabel(1);
461 // Check on the track index to find the right primaries
462 if ((track0 > fPIDindexMin) &&
463 (track0 <= fPIDindexMax)) {
465 // Check whether it is a pure cluster, i.e. not overlapping
466 Bool_t accept = kTRUE;
467 if ((fPIDpurePoints) && (track1 != -1)) {
472 particle = gAlice->GetMCApp()->Particle(track0);
473 if (particle->GetFirstMother() == -1) {
474 switch (TMath::Abs(particle->GetPdgCode())) {
497 ratioEl = ((Float_t) nClusterEl) / ((Float_t) nCluster);
498 ratioPi = ((Float_t) nClusterPi) / ((Float_t) nCluster);
499 if (ratioEl > fPIDratioMin) {
502 else if (ratioPi > fPIDratioMin) {
506 // printf("AliTRDpid::MCpid -- ");
507 // printf("nCluster = %d, nClusterEl = %d, nClusterPi = %d\n"
508 // ,nCluster,nClusterEl,nClusterPi);
509 // printf("AliTRDpid::MCpid -- ");
510 // printf("nClusterPure = %d, nClusterMix = %d\n"
511 // ,nClusterPure,nClusterMix);
512 // printf("AliTRDpid::MCpid -- ");
513 // printf("ratioEl = %f, ratioPi = %f, ipid = %d\n",ratioEl,ratioPi,ipid);
519 //_____________________________________________________________________________
520 Int_t AliTRDpid::MCpid(const AliTRDtrack *t, Int_t *pdg
525 // Determines the pid of the MC track <t>.
526 // Returns the number of different MC particles that contribute to this
527 // track. <pdg> contains their PDG code, ordered by their frequency.
528 // <nFound> contains the corresponding number of cluster that contain
529 // contributions to this cluster.
532 const Int_t kNtrack = 3;
533 const Int_t kNpart = 200;
535 AliTRDcluster *cluster;
543 if (!fClusterArray) {
544 printf("AliTRDpid::MCpid -- ");
545 printf("ClusterArray not defined. Initialize the PID object first\n");
550 Int_t pdgTmp[kNpart];
551 Int_t nFoundTmp[kNpart];
552 Int_t indicesTmp[kNpart];
553 for (iPart = 0; iPart < kNpart; iPart++) {
558 nFoundTmp[iPart] = 0;
559 indicesTmp[iPart] = 0;
563 // Loop through all clusters associated to this track
564 Int_t nCluster = t->GetNumberOfClusters();
565 for (iCluster = 0; iCluster < nCluster; iCluster++) {
568 Int_t index = t->GetClusterIndex(iCluster);
569 if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
573 // Get the MC track indices
574 for (iTrack = 0; iTrack < kNtrack; iTrack++) {
576 Int_t trackIndex = cluster->GetLabel(iTrack);
577 if (trackIndex >= 0) {
578 particle = gAlice->GetMCApp()->Particle(trackIndex);
579 Int_t pdgCode = particle->GetPdgCode();
580 Bool_t newPart = kTRUE;
581 for (iPart = 0; iPart < nPart; iPart++) {
582 if (trackIndex == indicesTmp[iPart]) {
588 if ((newPart) && (nPart < kNpart)) {
589 indicesTmp[nPart] = trackIndex;
590 pdgTmp[nPart] = pdgCode;
600 // Sort the arrays by the frequency of the appearance of a MC track
601 TMath::Sort(nPart,nFoundTmp,sort,kTRUE);
602 for (iPart = 0; iPart < nPart; iPart++) {
603 pdg[iPart] = pdgTmp[sort[iPart]];
604 nFound[iPart] = nFoundTmp[sort[iPart]];
605 indices[iPart] = indicesTmp[sort[iPart]];
606 // printf(" iPart = %d, pdg = %d, nFound = %d, indices = %d\n"
607 // ,iPart,pdg[iPart],nFound[iPart],indices[iPart]);
614 //_____________________________________________________________________________
615 Bool_t AliTRDpid::SumCharge(const AliTRDtrack *t
620 // Sums up the charge in each plane for track <t>.
623 Bool_t status = kTRUE;
625 AliTRDcluster *cluster;
627 const Int_t kNplane = AliTRDgeometry::Nplan();
628 for (Int_t iPlane = 0; iPlane < kNplane; iPlane++) {
629 charge[iPlane] = 0.0;
630 nCluster[iPlane] = 0;
633 if (!fClusterArray) {
634 printf("AliTRDpid::SumCharge -- ");
635 printf("ClusterArray not defined. Initialize the PID object first\n");
639 printf("AliTRDpid::SumCharge -- ");
640 printf("TRD geometry not defined. Initialize the PID object first\n");
644 // Loop through all clusters associated to this track
645 Int_t nClus = t->GetNumberOfClusters();
646 for (Int_t iClus = 0; iClus < nClus; iClus++) {
649 Int_t index = t->GetClusterIndex(iClus);
650 if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
656 if ((!cluster->From2pad()) &&
657 (!cluster->From3pad())) continue;
661 Int_t plane = fGeometry->GetPlane(cluster->GetDetector());
662 //charge[plane] += cluster->GetQ();
664 charge[plane] += t->GetClusterdQdl(iClus);
668 // printf("AliTRDpid::SumCharge -- ");
669 // printf("charge = %d, %d, %d, %d, %d, %d\n"
670 // ,(Int_t) charge[0]
671 // ,(Int_t) charge[1]
672 // ,(Int_t) charge[2]
673 // ,(Int_t) charge[3]
674 // ,(Int_t) charge[4]
675 // ,(Int_t) charge[5]);
676 // printf("AliTRDpid::SumCharge -- ");
677 // printf("nCluster = %d, %d, %d, %d, %d, %d\n"