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.3  2001/11/14 10:50:46  cblume
  19 Changes in digits IO. Add merging of summable digits
  20
  21 Revision 1.2  2001/11/06 17:19:41  cblume
  22 Add detailed geometry and simple simulator
  23
  24 Revision 1.1  2001/05/07 08:08:05  cblume
  25 Update of TRD code
  26
  27 */
  28
  29 ///////////////////////////////////////////////////////////////////////////////
  30 //                                                                           //
  31 //   The TRD particle identification base class                              //
  32 //                                                                           //
  33 //   Its main purposes are:                                                  //
  34 //      - Provide I/O framework for all neccessary files                     //
  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
  57 ClassImp(AliTRDpid)
  58
  59 //_____________________________________________________________________________
  60 AliTRDpid::AliTRDpid():TNamed()
  61 {
  62   //
  63   // AliTRDpid default constructor
  64   //
  65
  66   fTrackArray    = NULL;
  67   fClusterArray  = NULL;
  68   fGeometry      = NULL;
  69   fFileKine      = NULL;
  70
  71   fPIDratioMin   = 0.0;
  72   fPIDpurePoints = kFALSE;
  73   fPIDindexMin   = 0;
  74   fPIDindexMax   = 0;
  75
  76   fEvent         = 0;
  77
  78   fThreePadOnly  = kFALSE;
  79
  80 }
  81
  82 //_____________________________________________________________________________
  83 AliTRDpid::AliTRDpid(const char* name, const char* title):TNamed(name,title)
  84 {
  85   //
  86   // AliTRDpid constructor
  87   //
  88
  89   fTrackArray   = NULL;
  90   fClusterArray = NULL;
  91   fGeometry     = NULL;
  92   fFileKine     = NULL;
  93
  94   fEvent        = 0;
  95
  96   Init();
  97
  98 }
  99
 100 //_____________________________________________________________________________
 101 AliTRDpid::AliTRDpid(const AliTRDpid &p)
 102 {
 103   //
 104   // AliTRDpid copy constructor
 105   //
 106
 107   ((AliTRDpid &) p).Copy(*this);
 108
 109 }
 110
 111 //_____________________________________________________________________________
 112 AliTRDpid::~AliTRDpid()
 113 {
 114   //
 115   // AliTRDpid destructor
 116   //
 117
 118   if (fClusterArray) {
 119     fClusterArray->Delete();
 120     delete fClusterArray;
 121     fClusterArray = NULL;
 122   }
 123
 124   if (fTrackArray) {
 125     fTrackArray->Delete();
 126     delete fTrackArray;
 127     fTrackArray = NULL;
 128   }
 129
 130   if (fFileKine) {
 131     delete fFileKine;
 132     fFileKine = NULL;
 133   }
 134
 135 }
 136
 137 //_____________________________________________________________________________
 138 AliTRDpid &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 //_____________________________________________________________________________
 150 void AliTRDpid::Copy(TObject &p)
 151 {
 152   //
 153   // Copy function
 154   //
 155
 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;
 165   ((AliTRDpid &) p).fEvent         = fEvent;
 166
 167 }
 168
 169 //_____________________________________________________________________________
 170 Bool_t AliTRDpid::Init()
 171 {
 172   //
 173   // Initializes the PID object
 174   //
 175
 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;
 185
 186   return kTRUE;
 187
 188 }
 189
 190 //_____________________________________________________________________________
 191 Bool_t AliTRDpid::AssignLikelihood()
 192 {
 193   //
 194   // Assigns the e / pi likelihood to all tracks.
 195   //
 196
 197   return AssignLikelihood(fTrackArray);
 198
 199 }
 200
 201 //_____________________________________________________________________________
 202 Bool_t AliTRDpid::AssignLikelihood(TObjArray *tarray)
 203 {
 204   //
 205   // Assigns the e / pi likelihood to all tracks in the array <tarray>.
 206   //
 207
 208   Bool_t status = kTRUE;
 209
 210   AliTRDtrack *track;
 211
 212   TIter nextTrack(tarray);
 213   while ((track = (AliTRDtrack *) nextTrack())) {
 214     if (!AssignLikelihood(track)) {
 215       status = kFALSE;
 216       continue;
 217     }
 218   }
 219
 220   return status;
 221
 222 }
 223
 224 //_____________________________________________________________________________
 225 Bool_t AliTRDpid::FillSpectra()
 226 {
 227   //
 228   // Fills the energy loss distributions with all tracks.
 229   //
 230
 231   return FillSpectra(fTrackArray);
 232
 233 }
 234
 235 //_____________________________________________________________________________
 236 Bool_t AliTRDpid::FillSpectra(TObjArray *tarray)
 237 {
 238   //
 239   // Fills the energy loss distributions with all tracks in <tarray>.
 240   //
 241
 242   Bool_t status = kTRUE;
 243
 244   AliTRDtrack *track;
 245
 246   TIter nextTrack(tarray);
 247   while ((track = (AliTRDtrack *) nextTrack())) {
 248     if (!FillSpectra(track)) {
 249       status = kFALSE;
 250       continue;
 251     }
 252   }
 253
 254   return kTRUE;
 255
 256 }
 257
 258 //_____________________________________________________________________________
 259 Bool_t AliTRDpid::Open(const Char_t *name, Int_t event)
 260 {
 261   //
 262   // Opens and reads the kine tree, the geometry, the cluster array
 263   // and the track array from the file <name>.
 264   //
 265
 266   Bool_t status = kTRUE;
 267
 268   status = ReadCluster(name);
 269   status = ReadTracks(name);
 270   status = ReadKine(name,event);
 271
 272   return status;
 273
 274 }
 275
 276 //_____________________________________________________________________________
 277 Bool_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>,
 285   // and the track array from the file <nametracks>.
 286   //
 287
 288   Bool_t status = kTRUE;
 289
 290   status = ReadCluster(namecluster);
 291   status = ReadTracks(nametracks);
 292   status = ReadKine(namekine,event);
 293
 294   return status;
 295
 296 }
 297
 298 //_____________________________________________________________________________
 299 Bool_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
 305   TFile *fFileKine = (TFile *) gROOT->GetListOfFiles()->FindObject(name);
 306   if (!fFileKine) {
 307     printf("AliTRDpid::ReadKine -- ");
 308     printf("Open file %s\n",name);
 309     fFileKine = new TFile(name);
 310     if (!fFileKine) {
 311       printf("AliTRDpid::ReadKine -- ");
 312       printf("Cannot open file %s\n",name);
 313       return kFALSE;
 314     }
 315   }
 316
 317   gAlice = (AliRun *) fFileKine->Get("gAlice");
 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
 339   fEvent = event;
 340
 341   return kTRUE;
 342
 343 }
 344
 345 //_____________________________________________________________________________
 346 Bool_t AliTRDpid::ReadCluster(const Char_t *name)
 347 {
 348   //
 349   // Opens and reads the cluster array from the file <name>.
 350   //
 351
 352   TDirectory *savedir = gDirectory;
 353
 354   if (fClusterArray) {
 355     fClusterArray->Delete();
 356   }
 357   else {
 358     fClusterArray = new TObjArray();
 359   }
 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
 375   savedir->cd();
 376
 377   return kTRUE;
 378
 379 }
 380
 381 //_____________________________________________________________________________
 382 Bool_t AliTRDpid::ReadTracks(const Char_t *name)
 383 {
 384   //
 385   // Opens and reads the track array from the file <name>.
 386   //
 387
 388   TDirectory *savedir = gDirectory;
 389
 390   if (fTrackArray) {
 391     fTrackArray->Delete();
 392   }
 393   else {
 394     fTrackArray = new TObjArray();
 395   }
 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
 409   Char_t treeName[12];
 410   sprintf(treeName,"TreeT%d_TRD",fEvent);
 411   TTree   *trackTree   = (TTree *) file->Get(treeName);
 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
 424   savedir->cd();
 425
 426   return kTRUE;
 427
 428 }
 429
 430 //_____________________________________________________________________________
 431 Int_t AliTRDpid::MCpid(const AliTRDtrack *t)
 432 {
 433   //
 434   // Determines the pid of the MC track <t>.
 435   //
 436
 437   // PDG codes
 438   const Int_t kPdgEl =  11;
 439   const Int_t kPdgPi = 211;
 440
 441   AliTRDcluster *cluster;
 442   TParticle     *particle;
 443
 444   Int_t   nClusterEl   = 0;
 445   Int_t   nClusterPi   = 0;
 446   Int_t   nClusterPure = 0;
 447   Int_t   nClusterMix  = 0;
 448
 449   Float_t ratioEl      = 0;
 450   Float_t ratioPi      = 0;
 451
 452   Int_t   ipid         = -1;
 453
 454   if (!fClusterArray) {
 455     printf("AliTRDpid::MCpid -- ");
 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);
 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       }
 504
 505     }
 506
 507   }
 508
 509   if (nCluster) {
 510     ratioEl = ((Float_t) nClusterEl) / ((Float_t) nCluster);
 511     ratioPi = ((Float_t) nClusterPi) / ((Float_t) nCluster);
 512     if      (ratioEl > fPIDratioMin) {
 513       ipid = kElectron;
 514     }
 515     else if (ratioPi > fPIDratioMin) {
 516       ipid = kPion;
 517     }
 518   }
 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);
 527
 528   return ipid;
 529
 530 }
 531
 532 //_____________________________________________________________________________
 533 Int_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 //_____________________________________________________________________________
 628 Bool_t AliTRDpid::SumCharge(const AliTRDtrack *t
 629                            , Float_t *charge
 630                            , Int_t   *nCluster)
 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++) {
 642     charge[iPlane]   = 0.0;
 643     nCluster[iPlane] = 0;
 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
 658   Int_t nClus = t->GetNclusters();
 659   for (Int_t iClus = 0; iClus < nClus; iClus++) {
 660
 661     // Get a cluster
 662     Int_t index = t->GetClusterIndex(iClus);
 663     if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(index))) {
 664       status = kFALSE;
 665       break;
 666     }
 667
 668     if (fThreePadOnly) {
 669       if ((!cluster->From2pad()) &&
 670           (!cluster->From3pad())) continue;
 671     }
 672
 673     // Sum up the charge
 674     Int_t plane    = fGeometry->GetPlane(cluster->GetDetector());
 675     //charge[plane] += cluster->GetQ();
 676     // Corrected charge
 677     charge[plane] += t->GetClusterdQdl(iClus);
 678     nCluster[plane]++;
 679
 680   }
 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]);
 697
 698   return status;
 699
 700 }