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