]> 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
New classes added
[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
16/*
17$Log$
a3d851c2 18Revision 1.1 2001/11/06 17:19:41 cblume
19Add detailed geometry and simple simulator
20
16bf9884 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
52ClassImp(AliTRDpidLQ)
53
54//_____________________________________________________________________________
55AliTRDpidLQ::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//_____________________________________________________________________________
78AliTRDpidLQ::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//_____________________________________________________________________________
101AliTRDpidLQ::AliTRDpidLQ(const AliTRDpidLQ &p)
102{
103 //
104 // AliTRDpidLQ copy constructor
105 //
106
107 ((AliTRDpidLQ &) p).Copy(*this);
108
109}
110
111//_____________________________________________________________________________
112AliTRDpidLQ::~AliTRDpidLQ()
113{
114 //
115 // AliTRDpidLQ destructor
116 //
117
118 if (fHist) {
119 fHist->Delete();
120 delete fHist;
121 }
122
123}
124
125//_____________________________________________________________________________
126AliTRDpidLQ &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//_____________________________________________________________________________
138void 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//_____________________________________________________________________________
161Bool_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//_____________________________________________________________________________
179Bool_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();
a3d851c2 186 Float_t * charge = new Float_t[kNpla];
187 Int_t * nCluster = new Int_t[kNpla];
16bf9884 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 {
a3d851c2 241 delete [] charge;
242 delete [] nCluster;
16bf9884 243 return kTRUE;
244 }
245
246 pSum = pEl + pPi;
a3d851c2 247 if (pSum <= 0) {
248 delete [] charge;
249 delete [] nCluster;
250 return kFALSE;
251 }
16bf9884 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
a3d851c2 260 delete [] charge;
261 delete [] nCluster;
16bf9884 262 return kTRUE;
263
264}
265
266//_____________________________________________________________________________
267Bool_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//_____________________________________________________________________________
344Bool_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
a3d851c2 354 Float_t * charge = new Float_t[kNpla];
355 Int_t * nCluster = new Int_t[kNpla];
16bf9884 356 Float_t mom = t->GetP();
357 Int_t ipid = MCpid(t);
358 TH1F *hTmp = NULL;
359
a3d851c2 360 if (!SumCharge(t,charge,nCluster)) {
361 delete [] charge;
362 delete [] nCluster;
363 return kFALSE;
364 }
16bf9884 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
a3d851c2 377 delete [] charge;
378 delete [] nCluster;
16bf9884 379 return kTRUE;
380
381}
382
383//_____________________________________________________________________________
384Int_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//_____________________________________________________________________________
399Int_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//_____________________________________________________________________________
412Int_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}