1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // The TRD particle identification class //
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 //
27 ///////////////////////////////////////////////////////////////////////////////
34 #include <TObjArray.h>
37 #include <TParticle.h>
41 #include "AliTRDpidLQ.h"
42 #include "AliTRDcluster.h"
43 #include "AliTRDtrack.h"
44 #include "AliTRDtracker.h"
45 #include "AliTRDgeometry.h"
49 //_____________________________________________________________________________
50 AliTRDpidLQ::AliTRDpidLQ():AliTRDpid()
53 // AliTRDpidLQ default constructor
72 //_____________________________________________________________________________
73 AliTRDpidLQ::AliTRDpidLQ(const char* name, const char* title)
74 :AliTRDpid(name,title)
77 // AliTRDpidLQ constructor
95 //_____________________________________________________________________________
96 AliTRDpidLQ::AliTRDpidLQ(const AliTRDpidLQ &p):AliTRDpid(p)
99 // AliTRDpidLQ copy constructor
102 ((AliTRDpidLQ &) p).Copy(*this);
106 //_____________________________________________________________________________
107 AliTRDpidLQ::~AliTRDpidLQ()
110 // AliTRDpidLQ destructor
120 //_____________________________________________________________________________
121 AliTRDpidLQ &AliTRDpidLQ::operator=(const AliTRDpidLQ &p)
124 // Assignment operator
127 if (this != &p) ((AliTRDpidLQ &) p).Copy(*this);
132 //_____________________________________________________________________________
133 void AliTRDpidLQ::Copy(TObject &p)
139 fHist->Copy(*((AliTRDpidLQ &) p).fHist);
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;
155 //_____________________________________________________________________________
156 Bool_t AliTRDpidLQ::Init()
159 // Initializes the PID object
173 //_____________________________________________________________________________
174 Bool_t AliTRDpidLQ::AssignLikelihood(AliTRDtrack *t)
177 // Assigns the e / pi likelihood to a given track
180 const Int_t kNpla = AliTRDgeometry::Nplan();
181 Float_t * charge = new Float_t[kNpla];
182 Int_t * nCluster = new Int_t[kNpla];
196 t->SetLikelihoodElectron(-1.);
197 if (TMath::IsNaN(t->GetP())) return kFALSE;
198 Float_t mom = t->GetP();
200 // Calculate the total charge in each plane
201 if (!SumCharge(t,charge,nCluster)) return kFALSE;
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);
218 pEl = pEl * pElPlane;
220 Int_t ibinPi = hTmpPi->FindBin(chargeNorm);
221 Float_t pPiPlane = hTmpPi->GetBinContent(ibinPi);
223 pPi = pPi * pPiPlane;
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);
250 // printf("---- mom = %f, lhEl = %f, lhPi = %f\n",mom,lhEl,lhPi);
252 // Assign the likelihoods
253 t->SetLikelihoodElectron(lhEl);
261 //_____________________________________________________________________________
262 Bool_t AliTRDpidLQ::CreateHistograms(const Int_t nmom
263 , const Float_t minmom
264 , const Float_t maxmom)
267 // Creates the histograms
276 fWidMom = (maxmom - minmom) / ((Float_t) nmom);
278 fHist = new TObjArray(kNpid * nmom);
279 for (imom = 0; imom < nmom; imom++) {
280 for (ipid = 0; ipid < kNpid; ipid++) {
282 Int_t index = GetIndex(imom,ipid);
285 sprintf(name ,"hQ%03d",index);
286 if (ipid == kElectron) {
287 sprintf(title,"Q electron p-bin %03d",imom);
290 sprintf(title,"Q pion p-bin %03d",imom);
292 TH1F *hTmp = new TH1F(name,title,fNLh,fMinLh,fMaxLh);
293 fHist->AddAt(hTmp,index);
302 // //_____________________________________________________________________________
303 // Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
306 // // Fills the energy loss distributions with track <t>.
309 // Bool_t status = kTRUE;
311 // if (TMath::IsNaN(t->GetP())) return kFALSE;
313 // Float_t mom = t->GetP();
314 // Int_t ipid = MCpid(t);
315 // TH1F *hTmp = NULL;
316 // AliTRDcluster *cluster = NULL;
318 // Int_t index = GetIndex(mom,ipid);
320 // hTmp = (TH1F *) fHist->UncheckedAt(index);
321 // // Loop through all clusters associated to this track
322 // Int_t nClus = t->GetNclusters();
323 // for (Int_t iClus = 0; iClus < nClus; iClus++) {
325 // Int_t idxClus = t->GetClusterIndex(iClus);
326 // if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(idxClus))) {
330 // hTmp->Fill(cluster->GetQ());
338 //_____________________________________________________________________________
339 Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
342 // Fills the energy loss distributions with track <t>.
345 const Int_t kNpla = AliTRDgeometry::Nplan();
347 if (TMath::IsNaN(t->GetP())) return kFALSE;
349 Float_t * charge = new Float_t[kNpla];
350 Int_t * nCluster = new Int_t[kNpla];
351 Float_t mom = t->GetP();
352 Int_t ipid = MCpid(t);
355 if (!SumCharge(t,charge,nCluster)) {
361 Int_t index = GetIndex(mom,ipid);
363 hTmp = (TH1F *) fHist->UncheckedAt(index);
364 for (Int_t ipla = 0; ipla < kNpla; ipla++) {
365 if ((charge[ipla] > fChargeMin) &&
366 (nCluster[ipla] > fNClusterMin)){
367 hTmp->Fill(charge[ipla] / ((Float_t) nCluster[ipla]));
378 //_____________________________________________________________________________
379 Int_t AliTRDpidLQ::GetIndex(const AliTRDtrack *t)
382 // Returns the histogram index
385 if (TMath::IsNaN(t->GetP())) return -1;
386 Float_t mom = t->GetP();
387 Int_t ipid = MCpid(t);
389 return GetIndex(mom,ipid);
393 //_____________________________________________________________________________
394 Int_t AliTRDpidLQ::GetIndex(const Float_t mom, const Int_t ipid)
397 // Returns the histogram index
400 if ((mom < fMinMom) || (mom > fMaxMom)) return -1;
401 Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
402 return GetIndex(imom,ipid);
406 //_____________________________________________________________________________
407 Int_t AliTRDpidLQ::GetIndex(const Int_t imom, const Int_t ipid)
410 // Returns the histogram index
413 if ((ipid < 0) || (ipid >= kNpid)) return -1;
414 return imom * kNpid + ipid;