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(Int_t nmom, Float_t minmom, Float_t maxmom)
265 // Creates the histograms
274 fWidMom = (maxmom - minmom) / ((Float_t) nmom);
276 fHist = new TObjArray(kNpid * nmom);
277 for (imom = 0; imom < nmom; imom++) {
278 for (ipid = 0; ipid < kNpid; ipid++) {
280 Int_t index = GetIndex(imom,ipid);
283 sprintf(name ,"hQ%03d",index);
284 if (ipid == kElectron) {
285 sprintf(title,"Q electron p-bin %03d",imom);
288 sprintf(title,"Q pion p-bin %03d",imom);
290 TH1F *hTmp = new TH1F(name,title,fNLh,fMinLh,fMaxLh);
291 fHist->AddAt(hTmp,index);
300 // //_____________________________________________________________________________
301 // Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
304 // // Fills the energy loss distributions with track <t>.
307 // Bool_t status = kTRUE;
309 // if (TMath::IsNaN(t->GetP())) return kFALSE;
311 // Float_t mom = t->GetP();
312 // Int_t ipid = MCpid(t);
313 // TH1F *hTmp = NULL;
314 // AliTRDcluster *cluster = NULL;
316 // Int_t index = GetIndex(mom,ipid);
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++) {
323 // Int_t idxClus = t->GetClusterIndex(iClus);
324 // if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(idxClus))) {
328 // hTmp->Fill(cluster->GetQ());
336 //_____________________________________________________________________________
337 Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
340 // Fills the energy loss distributions with track <t>.
343 const Int_t kNpla = AliTRDgeometry::Nplan();
345 if (TMath::IsNaN(t->GetP())) return kFALSE;
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);
353 if (!SumCharge(t,charge,nCluster)) {
359 Int_t index = GetIndex(mom,ipid);
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]));
376 //_____________________________________________________________________________
377 Int_t AliTRDpidLQ::GetIndex(const AliTRDtrack *t)
380 // Returns the histogram index
383 if (TMath::IsNaN(t->GetP())) return -1;
384 Float_t mom = t->GetP();
385 Int_t ipid = MCpid(t);
387 return GetIndex(mom,ipid);
391 //_____________________________________________________________________________
392 Int_t AliTRDpidLQ::GetIndex(Float_t mom, Int_t ipid)
395 // Returns the histogram index
398 if ((mom < fMinMom) || (mom > fMaxMom)) return -1;
399 Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
400 return GetIndex(imom,ipid);
404 //_____________________________________________________________________________
405 Int_t AliTRDpidLQ::GetIndex(Int_t imom, Int_t ipid)
408 // Returns the histogram index
411 if ((ipid < 0) || (ipid >= kNpid)) return -1;
412 return imom * kNpid + ipid;