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