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