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