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