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