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 Revision 1.3 2001/11/14 10:50:46 cblume
19 Changes in digits IO. Add merging of summable digits
21 Revision 1.2 2001/11/06 17:19:41 cblume
22 Add detailed geometry and simple simulator
24 Revision 1.1 2001/05/07 08:08:05 cblume
29 ///////////////////////////////////////////////////////////////////////////////
31 // The TRD particle identification base class //
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 //
37 ///////////////////////////////////////////////////////////////////////////////
44 #include <TObjArray.h>
47 #include <TParticle.h>
51 #include "AliTRDpid.h"
52 #include "AliTRDcluster.h"
53 #include "AliTRDtrack.h"
54 #include "AliTRDtracker.h"
55 #include "AliTRDgeometry.h"
59 //_____________________________________________________________________________
60 AliTRDpid::AliTRDpid():TNamed()
63 // AliTRDpid default constructor
72 fPIDpurePoints = kFALSE;
78 fThreePadOnly = kFALSE;
82 //_____________________________________________________________________________
83 AliTRDpid::AliTRDpid(const char* name, const char* title):TNamed(name,title)
86 // AliTRDpid constructor
100 //_____________________________________________________________________________
101 AliTRDpid::AliTRDpid(const AliTRDpid &p)
104 // AliTRDpid copy constructor
107 ((AliTRDpid &) p).Copy(*this);
111 //_____________________________________________________________________________
112 AliTRDpid::~AliTRDpid()
115 // AliTRDpid destructor
119 fClusterArray->Delete();
120 delete fClusterArray;
121 fClusterArray = NULL;
125 fTrackArray->Delete();
137 //_____________________________________________________________________________
138 AliTRDpid &AliTRDpid::operator=(const AliTRDpid &p)
141 // Assignment operator
144 if (this != &p) ((AliTRDpid &) p).Copy(*this);
149 //_____________________________________________________________________________
150 void AliTRDpid::Copy(TObject &p)
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;
169 //_____________________________________________________________________________
170 Bool_t AliTRDpid::Init()
173 // Initializes the PID object
176 fClusterArray = new TObjArray();
177 fTrackArray = new TObjArray();
180 fPIDpurePoints = kTRUE;
184 fThreePadOnly = kFALSE;
190 //_____________________________________________________________________________
191 Bool_t AliTRDpid::AssignLikelihood()
194 // Assigns the e / pi likelihood to all tracks.
197 return AssignLikelihood(fTrackArray);
201 //_____________________________________________________________________________
202 Bool_t AliTRDpid::AssignLikelihood(TObjArray *tarray)
205 // Assigns the e / pi likelihood to all tracks in the array <tarray>.
208 Bool_t status = kTRUE;
212 TIter nextTrack(tarray);
213 while ((track = (AliTRDtrack *) nextTrack())) {
214 if (!AssignLikelihood(track)) {
224 //_____________________________________________________________________________
225 Bool_t AliTRDpid::FillSpectra()
228 // Fills the energy loss distributions with all tracks.
231 return FillSpectra(fTrackArray);
235 //_____________________________________________________________________________
236 Bool_t AliTRDpid::FillSpectra(TObjArray *tarray)
239 // Fills the energy loss distributions with all tracks in <tarray>.
242 Bool_t status = kTRUE;
246 TIter nextTrack(tarray);
247 while ((track = (AliTRDtrack *) nextTrack())) {
248 if (!FillSpectra(track)) {
258 //_____________________________________________________________________________
259 Bool_t AliTRDpid::Open(const Char_t *name, Int_t event)
262 // Opens and reads the kine tree, the geometry, the cluster array
263 // and the track array from the file <name>.
266 Bool_t status = kTRUE;
268 status = ReadCluster(name);
269 status = ReadTracks(name);
270 status = ReadKine(name,event);
276 //_____________________________________________________________________________
277 Bool_t AliTRDpid::Open(const Char_t *namekine
278 , const Char_t *namecluster
279 , const Char_t *nametracks
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>.
288 Bool_t status = kTRUE;
290 status = ReadCluster(namecluster);
291 status = ReadTracks(nametracks);
292 status = ReadKine(namekine,event);
298 //_____________________________________________________________________________
299 Bool_t AliTRDpid::ReadKine(const Char_t *name, Int_t event)
302 // Opens and reads the kine tree and the geometry from the file <name>.
305 TFile *fFileKine = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
307 printf("AliTRDpid::ReadKine -- ");
308 printf("Open file %s\n",name);
309 fFileKine = new TFile(name);
311 printf("AliTRDpid::ReadKine -- ");
312 printf("Cannot open file %s\n",name);
317 gAlice = (AliRun *) fFileKine->Get("gAlice");
319 printf("AliTRDpid::ReadKine -- ");
320 printf("No AliRun object found\n");
323 gAlice->GetEvent(event);
325 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
327 printf("AliTRDpid::ReadKine -- ");
328 printf("No TRD object found\n");
332 fGeometry = trd->GetGeometry();
334 printf("AliTRDpid::ReadKine -- ");
335 printf("No TRD geometry found\n");
345 //_____________________________________________________________________________
346 Bool_t AliTRDpid::ReadCluster(const Char_t *name)
349 // Opens and reads the cluster array from the file <name>.
352 TDirectory *savedir = gDirectory;
355 fClusterArray->Delete();
358 fClusterArray = new TObjArray();
361 printf("AliTRDpid::ReadCluster -- ");
362 printf("Open file %s\n",name);
364 AliTRDtracker *tracker = new AliTRDtracker("dummy","dummy");
365 tracker->ReadClusters(fClusterArray,name);
367 if (!fClusterArray) {
368 printf("AliTRDpid::ReadCluster -- ");
369 printf("Error reading the cluster array from file %s\n",name);
381 //_____________________________________________________________________________
382 Bool_t AliTRDpid::ReadTracks(const Char_t *name)
385 // Opens and reads the track array from the file <name>.
388 TDirectory *savedir = gDirectory;
391 fTrackArray->Delete();
394 fTrackArray = new TObjArray();
397 TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
399 printf("AliTRDpid::ReadTracks -- ");
400 printf("Open file %s\n",name);
401 file = new TFile(name);
403 printf("AliTRDpid::ReadTracks -- ");
404 printf("Cannot open file %s\n",name);
410 sprintf(treeName,"TreeT%d_TRD",fEvent);
411 TTree *trackTree = (TTree *) file->Get(treeName);
412 TBranch *trackBranch = trackTree->GetBranch("tracks");
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);
430 //_____________________________________________________________________________
431 Int_t AliTRDpid::MCpid(const AliTRDtrack *t)
434 // Determines the pid of the MC track <t>.
438 const Int_t kPdgEl = 11;
439 const Int_t kPdgPi = 211;
441 AliTRDcluster *cluster;
444 Int_t nClusterEl = 0;
445 Int_t nClusterPi = 0;
446 Int_t nClusterPure = 0;
447 Int_t nClusterMix = 0;
454 if (!fClusterArray) {
455 printf("AliTRDpid::MCpid -- ");
456 printf("ClusterArray not defined. Initialize the PID object first\n");
460 // Loop through all clusters associated to this track
461 Int_t nCluster = t->GetNclusters();
462 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
465 Int_t index = t->GetClusterIndex(iCluster);
466 if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
470 // Get the first two MC track indices
471 Int_t track0 = cluster->GetTrackIndex(0);
472 Int_t track1 = cluster->GetTrackIndex(1);
474 // Check on the track index to find the right primaries
475 if ((track0 > fPIDindexMin) &&
476 (track0 <= fPIDindexMax)) {
478 // Check whether it is a pure cluster, i.e. not overlapping
479 Bool_t accept = kTRUE;
480 if ((fPIDpurePoints) && (track1 != -1)) {
485 particle = gAlice->Particle(track0);
486 if (particle->GetFirstMother() == -1) {
487 switch (TMath::Abs(particle->GetPdgCode())) {
510 ratioEl = ((Float_t) nClusterEl) / ((Float_t) nCluster);
511 ratioPi = ((Float_t) nClusterPi) / ((Float_t) nCluster);
512 if (ratioEl > fPIDratioMin) {
515 else if (ratioPi > fPIDratioMin) {
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);
532 //_____________________________________________________________________________
533 Int_t AliTRDpid::MCpid(const AliTRDtrack *t, Int_t *pdg
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.
545 const Int_t kNtrack = 3;
546 const Int_t kNpart = 200;
548 AliTRDcluster *cluster;
556 if (!fClusterArray) {
557 printf("AliTRDpid::MCpid -- ");
558 printf("ClusterArray not defined. Initialize the PID object first\n");
563 Int_t pdgTmp[kNpart];
564 Int_t nFoundTmp[kNpart];
565 Int_t indicesTmp[kNpart];
566 for (iPart = 0; iPart < kNpart; iPart++) {
571 nFoundTmp[iPart] = 0;
572 indicesTmp[iPart] = 0;
576 // Loop through all clusters associated to this track
577 Int_t nCluster = t->GetNclusters();
578 for (iCluster = 0; iCluster < nCluster; iCluster++) {
581 Int_t index = t->GetClusterIndex(iCluster);
582 if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
586 // Get the MC track indices
587 for (iTrack = 0; iTrack < kNtrack; iTrack++) {
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]) {
601 if ((newPart) && (nPart < kNpart)) {
602 indicesTmp[nPart] = trackIndex;
603 pdgTmp[nPart] = pdgCode;
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]);
627 //_____________________________________________________________________________
628 Bool_t AliTRDpid::SumCharge(const AliTRDtrack *t
633 // Sums up the charge in each plane for track <t>.
636 Bool_t status = kTRUE;
638 AliTRDcluster *cluster;
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;
646 if (!fClusterArray) {
647 printf("AliTRDpid::SumCharge -- ");
648 printf("ClusterArray not defined. Initialize the PID object first\n");
652 printf("AliTRDpid::SumCharge -- ");
653 printf("TRD geometry not defined. Initialize the PID object first\n");
657 // Loop through all clusters associated to this track
658 Int_t nClus = t->GetNclusters();
659 for (Int_t iClus = 0; iClus < nClus; iClus++) {
662 Int_t index = t->GetClusterIndex(iClus);
663 if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
669 if ((!cluster->From2pad()) &&
670 (!cluster->From3pad())) continue;
674 Int_t plane = fGeometry->GetPlane(cluster->GetDetector());
675 //charge[plane] += cluster->GetQ();
677 charge[plane] += t->GetClusterdQdl(iClus);
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"