Memory leak fixed
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALPID.cxx
CommitLineData
dc293ae9 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/* $Id$ */
17/* History of cvs commits:
18 *
19 * $Log$
7592dfc4 20 * Revision 1.11 2007/03/30 13:50:34 gustavo
21 * PID for particles with E < 5 GeV was not done, temporal solution found (Guenole)
22 *
4ce73a00 23 * Revision 1.10 2007/03/09 14:34:11 gustavo
24 * Correct probability calculation, added missing initialization of data members
25 *
4cfd8ef9 26 * Revision 1.9 2007/02/20 20:17:43 hristov
27 * Corrected array size, removed warnings (icc)
28 *
d69ab345 29 * Revision 1.8 2006/12/19 08:49:35 gustavo
30 * New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled when calling AliEMCALPID in AliEMCALReconstructor::FillESD()
31 *
dc293ae9 32 *
33 */
34// to compute PID for all the clusters in ESDs.root file
35// the ESDs.root have to be in the same directory as the class
36//
37//
38//
39//
40//
41// AliEMCALPID::CalculPID(Energy,Lambda0)
42// Calcul PID for all clusters in AliESDs.root file
43// keep this function for the moment for a simple verification, could be removed
44//
45//
46//
47// AliEMCALPID::CalculPID(Energy,Lambda0)
48// calcul PID Weght for a cluster with Energy, Lambda0 .
49// Double_t PIDFinal[AliPID::kSPECIESN] is the standard PID for :
50//
51//
52//
53// kElectron : fPIDFinal[0]
54// kMuon : fPIDFinal[1]
55// kPion : fPIDFinal[2]
56// kKaon : fPIDFinal[3]
57// kProton : fPIDFinal[4]
58// kPhoton : fPIDFinal[5]
59// kPi0 : fPIDFinal[6]
60// kNeutron : fPIDFinal[7]
61// kKaon0 : fPIDFinal[8]
62// kEleCon : fPIDFinal[9]
63// kUnknown : fPIDFinal[10]
64//
65//
66// PID[3] is a simple PID for
67// Electron & Photon PID[0]
68// Pi0 PID[1]
69// Hadron PID[2]
70//
71//
72//
73//
74//
75// --- ROOT system ---
76
77// standard C++ includes
78#include <Riostream.h>
79// #include <cstdlib>
80// #include <cmath>
81
82// ROOT includes
83#include "TTree.h"
84#include "TStyle.h"
85#include "TVector3.h"
86#include "TBranch.h"
87#include "TClonesArray.h"
88#include "TCanvas.h"
89#include "TLorentzVector.h"
90#include "TMath.h"
91#include "TFile.h"
92#include "TH1.h"
93#include "TH2.h"
94#include "TParticle.h"
95
96// // STEER includes
97// #include "AliRun.h"
98// #include "AliRunLoader.h"
99// #include "AliHeader.h"
100// #include "AliLoader.h"
101// #include "AliStack.h"
102// #include "AliESDtrack.h"
103// #include "AliESD.h"
104#include "AliLog.h"
105#include "AliEMCALPID.h"
106
107ClassImp(AliEMCALPID)
108
4cfd8ef9 109//______________________________________________
110 AliEMCALPID::AliEMCALPID():
111 fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE)
dc293ae9 112{
4cfd8ef9 113 //
114 // Constructor.
115 // Initialize all constant values which have to be used
116 // during PID algorithm execution
117 //
118
119 // set flag for printing to FALSE by default
120 fPrintInfo = kFALSE;
121
122 // as a first step, all array elements are initialized to 0.0
123 Int_t i, j;
124 for (i = 0; i < 6; i++) {
125 for (j = 0; j < 6; j++) {
126 fGamma[i][j] = fHadron[i][j] = fPiZero5to10[i][j] = fPiZero10to60[i][j] = 0.;
127 }
128 }
129
130 // then, only the ones which must be not zero are initialized
131 // while the others will remain to the value 0.0
132
133 fGamma[0][0] = 0.038022;
134 fGamma[0][1] = -0.0001883;
135 fGamma[0][2] = 5.449e-06;
136
137 fGamma[1][0] = 0.207313;
138 fGamma[1][1] = -0.000978;
139 fGamma[1][2] = 0.00001634;
140
141 fGamma[2][0] = 0.043364;
142 fGamma[2][1] = -0.0002048;
143 fGamma[2][2] = 8.661e-06;
144 fGamma[2][3] = -1.353e-07;
145
146 fGamma[3][0] = 0.265004;
147 fGamma[3][1] = 0.061298;
148 fGamma[3][2] = -0.003203;
149 fGamma[3][3] = 4.73e-05;
150
151 fGamma[4][0] = 0.243579;
152 fGamma[4][1] = -1.614e-05;
153
154 fGamma[5][0] = 0.002942;
155 fGamma[5][1] = -3.976e-05;
156
157 fHadron[0][0] = 0.011945 / 3.;
158 fHadron[0][1] = 0.000386 / 3.;
159 fHadron[0][2] = -0.000014 / 3.;
160 fHadron[0][3] = 1.336e-07 / 3.;
161
162 fHadron[1][0] = 0.496544;
163 fHadron[1][1] = -0.003226;
164 fHadron[1][2] = 0.00001678;
165
166 fHadron[2][0] = 0.144838;
167 fHadron[2][1] = -0.002954;
168 fHadron[2][2] = 0.00008754;
169 fHadron[2][3] = -7.587e-07;
170
171 fHadron[3][0] = 1.264461 / 7.;
172 fHadron[3][1] = 0.002097 / 7.;
173
174 fHadron[4][0] = 0.261950;
175 fHadron[4][1] = -0.001078;
176 fHadron[4][2] = 0.00003237;
177 fHadron[4][3] = -3.241e-07;
178 fHadron[4][4] = 0.;
179 fHadron[4][5] = 0.;
180 fHadron[5][0] = 0.010317;
181 fHadron[5][1] = 0.;
182 fHadron[5][2] = 0.;
183 fHadron[5][3] = 0.;
184 fHadron[5][4] = 0.;
185 fHadron[5][5] = 0.;
186
187 fPiZero5to10[0][0] = 0.009138;
188 fPiZero5to10[0][1] = 0.0006377;
189
190 fPiZero5to10[1][0] = 0.08;
191
192 fPiZero5to10[2][0] = -0.061119;
193 fPiZero5to10[2][1] = 0.019013;
194
195 fPiZero5to10[3][0] = 0.2;
196
197 fPiZero5to10[4][0] = 0.252044;
198 fPiZero5to10[4][1] = -0.002315;
199
200 fPiZero5to10[5][0] = 0.002942;
201 fPiZero5to10[5][1] = -3.976e-05;
202
203 fPiZero10to60[0][0] = 0.009138;
204 fPiZero10to60[0][1] = 0.0006377;
205
206 fPiZero10to60[1][0] = 1.272837;
207 fPiZero10to60[1][1] = -0.069708;
208 fPiZero10to60[1][2] = 0.001568;
209 fPiZero10to60[1][3] = -1.162e-05;
210
211 fPiZero10to60[2][0] = 0.139703;
212 fPiZero10to60[2][1] = 0.003687;
213 fPiZero10to60[2][2] = -0.000568;
214 fPiZero10to60[2][3] = 1.498e-05;
215 fPiZero10to60[2][4] = -1.174e-07;
216
217 fPiZero10to60[3][0] = -0.826367;
218 fPiZero10to60[3][1] = 0.096951;
219 fPiZero10to60[3][2] = -0.002215;
220 fPiZero10to60[3][3] = 2.523e-05;
221
222 fPiZero10to60[4][0] = 0.249890;
223 fPiZero10to60[4][1] = -0.000063;
224
225 fPiZero10to60[5][0] = 0.002942;
226 fPiZero10to60[5][1] = -3.976e-05;
227
228 fPIDWeight[0] = -1;
229 fPIDWeight[1] = -1;
230 fPIDWeight[2] = -1;
231 fReconstructor = kFALSE;
dc293ae9 232}
4cfd8ef9 233
234//______________________________________________
dc293ae9 235void AliEMCALPID::RunPID(AliESD *esd)
236{
237//
238// Make the PID for all the EMCAL clusters containedin the ESDs File
239// but just gamma/PiO/Hadron
240//
241 // trivial check against NULL object passed
4cfd8ef9 242
dc293ae9 243 if (esd == 0x0) {
4ce73a00 244 AliInfo("NULL ESD object passed !!" );
dc293ae9 245 return ;
246 }
4cfd8ef9 247
dc293ae9 248 Int_t nClusters = esd->GetNumberOfEMCALClusters();
249 Int_t firstCluster = esd->GetFirstEMCALCluster();
250 Double_t energy, lambda0;
251 for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
252
253 AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
7592dfc4 254 energy = clust->E();
dc293ae9 255 lambda0 = clust->GetM02();
256 // verify cluster type
257 Int_t clusterType= clust->GetClusterType();
4ce73a00 258 if (clusterType == AliESDCaloCluster::kClusterv1 && lambda0 != 0 && energy < 1000) {
259
260
dc293ae9 261 // reject clusters with lambda0 = 0
4ce73a00 262
263
dc293ae9 264 ComputePID(energy, lambda0);
4ce73a00 265
266
dc293ae9 267 if (fPrintInfo) {
268 AliInfo("___________________________________________________");
269 AliInfo(Form( "Particle Energy = %f",energy));
270 AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
271 AliInfo("PIDWeight of the particle :" );
272 AliInfo(Form( " GAMMA : %f",fPID[0] ));
273 AliInfo(Form( " PiZero : %f",fPID[1] ));
274 AliInfo(Form( " HADRON : %f", fPID[2] ));
275 AliInfo("_________________________________________");
276 AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
277 AliInfo(Form( " kMuon : %f", fPIDFinal[1] ));
278 AliInfo(Form( " kPion : %f", fPIDFinal[2] ));
279 AliInfo(Form( " kKaon : %f", fPIDFinal[3] ));
280 AliInfo(Form( " kProton : %f", fPIDFinal[4] ));
281 AliInfo(Form( " kPhoton : %f", fPIDFinal[5] ));
282 AliInfo(Form( " kPi0 : %f", fPIDFinal[6] ));
283 AliInfo(Form( " kNeutron : %f", fPIDFinal[7] ));
284 AliInfo(Form( " kKaon0 : %f", fPIDFinal[8] ));
285 AliInfo(Form( " kEleCon : %f", fPIDFinal[9] ));
286 AliInfo(Form( " kUnknown : %f", fPIDFinal[10] ));
287 AliInfo("___________________________________________________");
288 }
4cfd8ef9 289
dc293ae9 290 if(fReconstructor) // In case it is called during reconstruction.
291 clust->SetPid(fPIDFinal);
292 } // end if (clusterType...)
293 } // end for (iCluster...)
294}
4cfd8ef9 295
296//__________________________________________________________
dc293ae9 297void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
298{
299//
300// This is the main command, which uses the distributions computed and parametrised,
301// and gives the PID by the bayesian method.
302//
4ce73a00 303
304if (energy<5){energy =6;}
305
306
4cfd8ef9 307 TArrayD paramDistribGamma = DistLambda0(energy, 1);
308 TArrayD paramDistribPiZero = DistLambda0(energy, 2);
309 TArrayD paramDistribHadron = DistLambda0(energy, 3);
310
311 Bool_t norm = kFALSE;
312
313 fProbGamma = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0];
314 fProbGamma += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3];
315 fProbPiZero = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
316 fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
317 fProbHadron = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
318 fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
319
320 // compute PID Weight
321 fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
322 fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
323 fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
324
325 SetPID(fPIDWeight[0], 0);
326 SetPID(fPIDWeight[1], 1);
327 SetPID(fPIDWeight[2], 2);
328
329 // sortie ecran pid Weight only for control (= in english ???)
330 if (fPrintInfo) {
331 AliInfo(Form( "Energy in loop = %f", energy) );
332 AliInfo(Form( "Lambda0 in loop = %f", lambda0) );
333 AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) );
334 // AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) );
335 AliInfo(Form( "fProbaPiZero = %f", fProbPiZero ));
336 AliInfo(Form( "fProbaHadron = %f", fProbHadron) );
337 AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f", fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
338 AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] ));
339 AliInfo("********************************************************" );
340 }
341
342 fPIDFinal[0] = fPIDWeight[0]/2;
343 fPIDFinal[1] = fPIDWeight[2]/8;
344 fPIDFinal[2] = fPIDWeight[2]/8;
345 fPIDFinal[3] = fPIDWeight[2]/8;
346 fPIDFinal[4] = fPIDWeight[2]/8;
347 fPIDFinal[5] = fPIDWeight[0]/2;
348 fPIDFinal[6] = fPIDWeight[1] ;
349 fPIDFinal[7] = fPIDWeight[2]/8;
350 fPIDFinal[8] = fPIDWeight[2]/8;
351 fPIDFinal[9] = fPIDWeight[2]/8;
352 fPIDFinal[10] = fPIDWeight[2]/8;
dc293ae9 353}
4cfd8ef9 354
355//________________________________________________________
dc293ae9 356TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
357{
4cfd8ef9 358 //
359 // Compute the values of the parametrised distributions using the data initialised before.
360 //
361 Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
362 Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
363 TArrayD distributionParam(6);
364
365 switch (type) {
366 case 1:
367 constGauss = Polynomial(energy, fGamma[0]);
368 meanGauss = Polynomial(energy, fGamma[1]);
369 sigmaGauss = Polynomial(energy, fGamma[2]);
370 constLandau = Polynomial(energy, fGamma[3]);
371 mpvLandau = Polynomial(energy, fGamma[4]);
372 sigmaLandau = Polynomial(energy, fGamma[5]);
373 break;
374 case 2:
375 if (energy < 10) {
376 constGauss = Polynomial(energy, fPiZero5to10[0]);
377 meanGauss = Polynomial(energy, fPiZero5to10[1]);
378 sigmaGauss = Polynomial(energy, fPiZero5to10[2]);
379 constLandau = Polynomial(energy, fPiZero5to10[3]);
380 mpvLandau = Polynomial(energy, fPiZero5to10[4]);
381 sigmaLandau = Polynomial(energy, fPiZero5to10[5]);
382 }
383 else {
384 constGauss = Polynomial(energy, fPiZero10to60[0]);
385 meanGauss = Polynomial(energy, fPiZero10to60[1]);
386 sigmaGauss = Polynomial(energy, fPiZero10to60[2]);
387 constLandau = Polynomial(energy, fPiZero10to60[3]);
388 mpvLandau = Polynomial(energy, fPiZero10to60[4]);
389 sigmaLandau = Polynomial(energy, fPiZero10to60[5]);
390 }
391 break;
392 case 3:
393 constGauss = Polynomial(energy, fHadron[0]);
394 meanGauss = Polynomial(energy, fHadron[1]);
395 sigmaGauss = Polynomial(energy, fHadron[2]);
396 constLandau = Polynomial(energy, fHadron[3]);
397 mpvLandau = Polynomial(energy, fHadron[4]);
398 sigmaLandau = Polynomial(energy, fHadron[5]);
399 break;
400 }
401
402 distributionParam[0] = constGauss;
403 distributionParam[1] = meanGauss;
404 distributionParam[2] = sigmaGauss;
405 distributionParam[3] = constLandau;
406 distributionParam[4] = mpvLandau;
407 distributionParam[5] = sigmaLandau;
408
409 return distributionParam;
dc293ae9 410}
4cfd8ef9 411
412//_______________________________________________________
dc293ae9 413Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
414{
4cfd8ef9 415 //
416 // Compute a polynomial for a given value of 'x'
417 // with the array of parameters passed as the second arg
418 //
419 Double_t y;
420 y = params[0];
421 y += params[1] * x;
422 y += params[2] * x * x;
423 y += params[3] * x * x * x;
424 y += params[4] * x * x * x * x;
425 y += params[5] * x * x * x * x * x;
426
427 return y;
dc293ae9 428}