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.14 2007/07/26 16:54:53 morsch
21 * Changes in AliESDEvent fwd declarartions.
23 * Revision 1.13 2007/07/11 13:43:29 hristov
24 * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
26 * Revision 1.12 2007/06/11 20:43:06 hristov
27 * Changes required by the updated AliESDCaloCluster (Gustavo)
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)
32 * Revision 1.10 2007/03/09 14:34:11 gustavo
33 * Correct probability calculation, added missing initialization of data members
35 * Revision 1.9 2007/02/20 20:17:43 hristov
36 * Corrected array size, removed warnings (icc)
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()
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
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
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 :
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]
75 // PID[3] is a simple PID for
76 // Electron & Photon PID[0]
84 // --- ROOT system ---
86 // standard C++ includes
87 #include <Riostream.h>
96 #include "TClonesArray.h"
98 #include "TLorentzVector.h"
103 #include "TParticle.h"
106 // #include "AliRun.h"
107 // #include "AliRunLoader.h"
108 // #include "AliHeader.h"
109 // #include "AliLoader.h"
110 // #include "AliStack.h"
111 // #include "AliESDtrack.h"
112 // #include "AliESDEvent.h"
114 #include "AliEMCALPID.h"
115 #include "AliESDCaloCluster.h"
117 ClassImp(AliEMCALPID)
119 //______________________________________________
120 AliEMCALPID::AliEMCALPID():
121 fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE)
125 // Initialize all constant values which have to be used
126 // during PID algorithm execution
129 // set flag for printing to FALSE by default
132 // as a first step, all array elements are initialized to 0.0
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.;
140 // then, only the ones which must be not zero are initialized
141 // while the others will remain to the value 0.0
143 fGamma[0][0] = 0.038022;
144 fGamma[0][1] = -0.0001883;
145 fGamma[0][2] = 5.449e-06;
147 fGamma[1][0] = 0.207313;
148 fGamma[1][1] = -0.000978;
149 fGamma[1][2] = 0.00001634;
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;
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;
161 fGamma[4][0] = 0.243579;
162 fGamma[4][1] = -1.614e-05;
164 fGamma[5][0] = 0.002942;
165 fGamma[5][1] = -3.976e-05;
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.;
172 fHadron[1][0] = 0.496544;
173 fHadron[1][1] = -0.003226;
174 fHadron[1][2] = 0.00001678;
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;
181 fHadron[3][0] = 1.264461 / 7.;
182 fHadron[3][1] = 0.002097 / 7.;
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;
190 fHadron[5][0] = 0.010317;
197 fPiZero5to10[0][0] = 0.009138;
198 fPiZero5to10[0][1] = 0.0006377;
200 fPiZero5to10[1][0] = 0.08;
202 fPiZero5to10[2][0] = -0.061119;
203 fPiZero5to10[2][1] = 0.019013;
205 fPiZero5to10[3][0] = 0.2;
207 fPiZero5to10[4][0] = 0.252044;
208 fPiZero5to10[4][1] = -0.002315;
210 fPiZero5to10[5][0] = 0.002942;
211 fPiZero5to10[5][1] = -3.976e-05;
213 fPiZero10to60[0][0] = 0.009138;
214 fPiZero10to60[0][1] = 0.0006377;
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;
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;
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;
232 fPiZero10to60[4][0] = 0.249890;
233 fPiZero10to60[4][1] = -0.000063;
235 fPiZero10to60[5][0] = 0.002942;
236 fPiZero10to60[5][1] = -3.976e-05;
241 fReconstructor = kFALSE;
244 //______________________________________________
245 void AliEMCALPID::RunPID(AliESDEvent *esd)
248 // Make the PID for all the EMCAL clusters containedin the ESDs File
249 // but just gamma/PiO/Hadron
251 // trivial check against NULL object passed
254 AliInfo("NULL ESD object passed !!" );
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++) {
263 AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
265 lambda0 = clust->GetM02();
266 // verify cluster type
267 Int_t clusterType= clust->GetClusterType();
268 if (clusterType == AliESDCaloCluster::kEMCALClusterv1 && lambda0 != 0 && energy < 1000) {
271 // reject clusters with lambda0 = 0
274 ComputePID(energy, lambda0);
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("___________________________________________________");
300 if(fReconstructor) // In case it is called during reconstruction.
301 clust->SetPid(fPIDFinal);
302 } // end if (clusterType...)
303 } // end for (iCluster...)
306 //__________________________________________________________
307 void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
310 // This is the main command, which uses the distributions computed and parametrised,
311 // and gives the PID by the bayesian method.
314 if (energy<5){energy =6;}
317 TArrayD paramDistribGamma = DistLambda0(energy, 1);
318 TArrayD paramDistribPiZero = DistLambda0(energy, 2);
319 TArrayD paramDistribHadron = DistLambda0(energy, 3);
321 Bool_t norm = kFALSE;
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];
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);
335 SetPID(fPIDWeight[0], 0);
336 SetPID(fPIDWeight[1], 1);
337 SetPID(fPIDWeight[2], 2);
339 // sortie ecran pid Weight only for control (= in english ???)
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("********************************************************" );
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;
365 //________________________________________________________
366 TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
369 // Compute the values of the parametrised distributions using the data initialised before.
371 Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
372 Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
373 TArrayD distributionParam(6);
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]);
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]);
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]);
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]);
412 distributionParam[0] = constGauss;
413 distributionParam[1] = meanGauss;
414 distributionParam[2] = sigmaGauss;
415 distributionParam[3] = constLandau;
416 distributionParam[4] = mpvLandau;
417 distributionParam[5] = sigmaLandau;
419 return distributionParam;
422 //_______________________________________________________
423 Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
426 // Compute a polynomial for a given value of 'x'
427 // with the array of parameters passed as the second arg
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;