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.4.6.1 2002/07/24 10:09:30 alibrary
21 Revision 1.5 2002/06/12 09:54:35 cblume
22 Update of tracking code provided by Sergei
24 Revision 1.4 2001/11/19 08:44:08 cblume
25 Fix bugs reported by Rene
27 Revision 1.3 2001/11/14 10:50:46 cblume
28 Changes in digits IO. Add merging of summable digits
30 Revision 1.2 2001/11/06 17:19:41 cblume
31 Add detailed geometry and simple simulator
33 Revision 1.1 2001/05/07 08:08:05 cblume
38 ///////////////////////////////////////////////////////////////////////////////
40 // The TRD particle identification base class //
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 //
46 ///////////////////////////////////////////////////////////////////////////////
53 #include <TObjArray.h>
56 #include <TParticle.h>
60 #include "AliTRDpid.h"
61 #include "AliTRDcluster.h"
62 #include "AliTRDtrack.h"
63 #include "AliTRDtracker.h"
64 #include "AliTRDgeometry.h"
68 //_____________________________________________________________________________
69 AliTRDpid::AliTRDpid():TNamed()
72 // AliTRDpid default constructor
81 fPIDpurePoints = kFALSE;
87 fThreePadOnly = kFALSE;
91 //_____________________________________________________________________________
92 AliTRDpid::AliTRDpid(const char* name, const char* title):TNamed(name,title)
95 // AliTRDpid constructor
109 //_____________________________________________________________________________
110 AliTRDpid::AliTRDpid(const AliTRDpid &p)
113 // AliTRDpid copy constructor
116 ((AliTRDpid &) p).Copy(*this);
120 //_____________________________________________________________________________
121 AliTRDpid::~AliTRDpid()
124 // AliTRDpid destructor
128 fClusterArray->Delete();
129 delete fClusterArray;
130 fClusterArray = NULL;
134 fTrackArray->Delete();
146 //_____________________________________________________________________________
147 AliTRDpid &AliTRDpid::operator=(const AliTRDpid &p)
150 // Assignment operator
153 if (this != &p) ((AliTRDpid &) p).Copy(*this);
158 //_____________________________________________________________________________
159 void AliTRDpid::Copy(TObject &p)
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;
178 //_____________________________________________________________________________
179 Bool_t AliTRDpid::Init()
182 // Initializes the PID object
185 fClusterArray = new TObjArray();
186 fTrackArray = new TObjArray();
189 fPIDpurePoints = kTRUE;
193 fThreePadOnly = kFALSE;
199 //_____________________________________________________________________________
200 Bool_t AliTRDpid::AssignLikelihood()
203 // Assigns the e / pi likelihood to all tracks.
206 return AssignLikelihood(fTrackArray);
210 //_____________________________________________________________________________
211 Bool_t AliTRDpid::AssignLikelihood(TObjArray *tarray)
214 // Assigns the e / pi likelihood to all tracks in the array <tarray>.
217 Bool_t status = kTRUE;
221 TIter nextTrack(tarray);
222 while ((track = (AliTRDtrack *) nextTrack())) {
223 if (!AssignLikelihood(track)) {
233 //_____________________________________________________________________________
234 Bool_t AliTRDpid::FillSpectra()
237 // Fills the energy loss distributions with all tracks.
240 return FillSpectra(fTrackArray);
244 //_____________________________________________________________________________
245 Bool_t AliTRDpid::FillSpectra(TObjArray *tarray)
248 // Fills the energy loss distributions with all tracks in <tarray>.
251 Bool_t status = kTRUE;
255 TIter nextTrack(tarray);
256 while ((track = (AliTRDtrack *) nextTrack())) {
257 if (!FillSpectra(track)) {
267 //_____________________________________________________________________________
268 Bool_t AliTRDpid::Open(const Char_t *name, Int_t event)
271 // Opens and reads the kine tree, the geometry, the cluster array
272 // and the track array from the file <name>.
275 Bool_t status = kTRUE;
277 status = ReadCluster(name);
278 status = ReadTracks(name);
279 status = ReadKine(name,event);
285 //_____________________________________________________________________________
286 Bool_t AliTRDpid::Open(const Char_t *namekine
287 , const Char_t *namecluster
288 , const Char_t *nametracks
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>.
297 Bool_t status = kTRUE;
299 status = ReadCluster(namecluster);
300 status = ReadTracks(nametracks);
301 status = ReadKine(namekine,event);
307 //_____________________________________________________________________________
308 Bool_t AliTRDpid::ReadKine(const Char_t *name, Int_t event)
311 // Opens and reads the kine tree and the geometry from the file <name>.
314 TFile *fFileKine = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
316 printf("AliTRDpid::ReadKine -- ");
317 printf("Open file %s\n",name);
318 fFileKine = new TFile(name);
320 printf("AliTRDpid::ReadKine -- ");
321 printf("Cannot open file %s\n",name);
326 gAlice = (AliRun *) fFileKine->Get("gAlice");
328 printf("AliTRDpid::ReadKine -- ");
329 printf("No AliRun object found\n");
332 gAlice->GetEvent(event);
334 AliTRD *trd = (AliTRD *) gAlice->GetDetector("TRD");
336 printf("AliTRDpid::ReadKine -- ");
337 printf("No TRD object found\n");
341 fGeometry = trd->GetGeometry();
343 printf("AliTRDpid::ReadKine -- ");
344 printf("No TRD geometry found\n");
354 //_____________________________________________________________________________
355 Bool_t AliTRDpid::ReadCluster(const Char_t *name)
358 // Opens and reads the cluster array from the file <name>.
361 TDirectory *savedir = gDirectory;
364 fClusterArray->Delete();
367 fClusterArray = new TObjArray();
370 printf("AliTRDpid::ReadCluster -- ");
371 printf("Open file %s\n",name);
373 AliTRDtracker *tracker = new AliTRDtracker();
374 tracker->ReadClusters(fClusterArray,name);
376 if (!fClusterArray) {
377 printf("AliTRDpid::ReadCluster -- ");
378 printf("Error reading the cluster array from file %s\n",name);
390 //_____________________________________________________________________________
391 Bool_t AliTRDpid::ReadTracks(const Char_t *name)
394 // Opens and reads the track array from the file <name>.
397 TDirectory *savedir = gDirectory;
400 fTrackArray->Delete();
403 fTrackArray = new TObjArray();
406 TFile *file = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
408 printf("AliTRDpid::ReadTracks -- ");
409 printf("Open file %s\n",name);
410 file = new TFile(name);
412 printf("AliTRDpid::ReadTracks -- ");
413 printf("Cannot open file %s\n",name);
419 sprintf(treeName,"TreeT%d_TRD",fEvent);
420 TTree *trackTree = (TTree *) file->Get(treeName);
421 TBranch *trackBranch = trackTree->GetBranch("tracks");
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);
439 //_____________________________________________________________________________
440 Int_t AliTRDpid::MCpid(const AliTRDtrack *t)
443 // Determines the pid of the MC track <t>.
447 const Int_t kPdgEl = 11;
448 const Int_t kPdgPi = 211;
450 AliTRDcluster *cluster;
453 Int_t nClusterEl = 0;
454 Int_t nClusterPi = 0;
455 Int_t nClusterPure = 0;
456 Int_t nClusterMix = 0;
463 if (!fClusterArray) {
464 printf("AliTRDpid::MCpid -- ");
465 printf("ClusterArray not defined. Initialize the PID object first\n");
469 // Loop through all clusters associated to this track
470 Int_t nCluster = t->GetNumberOfClusters();
471 for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
474 Int_t index = t->GetClusterIndex(iCluster);
475 if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
479 // Get the first two MC track indices
480 Int_t track0 = cluster->GetLabel(0);
481 Int_t track1 = cluster->GetLabel(1);
483 // Check on the track index to find the right primaries
484 if ((track0 > fPIDindexMin) &&
485 (track0 <= fPIDindexMax)) {
487 // Check whether it is a pure cluster, i.e. not overlapping
488 Bool_t accept = kTRUE;
489 if ((fPIDpurePoints) && (track1 != -1)) {
494 particle = gAlice->Particle(track0);
495 if (particle->GetFirstMother() == -1) {
496 switch (TMath::Abs(particle->GetPdgCode())) {
519 ratioEl = ((Float_t) nClusterEl) / ((Float_t) nCluster);
520 ratioPi = ((Float_t) nClusterPi) / ((Float_t) nCluster);
521 if (ratioEl > fPIDratioMin) {
524 else if (ratioPi > fPIDratioMin) {
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);
541 //_____________________________________________________________________________
542 Int_t AliTRDpid::MCpid(const AliTRDtrack *t, Int_t *pdg
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.
554 const Int_t kNtrack = 3;
555 const Int_t kNpart = 200;
557 AliTRDcluster *cluster;
565 if (!fClusterArray) {
566 printf("AliTRDpid::MCpid -- ");
567 printf("ClusterArray not defined. Initialize the PID object first\n");
572 Int_t pdgTmp[kNpart];
573 Int_t nFoundTmp[kNpart];
574 Int_t indicesTmp[kNpart];
575 for (iPart = 0; iPart < kNpart; iPart++) {
580 nFoundTmp[iPart] = 0;
581 indicesTmp[iPart] = 0;
585 // Loop through all clusters associated to this track
586 Int_t nCluster = t->GetNumberOfClusters();
587 for (iCluster = 0; iCluster < nCluster; iCluster++) {
590 Int_t index = t->GetClusterIndex(iCluster);
591 if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
595 // Get the MC track indices
596 for (iTrack = 0; iTrack < kNtrack; iTrack++) {
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]) {
610 if ((newPart) && (nPart < kNpart)) {
611 indicesTmp[nPart] = trackIndex;
612 pdgTmp[nPart] = pdgCode;
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]);
636 //_____________________________________________________________________________
637 Bool_t AliTRDpid::SumCharge(const AliTRDtrack *t
642 // Sums up the charge in each plane for track <t>.
645 Bool_t status = kTRUE;
647 AliTRDcluster *cluster;
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;
655 if (!fClusterArray) {
656 printf("AliTRDpid::SumCharge -- ");
657 printf("ClusterArray not defined. Initialize the PID object first\n");
661 printf("AliTRDpid::SumCharge -- ");
662 printf("TRD geometry not defined. Initialize the PID object first\n");
666 // Loop through all clusters associated to this track
667 Int_t nClus = t->GetNumberOfClusters();
668 for (Int_t iClus = 0; iClus < nClus; iClus++) {
671 Int_t index = t->GetClusterIndex(iClus);
672 if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
678 if ((!cluster->From2pad()) &&
679 (!cluster->From3pad())) continue;
683 Int_t plane = fGeometry->GetPlane(cluster->GetDetector());
684 //charge[plane] += cluster->GetQ();
686 charge[plane] += t->GetClusterdQdl(iClus);
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"