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