TRD/AliTRDpid.cxx

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