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>
40 #include "AliTRDpidLQ.h"
41 #include "AliTRDcluster.h"
42 #include "AliTRDtrack.h"
43 #include "AliTRDtracker.h"
44 #include "AliTRDgeometry.h"
48 //_____________________________________________________________________________
49 AliTRDpidLQ::AliTRDpidLQ():AliTRDpid()
52 // AliTRDpidLQ default constructor
71 //_____________________________________________________________________________
72 AliTRDpidLQ::AliTRDpidLQ(const char* name, const char* title)
73 :AliTRDpid(name,title)
76 // AliTRDpidLQ constructor
94 //_____________________________________________________________________________
95 AliTRDpidLQ::AliTRDpidLQ(const AliTRDpidLQ &p):AliTRDpid(p)
98 // AliTRDpidLQ copy constructor
101 ((AliTRDpidLQ &) p).Copy(*this);
105 //_____________________________________________________________________________
106 AliTRDpidLQ::~AliTRDpidLQ()
109 // AliTRDpidLQ destructor
119 //_____________________________________________________________________________
120 AliTRDpidLQ &AliTRDpidLQ::operator=(const AliTRDpidLQ &p)
123 // Assignment operator
126 if (this != &p) ((AliTRDpidLQ &) p).Copy(*this);
131 //_____________________________________________________________________________
132 void AliTRDpidLQ::Copy(TObject &p) const
138 fHist->Copy(*((AliTRDpidLQ &) p).fHist);
140 ((AliTRDpidLQ &) p).fNMom = fNMom;
141 ((AliTRDpidLQ &) p).fMinMom = fMinMom;
142 ((AliTRDpidLQ &) p).fMaxMom = fMaxMom;
143 ((AliTRDpidLQ &) p).fWidMom = fWidMom;
144 ((AliTRDpidLQ &) p).fChargeMin = fChargeMin;
145 ((AliTRDpidLQ &) p).fNClusterMin = fNClusterMin;
146 ((AliTRDpidLQ &) p).fNLh = fNLh;
147 ((AliTRDpidLQ &) p).fMinLh = fMinLh;
148 ((AliTRDpidLQ &) p).fMaxLh = fMaxLh;
154 //_____________________________________________________________________________
155 Bool_t AliTRDpidLQ::Init()
158 // Initializes the PID object
172 //_____________________________________________________________________________
173 Bool_t AliTRDpidLQ::AssignLikelihood(AliTRDtrack *t)
176 // Assigns the e / pi likelihood to a given track
179 const Int_t kNpla = AliTRDgeometry::Nplan();
180 Float_t * charge = new Float_t[kNpla];
181 Int_t * nCluster = new Int_t[kNpla];
195 t->SetLikelihoodElectron(-1.);
196 if (TMath::IsNaN(t->GetP())) return kFALSE;
197 Float_t mom = t->GetP();
199 // Calculate the total charge in each plane
200 if (!SumCharge(t,charge,nCluster)) return kFALSE;
202 indexEl = GetIndex(mom,kElectron);
203 indexPi = GetIndex(mom,kPion);
204 if ((indexEl > -1) && (indexPi > -1)) {
205 hTmpEl = (TH1F *) fHist->UncheckedAt(indexEl);
206 hTmpPi = (TH1F *) fHist->UncheckedAt(indexPi);
207 for (Int_t ipla = 0; ipla < kNpla; ipla++) {
208 Float_t chargePlane = charge[ipla];
209 Int_t nClusterPlane = nCluster[ipla];
210 if ((chargePlane > fChargeMin) &&
211 (nClusterPlane > fNClusterMin)){
212 Float_t chargeNorm = chargePlane / ((Float_t) nClusterPlane);
213 if (chargeNorm < fMaxLh) {
214 Int_t ibinEl = hTmpEl->FindBin(chargeNorm);
215 Float_t pElPlane = hTmpEl->GetBinContent(ibinEl);
217 pEl = pEl * pElPlane;
219 Int_t ibinPi = hTmpPi->FindBin(chargeNorm);
220 Float_t pPiPlane = hTmpPi->GetBinContent(ibinPi);
222 pPi = pPi * pPiPlane;
224 // printf(" ipla = %d, nCluster = %d, charge = %f\n"
225 // ,ipla,nClusterPlane,chargeNorm);
226 // printf("electron: pElPlane = %f, ibinEl = %d, pEl = %f\n"
227 // ,pElPlane,ibinEl,pEl);
228 // printf(" pion: pPiPlane = %f, ibinPi = %d, pPi = %f\n"
229 // ,pPiPlane,ibinPi,pPi);
249 // printf("---- mom = %f, lhEl = %f, lhPi = %f\n",mom,lhEl,lhPi);
251 // Assign the likelihoods
252 t->SetLikelihoodElectron(lhEl);
260 //_____________________________________________________________________________
261 Bool_t AliTRDpidLQ::CreateHistograms(Int_t nmom, Float_t minmom, Float_t maxmom)
264 // Creates the histograms
273 fWidMom = (maxmom - minmom) / ((Float_t) nmom);
275 fHist = new TObjArray(kNpid * nmom);
276 for (imom = 0; imom < nmom; imom++) {
277 for (ipid = 0; ipid < kNpid; ipid++) {
279 Int_t index = GetIndex(imom,ipid);
282 sprintf(name ,"hQ%03d",index);
283 if (ipid == kElectron) {
284 sprintf(title,"Q electron p-bin %03d",imom);
287 sprintf(title,"Q pion p-bin %03d",imom);
289 TH1F *hTmp = new TH1F(name,title,fNLh,fMinLh,fMaxLh);
290 fHist->AddAt(hTmp,index);
299 // //_____________________________________________________________________________
300 // Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
303 // // Fills the energy loss distributions with track <t>.
306 // Bool_t status = kTRUE;
308 // if (TMath::IsNaN(t->GetP())) return kFALSE;
310 // Float_t mom = t->GetP();
311 // Int_t ipid = MCpid(t);
312 // TH1F *hTmp = NULL;
313 // AliTRDcluster *cluster = NULL;
315 // Int_t index = GetIndex(mom,ipid);
317 // hTmp = (TH1F *) fHist->UncheckedAt(index);
318 // // Loop through all clusters associated to this track
319 // Int_t nClus = t->GetNclusters();
320 // for (Int_t iClus = 0; iClus < nClus; iClus++) {
322 // Int_t idxClus = t->GetClusterIndex(iClus);
323 // if (!(cluster = (AliTRDcluster *) fClusterArray->UncheckedAt(idxClus))) {
327 // hTmp->Fill(cluster->GetQ());
335 //_____________________________________________________________________________
336 Bool_t AliTRDpidLQ::FillSpectra(const AliTRDtrack *t)
339 // Fills the energy loss distributions with track <t>.
342 const Int_t kNpla = AliTRDgeometry::Nplan();
344 if (TMath::IsNaN(t->GetP())) return kFALSE;
346 Float_t * charge = new Float_t[kNpla];
347 Int_t * nCluster = new Int_t[kNpla];
348 Float_t mom = t->GetP();
349 Int_t ipid = MCpid(t);
352 if (!SumCharge(t,charge,nCluster)) {
358 Int_t index = GetIndex(mom,ipid);
360 hTmp = (TH1F *) fHist->UncheckedAt(index);
361 for (Int_t ipla = 0; ipla < kNpla; ipla++) {
362 if ((charge[ipla] > fChargeMin) &&
363 (nCluster[ipla] > fNClusterMin)){
364 hTmp->Fill(charge[ipla] / ((Float_t) nCluster[ipla]));
375 //_____________________________________________________________________________
376 Int_t AliTRDpidLQ::GetIndex(const AliTRDtrack *t)
379 // Returns the histogram index
382 if (TMath::IsNaN(t->GetP())) return -1;
383 Float_t mom = t->GetP();
384 Int_t ipid = MCpid(t);
386 return GetIndex(mom,ipid);
390 //_____________________________________________________________________________
391 Int_t AliTRDpidLQ::GetIndex(Float_t mom, Int_t ipid)
394 // Returns the histogram index
397 if ((mom < fMinMom) || (mom > fMaxMom)) return -1;
398 Int_t imom = ((Int_t) ((mom - fMinMom) / fWidMom));
399 return GetIndex(imom,ipid);
403 //_____________________________________________________________________________
404 Int_t AliTRDpidLQ::GetIndex(Int_t imom, Int_t ipid)
407 // Returns the histogram index
410 if ((ipid < 0) || (ipid >= kNpid)) return -1;
411 return imom * kNpid + ipid;