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