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