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