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