TRD/AliTRDpidLQ.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 class                                   //
  21 //                                                                           //
  22 //   Its main purposes are:                                                  //
  23 //      - Creation and bookkeeping of the propability distributions          //
  24 //      - Assignment of a e/pi propability to a given track based on         //
  25 //        the LQ method                                                      //
  26 //                                                                           //
  27 ///////////////////////////////////////////////////////////////////////////////
  28
  29 #include <stdlib.h>
  30 #include <math.h>
  31
  32 #include <TROOT.h>
  33 #include <TH1.h>
  34 #include <TObjArray.h>
  35 #include <TTree.h>
  36 #include <TFile.h>
  37 #include <TParticle.h>
  38
  39 #include "AliRun.h"
  40 #include "AliTRDpidLQ.h"
  41 #include "AliTRDcluster.h"
  42 #include "AliTRDtrack.h"
  43 #include "AliTRDtracker.h"
  44 #include "AliTRDgeometry.h"
  45
  46 ClassImp(AliTRDpidLQ)
  47
  48 //_____________________________________________________________________________
  49 AliTRDpidLQ::AliTRDpidLQ():AliTRDpid()
  50 {
  51   //
  52   // AliTRDpidLQ default constructor
  53   //
  54
  55   fNMom   = 0;
  56   fMinMom = 0;
  57   fMaxMom = 0;
  58   fWidMom = 0;
  59
  60   fNLh    = 0;
  61   fMinLh  = 0;
  62   fMaxLh  = 0;
  63
  64   fHist   = NULL;
  65
  66   fChargeMin   = 0;
  67   fNClusterMin = 0;
  68
  69 }
  70
  71 //_____________________________________________________________________________
  72 AliTRDpidLQ::AliTRDpidLQ(const char* name, const char* title)
  73             :AliTRDpid(name,title)
  74 {
  75   //
  76   // AliTRDpidLQ constructor
  77   //
  78
  79   fNMom   = 0;
  80   fMinMom = 0;
  81   fMaxMom = 0;
  82   fWidMom = 0;
  83
  84   fNLh    = 0;
  85   fMinLh  = 0;
  86   fMaxLh  = 0;
  87
  88   fHist   = NULL;
  89
  90   Init();
  91
  92 }
  93
  94 //_____________________________________________________________________________
  95 AliTRDpidLQ::AliTRDpidLQ(const AliTRDpidLQ &p):AliTRDpid(p)
  96 {
  97   //
  98   // AliTRDpidLQ copy constructor
  99   //
 100
 101   ((AliTRDpidLQ &) p).Copy(*this);
 102
 103 }
 104
 105 //_____________________________________________________________________________
 106 AliTRDpidLQ::~AliTRDpidLQ()
 107 {
 108   //
 109   // AliTRDpidLQ destructor
 110   //
 111
 112   if (fHist) {
 113     fHist->Delete();
 114     delete fHist;
 115   }
 116
 117 }
 118
 119 //_____________________________________________________________________________
 120 AliTRDpidLQ &AliTRDpidLQ::operator=(const AliTRDpidLQ &p)
 121 {
 122   //
 123   // Assignment operator
 124   //
 125
 126   if (this != &p) ((AliTRDpidLQ &) p).Copy(*this);
 127   return *this;
 128
 129 }
 130
 131 //_____________________________________________________________________________
 132 void AliTRDpidLQ::Copy(TObject &p) const
 133 {
 134   //
 135   // Copy function
 136   //
 137
 138   fHist->Copy(*((AliTRDpidLQ &) p).fHist);
 139
 140   ((AliTRDpidLQ &) p).fNMom        = fNMom;
 141   ((AliTRDpidLQ &) p).fMinMom      = fMinMom;
 142   ((AliTRDpidLQ &) p).fMaxMom      = fMaxMom;
 143   ((AliTRDpidLQ &) p).fWidMom      = fWidMom;
 144   ((AliTRDpidLQ &) p).fChargeMin   = fChargeMin;
 145   ((AliTRDpidLQ &) p).fNClusterMin = fNClusterMin;
 146   ((AliTRDpidLQ &) p).fNLh         = fNLh;
 147   ((AliTRDpidLQ &) p).fMinLh       = fMinLh;
 148   ((AliTRDpidLQ &) p).fMaxLh       = fMaxLh;
 149
 150   AliTRDpid::Copy(p);
 151
 152 }
 153
 154 //_____________________________________________________________________________
 155 Bool_t AliTRDpidLQ::Init()
 156 {
 157   //
 158   // Initializes the PID object
 159   //
 160
 161   fChargeMin   = 0.0;
 162   fNClusterMin = 10;
 163
 164   fNLh         = 100;
 165   fMinLh       =   0.0;
 166   fMaxLh       = 500.0;
 167
 168   return kTRUE;
 169
 170 }
 171
 172 //_____________________________________________________________________________
 173 Bool_t AliTRDpidLQ::AssignLikelihood(AliTRDtrack *t)
 174 {
 175   //
 176   // Assigns the e / pi likelihood to a given track
 177   //
 178
 179   const Int_t kNpla = AliTRDgeometry::Nplan();
 180   Float_t * charge = new Float_t[kNpla];
 181   Int_t   * nCluster = new Int_t[kNpla];
 182
 183   Float_t  lhPi = 0;
 184   Float_t  lhEl = 0;
 185
 186   Double_t pPi  = 1;
 187   Double_t pEl  = 1;
 188   Double_t pSum = 0;
 189
 190   Int_t   indexEl;
 191   Int_t   indexPi;
 192   TH1F   *hTmpEl;
 193   TH1F   *hTmpPi;
 194
 195   t->SetLikelihoodElectron(-1.);
 196   if (TMath::IsNaN(t->GetP())) return kFALSE;
 197   Float_t mom = t->GetP();
 198
 199   // Calculate the total charge in each plane
 200   if (!SumCharge(t,charge,nCluster)) return kFALSE;
 201
 202   indexEl = GetIndex(mom,kElectron);
 203   indexPi = GetIndex(mom,kPion);
 204   if ((indexEl > -1) && (indexPi > -1)) {
 205     hTmpEl = (TH1F *) fHist->UncheckedAt(indexEl);
 206     hTmpPi = (TH1F *) fHist->UncheckedAt(indexPi);
 207     for (Int_t ipla = 0; ipla < kNpla; ipla++) {
 208       Float_t chargePlane   = charge[ipla];
 209       Int_t   nClusterPlane = nCluster[ipla];
 210       if ((chargePlane   >   fChargeMin) &&
 211           (nClusterPlane > fNClusterMin)){
 212         Float_t chargeNorm = chargePlane / ((Float_t) nClusterPlane);
 213         if (chargeNorm < fMaxLh) {
 214           Int_t   ibinEl     = hTmpEl->FindBin(chargeNorm);
 215           Float_t pElPlane   = hTmpEl->GetBinContent(ibinEl);
 216           if (pElPlane > 0) {
 217             pEl  = pEl * pElPlane;
 218           }
 219           Int_t   ibinPi     = hTmpPi->FindBin(chargeNorm);
 220           Float_t pPiPlane   = hTmpPi->GetBinContent(ibinPi);
 221           if (pPiPlane > 0) {
 222             pPi  = pPi * pPiPlane;
 223           }
 224 //           printf("          ipla = %d, nCluster = %d, charge = %f\n"
 225 //                  ,ipla,nClusterPlane,chargeNorm);
 226 //           printf("electron: pElPlane = %f, ibinEl = %d, pEl = %f\n"
 227 //               ,pElPlane,ibinEl,pEl);
 228 //           printf("    pion: pPiPlane = %f, ibinPi = %d, pPi = %f\n"
 229 //               ,pPiPlane,ibinPi,pPi);
 230         }
 231       }
 232     }
 233   }
 234   else {
 235     delete [] charge;
 236     delete [] nCluster;
 237     return kTRUE;
 238   }
 239
 240   pSum = pEl + pPi;
 241   if (pSum <= 0) {
 242     delete [] charge;
 243     delete [] nCluster;
 244     return kFALSE;
 245   }
 246   lhEl = pEl / pSum;
 247   lhPi = pPi / pSum;
 248
 249 //   printf("---- mom = %f, lhEl = %f, lhPi = %f\n",mom,lhEl,lhPi);
 250
 251   // Assign the likelihoods
 252   t->SetLikelihoodElectron(lhEl);
 253
 254   delete [] charge;
 255   delete [] nCluster;
 256   return kTRUE;
 257
 258 }
 259
 260 //_____________________________________________________________________________
 261 Bool_t AliTRDpidLQ::CreateHistograms(Int_t   nmom, Float_t minmom, Float_t maxmom)
 262 {
 263   //
 264   // Creates the histograms
 265   //
 266
 267   Int_t imom;
 268   Int_t ipid;
 269
 270   fNMom   = nmom;
 271   fMinMom = minmom;
 272   fMaxMom = maxmom;
 273   fWidMom = (maxmom - minmom) / ((Float_t) nmom);
 274
 275   fHist = new TObjArray(kNpid * nmom);
 276   for (imom = 0; imom < nmom;  imom++) {
 277     for (ipid = 0; ipid < kNpid; ipid++) {
 278
 279       Int_t index = GetIndex(imom,ipid);
 280       Char_t name[10];
 281       Char_t title[80];
 282       sprintf(name ,"hQ%03d",index);
 283       if (ipid == kElectron) {
 284         sprintf(title,"Q electron p-bin %03d",imom);
 285       }
 286       else {
 287         sprintf(title,"Q pion p-bin %03d",imom);
 288       }
 289       TH1F *hTmp = new TH1F(name,title,fNLh,fMinLh,fMaxLh);
 290       fHist->AddAt(hTmp,index);
 291
 292     }
 293   }
 294
 295   return kTRUE;
 296
 297 }
 298
 299 // //_____________________________________________________________________________
 300 // Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
 301 // {
 302 //   //
 303 //   // Fills the energy loss distributions with track <t>.
 304 //   //
 305
 306 //   Bool_t status = kTRUE;
 307
 308 //   if (TMath::IsNaN(t->GetP())) return kFALSE;
 309
 310 //   Float_t        mom     = t->GetP();
 311 //   Int_t          ipid    = MCpid(t);
 312 //   TH1F          *hTmp    = NULL;
 313 //   AliTRDcluster *cluster = NULL;
 314
 315 //   Int_t index = GetIndex(mom,ipid);
 316 //   if (index > -1) {
 317 //     hTmp = (TH1F *) fHist->UncheckedAt(index);
 318 //     // Loop through all clusters associated to this track
 319 //     Int_t nClus = t->GetNclusters();
 320 //     for (Int_t iClus = 0; iClus < nClus; iClus++) {
 321 //       // Get a cluster
 322 //       Int_t idxClus = t->GetClusterIndex(iClus);
 323 //       if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(idxClus))) {
 324 //         status = kFALSE;
 325 //         break;
 326 //       }
 327 //       hTmp->Fill(cluster->GetQ());
 328 //     }
 329 //   }
 330
 331 //   return status;
 332
 333 // }
 334
 335 //_____________________________________________________________________________
 336 Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
 337 {
 338   //
 339   // Fills the energy loss distributions with track <t>.
 340   //
 341
 342   const Int_t kNpla = AliTRDgeometry::Nplan();
 343
 344   if (TMath::IsNaN(t->GetP())) return kFALSE;
 345
 346   Float_t * charge = new Float_t[kNpla];
 347   Int_t   * nCluster = new Int_t[kNpla];
 348   Float_t mom  = t->GetP();
 349   Int_t   ipid = MCpid(t);
 350   TH1F   *hTmp = NULL;
 351
 352   if (!SumCharge(t,charge,nCluster)) {
 353     delete [] charge;
 354     delete [] nCluster;
 355     return kFALSE;
 356   }
 357
 358   Int_t index = GetIndex(mom,ipid);
 359   if (index > -1) {
 360     hTmp = (TH1F *) fHist->UncheckedAt(index);
 361     for (Int_t ipla = 0; ipla < kNpla; ipla++) {
 362       if ((charge[ipla]   >   fChargeMin) &&
 363           (nCluster[ipla] > fNClusterMin)){
 364         hTmp->Fill(charge[ipla] / ((Float_t) nCluster[ipla]));
 365       }
 366     }
 367   }
 368
 369   delete [] charge;
 370   delete [] nCluster;
 371   return kTRUE;
 372
 373 }
 374
 375 //_____________________________________________________________________________
 376 Int_t AliTRDpidLQ::GetIndex(const AliTRDtrack *t)
 377 {
 378   //
 379   // Returns the histogram index
 380   //
 381
 382   if (TMath::IsNaN(t->GetP())) return -1;
 383   Float_t mom  = t->GetP();
 384   Int_t   ipid = MCpid(t);
 385
 386   return GetIndex(mom,ipid);
 387
 388 }
 389
 390 //_____________________________________________________________________________
 391 Int_t AliTRDpidLQ::GetIndex(Float_t mom, Int_t ipid)
 392 {
 393   //
 394   // Returns the histogram index
 395   //
 396
 397   if ((mom < fMinMom) || (mom > fMaxMom))  return -1;
 398   Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
 399   return GetIndex(imom,ipid);
 400
 401 }
 402
 403 //_____________________________________________________________________________
 404 Int_t AliTRDpidLQ::GetIndex(Int_t imom, Int_t ipid)
 405 {
 406   //
 407   // Returns the histogram index
 408   //
 409
 410   if ((ipid < 0) || (ipid >= kNpid)) return -1;
 411   return imom * kNpid + ipid;
 412
 413 }