1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 /* History of cvs commits:
20 * Revision 1.13 2007/07/11 13:43:29 hristov
21 * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
23 * Revision 1.12 2007/06/11 20:43:06 hristov
24 * Changes required by the updated AliESDCaloCluster (Gustavo)
26 * Revision 1.11 2007/03/30 13:50:34 gustavo
27 * PID for particles with E < 5 GeV was not done, temporal solution found (Guenole)
29 * Revision 1.10 2007/03/09 14:34:11 gustavo
30 * Correct probability calculation, added missing initialization of data members
32 * Revision 1.9 2007/02/20 20:17:43 hristov
33 * Corrected array size, removed warnings (icc)
35 * Revision 1.8 2006/12/19 08:49:35 gustavo
36 * New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled when calling AliEMCALPID in AliEMCALReconstructor::FillESD()
40 // to compute PID for all the clusters in ESDs.root file
41 // the ESDs.root have to be in the same directory as the class
47 // AliEMCALPID::CalculPID(Energy,Lambda0)
48 // Calcul PID for all clusters in AliESDs.root file
49 // keep this function for the moment for a simple verification, could be removed
53 // AliEMCALPID::CalculPID(Energy,Lambda0)
54 // calcul PID Weght for a cluster with Energy, Lambda0 .
55 // Double_t PIDFinal[AliPID::kSPECIESN] is the standard PID for :
59 // kElectron : fPIDFinal[0]
60 // kMuon : fPIDFinal[1]
61 // kPion : fPIDFinal[2]
62 // kKaon : fPIDFinal[3]
63 // kProton : fPIDFinal[4]
64 // kPhoton : fPIDFinal[5]
65 // kPi0 : fPIDFinal[6]
66 // kNeutron : fPIDFinal[7]
67 // kKaon0 : fPIDFinal[8]
68 // kEleCon : fPIDFinal[9]
69 // kUnknown : fPIDFinal[10]
72 // PID[3] is a simple PID for
73 // Electron & Photon PID[0]
81 // --- ROOT system ---
83 // standard C++ includes
84 #include <Riostream.h>
93 #include "TClonesArray.h"
95 #include "TLorentzVector.h"
100 #include "TParticle.h"
103 // #include "AliRun.h"
104 // #include "AliRunLoader.h"
105 // #include "AliHeader.h"
106 // #include "AliLoader.h"
107 // #include "AliStack.h"
108 // #include "AliESDtrack.h"
109 // #include "AliESDEvent.h"
111 #include "AliEMCALPID.h"
112 #include "AliESDCaloCluster.h"
114 ClassImp(AliEMCALPID)
116 //______________________________________________
117 AliEMCALPID::AliEMCALPID():
118 fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE)
122 // Initialize all constant values which have to be used
123 // during PID algorithm execution
126 // set flag for printing to FALSE by default
129 // as a first step, all array elements are initialized to 0.0
131 for (i = 0; i < 6; i++) {
132 for (j = 0; j < 6; j++) {
133 fGamma[i][j] = fHadron[i][j] = fPiZero5to10[i][j] = fPiZero10to60[i][j] = 0.;
137 // then, only the ones which must be not zero are initialized
138 // while the others will remain to the value 0.0
140 fGamma[0][0] = 0.038022;
141 fGamma[0][1] = -0.0001883;
142 fGamma[0][2] = 5.449e-06;
144 fGamma[1][0] = 0.207313;
145 fGamma[1][1] = -0.000978;
146 fGamma[1][2] = 0.00001634;
148 fGamma[2][0] = 0.043364;
149 fGamma[2][1] = -0.0002048;
150 fGamma[2][2] = 8.661e-06;
151 fGamma[2][3] = -1.353e-07;
153 fGamma[3][0] = 0.265004;
154 fGamma[3][1] = 0.061298;
155 fGamma[3][2] = -0.003203;
156 fGamma[3][3] = 4.73e-05;
158 fGamma[4][0] = 0.243579;
159 fGamma[4][1] = -1.614e-05;
161 fGamma[5][0] = 0.002942;
162 fGamma[5][1] = -3.976e-05;
164 fHadron[0][0] = 0.011945 / 3.;
165 fHadron[0][1] = 0.000386 / 3.;
166 fHadron[0][2] = -0.000014 / 3.;
167 fHadron[0][3] = 1.336e-07 / 3.;
169 fHadron[1][0] = 0.496544;
170 fHadron[1][1] = -0.003226;
171 fHadron[1][2] = 0.00001678;
173 fHadron[2][0] = 0.144838;
174 fHadron[2][1] = -0.002954;
175 fHadron[2][2] = 0.00008754;
176 fHadron[2][3] = -7.587e-07;
178 fHadron[3][0] = 1.264461 / 7.;
179 fHadron[3][1] = 0.002097 / 7.;
181 fHadron[4][0] = 0.261950;
182 fHadron[4][1] = -0.001078;
183 fHadron[4][2] = 0.00003237;
184 fHadron[4][3] = -3.241e-07;
187 fHadron[5][0] = 0.010317;
194 fPiZero5to10[0][0] = 0.009138;
195 fPiZero5to10[0][1] = 0.0006377;
197 fPiZero5to10[1][0] = 0.08;
199 fPiZero5to10[2][0] = -0.061119;
200 fPiZero5to10[2][1] = 0.019013;
202 fPiZero5to10[3][0] = 0.2;
204 fPiZero5to10[4][0] = 0.252044;
205 fPiZero5to10[4][1] = -0.002315;
207 fPiZero5to10[5][0] = 0.002942;
208 fPiZero5to10[5][1] = -3.976e-05;
210 fPiZero10to60[0][0] = 0.009138;
211 fPiZero10to60[0][1] = 0.0006377;
213 fPiZero10to60[1][0] = 1.272837;
214 fPiZero10to60[1][1] = -0.069708;
215 fPiZero10to60[1][2] = 0.001568;
216 fPiZero10to60[1][3] = -1.162e-05;
218 fPiZero10to60[2][0] = 0.139703;
219 fPiZero10to60[2][1] = 0.003687;
220 fPiZero10to60[2][2] = -0.000568;
221 fPiZero10to60[2][3] = 1.498e-05;
222 fPiZero10to60[2][4] = -1.174e-07;
224 fPiZero10to60[3][0] = -0.826367;
225 fPiZero10to60[3][1] = 0.096951;
226 fPiZero10to60[3][2] = -0.002215;
227 fPiZero10to60[3][3] = 2.523e-05;
229 fPiZero10to60[4][0] = 0.249890;
230 fPiZero10to60[4][1] = -0.000063;
232 fPiZero10to60[5][0] = 0.002942;
233 fPiZero10to60[5][1] = -3.976e-05;
238 fReconstructor = kFALSE;
241 //______________________________________________
242 void AliEMCALPID::RunPID(AliESDEvent *esd)
245 // Make the PID for all the EMCAL clusters containedin the ESDs File
246 // but just gamma/PiO/Hadron
248 // trivial check against NULL object passed
251 AliInfo("NULL ESD object passed !!" );
255 Int_t nClusters = esd->GetNumberOfEMCALClusters();
256 Int_t firstCluster = esd->GetFirstEMCALCluster();
257 Double_t energy, lambda0;
258 for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
260 AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
262 lambda0 = clust->GetM02();
263 // verify cluster type
264 Int_t clusterType= clust->GetClusterType();
265 if (clusterType == AliESDCaloCluster::kClusterv1 && lambda0 != 0 && energy < 1000) {
268 // reject clusters with lambda0 = 0
271 ComputePID(energy, lambda0);
275 AliInfo("___________________________________________________");
276 AliInfo(Form( "Particle Energy = %f",energy));
277 AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
278 AliInfo("PIDWeight of the particle :" );
279 AliInfo(Form( " GAMMA : %f",fPID[0] ));
280 AliInfo(Form( " PiZero : %f",fPID[1] ));
281 AliInfo(Form( " HADRON : %f", fPID[2] ));
282 AliInfo("_________________________________________");
283 AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
284 AliInfo(Form( " kMuon : %f", fPIDFinal[1] ));
285 AliInfo(Form( " kPion : %f", fPIDFinal[2] ));
286 AliInfo(Form( " kKaon : %f", fPIDFinal[3] ));
287 AliInfo(Form( " kProton : %f", fPIDFinal[4] ));
288 AliInfo(Form( " kPhoton : %f", fPIDFinal[5] ));
289 AliInfo(Form( " kPi0 : %f", fPIDFinal[6] ));
290 AliInfo(Form( " kNeutron : %f", fPIDFinal[7] ));
291 AliInfo(Form( " kKaon0 : %f", fPIDFinal[8] ));
292 AliInfo(Form( " kEleCon : %f", fPIDFinal[9] ));
293 AliInfo(Form( " kUnknown : %f", fPIDFinal[10] ));
294 AliInfo("___________________________________________________");
297 if(fReconstructor) // In case it is called during reconstruction.
298 clust->SetPid(fPIDFinal);
299 } // end if (clusterType...)
300 } // end for (iCluster...)
303 //__________________________________________________________
304 void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
307 // This is the main command, which uses the distributions computed and parametrised,
308 // and gives the PID by the bayesian method.
311 if (energy<5){energy =6;}
314 TArrayD paramDistribGamma = DistLambda0(energy, 1);
315 TArrayD paramDistribPiZero = DistLambda0(energy, 2);
316 TArrayD paramDistribHadron = DistLambda0(energy, 3);
318 Bool_t norm = kFALSE;
320 fProbGamma = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0];
321 fProbGamma += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3];
322 fProbPiZero = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
323 fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
324 fProbHadron = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
325 fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
327 // compute PID Weight
328 fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
329 fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
330 fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
332 SetPID(fPIDWeight[0], 0);
333 SetPID(fPIDWeight[1], 1);
334 SetPID(fPIDWeight[2], 2);
336 // sortie ecran pid Weight only for control (= in english ???)
338 AliInfo(Form( "Energy in loop = %f", energy) );
339 AliInfo(Form( "Lambda0 in loop = %f", lambda0) );
340 AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) );
341 // AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) );
342 AliInfo(Form( "fProbaPiZero = %f", fProbPiZero ));
343 AliInfo(Form( "fProbaHadron = %f", fProbHadron) );
344 AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f", fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
345 AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] ));
346 AliInfo("********************************************************" );
349 fPIDFinal[0] = fPIDWeight[0]/2;
350 fPIDFinal[1] = fPIDWeight[2]/8;
351 fPIDFinal[2] = fPIDWeight[2]/8;
352 fPIDFinal[3] = fPIDWeight[2]/8;
353 fPIDFinal[4] = fPIDWeight[2]/8;
354 fPIDFinal[5] = fPIDWeight[0]/2;
355 fPIDFinal[6] = fPIDWeight[1] ;
356 fPIDFinal[7] = fPIDWeight[2]/8;
357 fPIDFinal[8] = fPIDWeight[2]/8;
358 fPIDFinal[9] = fPIDWeight[2]/8;
359 fPIDFinal[10] = fPIDWeight[2]/8;
362 //________________________________________________________
363 TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
366 // Compute the values of the parametrised distributions using the data initialised before.
368 Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
369 Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
370 TArrayD distributionParam(6);
374 constGauss = Polynomial(energy, fGamma[0]);
375 meanGauss = Polynomial(energy, fGamma[1]);
376 sigmaGauss = Polynomial(energy, fGamma[2]);
377 constLandau = Polynomial(energy, fGamma[3]);
378 mpvLandau = Polynomial(energy, fGamma[4]);
379 sigmaLandau = Polynomial(energy, fGamma[5]);
383 constGauss = Polynomial(energy, fPiZero5to10[0]);
384 meanGauss = Polynomial(energy, fPiZero5to10[1]);
385 sigmaGauss = Polynomial(energy, fPiZero5to10[2]);
386 constLandau = Polynomial(energy, fPiZero5to10[3]);
387 mpvLandau = Polynomial(energy, fPiZero5to10[4]);
388 sigmaLandau = Polynomial(energy, fPiZero5to10[5]);
391 constGauss = Polynomial(energy, fPiZero10to60[0]);
392 meanGauss = Polynomial(energy, fPiZero10to60[1]);
393 sigmaGauss = Polynomial(energy, fPiZero10to60[2]);
394 constLandau = Polynomial(energy, fPiZero10to60[3]);
395 mpvLandau = Polynomial(energy, fPiZero10to60[4]);
396 sigmaLandau = Polynomial(energy, fPiZero10to60[5]);
400 constGauss = Polynomial(energy, fHadron[0]);
401 meanGauss = Polynomial(energy, fHadron[1]);
402 sigmaGauss = Polynomial(energy, fHadron[2]);
403 constLandau = Polynomial(energy, fHadron[3]);
404 mpvLandau = Polynomial(energy, fHadron[4]);
405 sigmaLandau = Polynomial(energy, fHadron[5]);
409 distributionParam[0] = constGauss;
410 distributionParam[1] = meanGauss;
411 distributionParam[2] = sigmaGauss;
412 distributionParam[3] = constLandau;
413 distributionParam[4] = mpvLandau;
414 distributionParam[5] = sigmaLandau;
416 return distributionParam;
419 //_______________________________________________________
420 Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
423 // Compute a polynomial for a given value of 'x'
424 // with the array of parameters passed as the second arg
429 y += params[2] * x * x;
430 y += params[3] * x * x * x;
431 y += params[4] * x * x * x * x;
432 y += params[5] * x * x * x * x * x;