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