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