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