PID parameters tuned for low and high flux environments (Marie Germain)
[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
18 //   Compute PID weights for all the clusters that are in AliESDs.root file
19 //   the AliESDs.root have to be in the same directory as the class
20 //
21 //   and do:    
22 //   AliEMCALPID *pid = new AliEMCALPID(kFALSE); // this calls the constructor which avoids the call to recparam 
23 //   pid->SetReconstructor(kFALSE);
24 //   pid->SetPrintInfo(kTRUE);
25 //   pid->SetHighFluxParam(); //   pid->SetLowFluxParam(); 
26 //   
27 //   then in cluster loop do
28 //   pid->ComputePID(energy, lambda0);
29 //        
30 //        Compute PID Weight for all clusters in AliESDs.root file
31 //        keep this function for the moment for a simple verification, could be removed
32 //
33 //   pid->GetPIDFinal(idx) gives the probabilities
34 //
35 //   Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
36 //
37 //
38 //
39 //      kElectron :  fPIDFinal[0]
40 //      kMuon     :  fPIDFinal[1]
41 //      kPion     :  fPIDFinal[2]
42 //      kKaon     :  fPIDFinal[3]
43 //      kProton   :  fPIDFinal[4]
44 //      kPhoton   :  fPIDFinal[5]
45 //      kPi0      :  fPIDFinal[6]
46 //      kNeutron  :  fPIDFinal[7]
47 //      kKaon0    :  fPIDFinal[8]
48 //      kEleCon   :  fPIDFinal[9]
49 //      kUnknown  :  fPIDFinal[10]
50 //
51 //
52 //    PID[3] is a simple PID for
53 //      Electron & Photon  PID[0]
54 //                    Pi0  PID[1]
55 //                 Hadron  PID[2]
56 //
57 // --- standard c ---
58
59 // standard C++ includes
60 //#include <Riostream.h>
61
62 // ROOT includes
63 //#include "TTree.h"
64 //#include "TVector3.h"
65 //#include "TBranch.h"
66 //#include "TClonesArray.h"
67 //#include "TLorentzVector.h"
68 #include "TMath.h"
69 //#include "TRefArray.h"
70 #include "TArrayD.h"
71
72 // STEER includes
73 #include "AliESDEvent.h"
74 //#include "AliLog.h"
75 #include "AliEMCALPID.h"
76 #include "AliESDCaloCluster.h"
77 //#include "AliEMCALRecParam.h"
78 #include "AliEMCALReconstructor.h"
79
80   
81 ClassImp(AliEMCALPID)
82   
83 //______________________________________________
84   AliEMCALPID::AliEMCALPID():
85     fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.), fWeightHadronEnergy(1.), fWeightGammaEnergy(1.),fWeightPiZeroEnergy(1.),fReconstructor(kTRUE)
86 {
87   //
88   // Constructor.
89   // Initialize all constant values which have to be used
90   // during PID algorithm execution
91   //
92   
93   InitParameters(); 
94   
95   
96 }
97
98 //______________________________________________
99 AliEMCALPID::AliEMCALPID(Bool_t reconstructor):
100   fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.), fWeightHadronEnergy(1.), fWeightGammaEnergy(1.),fWeightPiZeroEnergy(1.),fReconstructor(reconstructor)
101 {
102   //
103   // Constructor.
104   // Initialize all constant values which have to be used
105   // during PID algorithm execution called when used in standalone mode 
106   //
107   
108   InitParameters(); 
109   
110 }
111
112 //______________________________________________
113 void AliEMCALPID::RunPID(AliESDEvent *esd)
114 {
115   //
116   // Make the PID for all the EMCAL clusters containedin the ESDs File
117   // but just gamma/PiO/Hadron
118   //
119   // trivial check against NULL object passed
120   
121   if (esd == 0x0) {
122     AliInfo("NULL ESD object passed !!" );
123     return ;
124   }
125   
126   Int_t nClusters = esd->GetNumberOfCaloClusters();
127   Int_t firstCluster = 0;
128   Double_t energy, lambda0;
129   for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
130     
131     AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
132     if (!clust->IsEMCAL()) continue ; 
133     
134     energy = clust->E();
135     lambda0 = clust->GetM02();
136     // verify cluster type
137     Int_t clusterType= clust->GetClusterType();
138     if (clusterType == AliESDCaloCluster::kEMCALClusterv1 && lambda0 != 0  && energy < 1000) {
139       
140       //      if (lambda0 != 0  && energy < 1000) {
141       
142       // reject clusters with lambda0 = 0
143       
144       
145       ComputePID(energy, lambda0);
146       
147       
148       if (fPrintInfo) {
149         AliInfo("___________________________________________________");
150         AliInfo(Form( "Particle Energy = %f",energy));
151         AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
152         AliInfo("PIDWeight of the particle :" );
153         AliInfo(Form( " GAMMA  : %f",fPID[0] ));
154         AliInfo(Form( " PiZero : %f",fPID[1] ));
155         AliInfo(Form( " HADRON : %f", fPID[2] ));
156         AliInfo("_________________________________________");
157         AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
158         AliInfo(Form( " kMuon     : %f", fPIDFinal[1] ));
159         AliInfo(Form( " kPion       : %f", fPIDFinal[2] ));
160         AliInfo(Form( " kKaon       : %f", fPIDFinal[3] ));
161         AliInfo(Form( " kProton   : %f", fPIDFinal[4] ));
162         AliInfo(Form( " kPhoton   : %f", fPIDFinal[5] ));
163         AliInfo(Form( " kPi0        : %f", fPIDFinal[6] ));
164         AliInfo(Form( " kNeutron  : %f", fPIDFinal[7] ));
165         AliInfo(Form( " kKaon0  : %f", fPIDFinal[8] ));
166         AliInfo(Form( " kEleCon   : %f", fPIDFinal[9] ));
167         AliInfo(Form( " kUnknown  : %f", fPIDFinal[10] ));
168         AliInfo("___________________________________________________");
169       }
170       
171       if(fReconstructor){ // In case it is called during reconstruction.
172         //      cout << "#############On remplit l esd avec les PIDWeight##########" << endl;
173         clust->SetPid(fPIDFinal);}
174     } // end if (clusterType...)
175   } // end for (iCluster...)
176 }
177
178 //__________________________________________________________
179 void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
180 {
181 //
182 // This is the main command, which uses the distributions computed and parametrised, 
183 // and gives the PID by the bayesian method.
184 //
185 //   cout << "ENERGY  " <<energy <<" lambda0 "<< lambda0<<  endl;
186   
187   Double_t weightGammaEnergy  = DistEnergy(energy, 1);
188   Double_t weightPiZeroEnergy = DistEnergy(energy, 2);
189   Double_t weightHadronEnergy = DistEnergy(energy, 3);
190   
191   //Double_t weightHadronEnergy = 1.;
192   
193   Double_t energyhadron=energy;
194   if(energyhadron<1.)energyhadron=1.; // no energy dependance of  parametrisation for hadrons below 1 GeV
195   if (energy<2){energy =2;} // no energy dependance of parametrisation for gamma and pi0 below 2 GeV
196   
197   if (energy>55){
198     energy =55.;
199     energyhadron=55.;
200   } // same parametrisation for gamma and hadrons above 55 GeV 
201   //   for the pi0 above 55GeV the 2 gammas supperposed no way to distinguish from real gamma  PIDWeight[1]=0
202   
203   TArrayD paramDistribGamma  = DistLambda0(energy, 1);
204   TArrayD paramDistribPiZero = DistLambda0(energy, 2);
205   TArrayD paramDistribHadron = DistLambda0(energyhadron, 3);
206   
207   Bool_t norm = kFALSE;
208   
209   
210   fProbGamma   = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0];
211   fProbGamma  += TMath::Landau(((1-paramDistribGamma[4])-lambda0),paramDistribGamma[4],paramDistribGamma[5],norm)* paramDistribGamma[3];
212   if(fProbGamma<0.)fProbGamma=0.;
213   
214   fProbGamma = fProbGamma*weightGammaEnergy;
215   
216   if(energy>10. || energy < 55.){
217     fProbPiZero  = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
218     fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
219     if(fProbPiZero<0. || energy<5.)fProbPiZero=0.;
220     fProbPiZero = fProbPiZero*weightPiZeroEnergy;
221   }
222   else {
223     fProbPiZero = 0.;
224   }
225   
226   fProbHadron  = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
227   fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
228   if(fProbHadron<0.)fProbHadron=0.;
229   fProbHadron = fProbHadron*weightHadronEnergy; // to take into account the probability for a hadron to have a given reconstructed energy 
230   
231   // compute PID Weight
232   if( (fProbGamma + fProbPiZero + fProbHadron)>0.){
233     fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
234     fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
235     fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
236   }
237   else{   
238 // cases where  energy and lambda0 large,  probably du to 2 clusters folded the clusters PID not assigned to hadron nor Pi0 nor gammas
239     fPIDWeight[0] = 0.;
240     fPIDWeight[1] = 0.;
241     fPIDWeight[2] = 0.;
242   }
243   
244   
245   // cout << " PID[0] "<<  fPIDWeight[0] <<  " PID[1] "<<  fPIDWeight[1] <<  " PID[2] "<<  fPIDWeight[2] << endl;
246   
247   SetPID(fPIDWeight[0], 0);
248   SetPID(fPIDWeight[1], 1);
249   SetPID(fPIDWeight[2], 2);
250   
251   // print  pid Weight only for control (= in english ???)
252   if (fPrintInfo) {
253     AliInfo(Form( "Energy in loop = %f", energy) );
254     AliInfo(Form( "Lambda0 in loop = %f", lambda0) );
255     AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) );
256     AliInfo(Form( "fProbaPiZero = %f", fProbPiZero ));
257     AliInfo(Form( "fProbaHadron = %f", fProbHadron) );
258     AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f",  fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
259     AliInfo("********************************************************" );
260   }
261   
262   fPIDFinal[0]  = fPIDWeight[0]/2; // photon
263   fPIDFinal[1]  = fPIDWeight[2]/8;
264   fPIDFinal[2]  = fPIDWeight[2]/8;
265   fPIDFinal[3]  = fPIDWeight[2]/8;
266   fPIDFinal[4]  = fPIDWeight[2]/8;
267   fPIDFinal[5]  = fPIDWeight[0]/2; // electron
268   fPIDFinal[6]  = fPIDWeight[1]  ; // Pi0
269   fPIDFinal[7]  = fPIDWeight[2]/8;
270   fPIDFinal[8]  = fPIDWeight[2]/8;
271   fPIDFinal[9]  = fPIDWeight[2]/8;
272   fPIDFinal[10] = fPIDWeight[2]/8;
273
274 }
275
276
277
278
279 //________________________________________________________
280 TArrayD AliEMCALPID::DistLambda0(const Double_t energy, const Int_t type) 
281 {
282   //
283   // Compute the values of the parametrised distributions using the data initialised before.
284   //
285   Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
286   Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
287   TArrayD  distributionParam(6);
288   
289   switch (type) {
290     
291   case 1:
292     
293     constGauss  = PolynomialMixed2(energy, fGamma[0]);
294     meanGauss   = PolynomialMixed2(energy, fGamma[1]);
295     sigmaGauss  = PolynomialMixed2(energy, fGamma[2]);
296     constLandau = PolynomialMixed2(energy, fGamma[3]);
297     mpvLandau   = PolynomialMixed2(energy, fGamma[4]);
298     sigmaLandau = PolynomialMixed2(energy, fGamma[5]);
299    break;
300
301   case 2:
302
303     constGauss  = PolynomialMixed2(energy, fPiZero[0]);
304     meanGauss   = PolynomialMixed2(energy, fPiZero[1]);
305     sigmaGauss  = PolynomialMixed2(energy, fPiZero[2]);
306     constLandau = PolynomialMixed2(energy, fPiZero[3]);
307     mpvLandau   = PolynomialMixed2(energy, fPiZero[4]);
308     sigmaLandau = PolynomialMixed2(energy, fPiZero[5]);
309     
310     break;
311   case 3:
312     
313     constGauss  = PolynomialMixed2(energy, fHadron[0]);
314     meanGauss   = PolynomialMixed2(energy, fHadron[1]);
315     sigmaGauss  = PolynomialMixed2(energy, fHadron[2]);
316     constLandau = PolynomialMixed2(energy, fHadron[3]);
317     mpvLandau   = PolynomialMixed2(energy, fHadron[4]);
318     sigmaLandau = PolynomialMixed2(energy, fHadron[5]);
319
320     break;
321   }
322   
323   distributionParam[0] = constGauss;
324   distributionParam[1] = meanGauss;
325   distributionParam[2] = sigmaGauss;
326   distributionParam[3] = constLandau;
327   distributionParam[4] = mpvLandau;
328   distributionParam[5] = sigmaLandau;
329   
330   return distributionParam;
331 }
332
333 //________________________________________________________
334 Double_t AliEMCALPID::DistEnergy(const Double_t energy, const Int_t type) 
335 {
336   //
337   // Compute the values of the weigh for a given energy the parametrised distribution using the data initialised before.
338   //
339   Double_t constante = 0.;
340   Double_t  energyParam;
341   
342   switch (type) {
343     
344   case 1:  
345     constante  = 1.;    
346     break;
347   case 2:
348       constante  = 1.;
349     break;
350   case 3:
351     constante  = PowerExp(energy, fHadronEnergyProb);
352     break;
353   }
354   
355   energyParam = constante;
356   
357   // //   cout << "Weight   " << constante << " for energy  "<< energy<< " GeV "<<  endl;
358   
359   return energyParam;
360 }
361
362
363 //_______________________________________________________
364 Double_t AliEMCALPID::Polynomial(const Double_t x, const Double_t *params) const
365 {
366   //
367   // Compute a polynomial for a given value of 'x'
368   // with the array of parameters passed as the second arg
369   //
370   
371   Double_t y;
372   y  = params[0];
373   y += params[1] * x;
374   y += params[2] * x * x;
375   y += params[3] * x * x * x;
376   y += params[4] * x * x * x * x;
377   y += params[5] * x * x * x * x * x;
378   
379   return y;
380 }
381 //_______________________________________________________
382 Double_t AliEMCALPID::Polynomial0(const Double_t *params) const 
383 {
384   //
385   // Compute a polynomial for a given value of 'x'
386   // with the array of parameters passed as the second arg
387   //
388   
389   Double_t y;
390   y  = params[0];
391   return y;
392 }
393
394 //_______________________________________________________
395 Double_t AliEMCALPID::Polynomialinv(const Double_t x, const Double_t *params) const
396 {
397   //
398   // Compute a polynomial for a given value of 'x'
399   // with the array of parameters passed as the second arg
400   //
401   
402   Double_t y;
403   if(x>0){
404   y  = params[0];
405   y += params[1] / x;
406   y += params[2] / (x * x);
407   y += params[3] / (x * x * x);
408   y += params[4] / (x * x * x * x);
409   y += params[5] / (x * x * x * x * x);
410   }  
411   else
412     y=0.;
413   return y;
414   
415 }
416 //_______________________________________________________
417 Double_t AliEMCALPID::PolynomialMixed1(const Double_t x, const Double_t *params) const 
418 {
419   //
420   // Compute a polynomial for a given value of 'x'
421   // with the array of parameters passed as the second arg
422   //
423   
424   Double_t y;
425   if(x>0){
426     y  = params[0] / x;
427     y += params[1] ;
428     y += params[2] * x ;
429     //   y += params[3] * 0.;
430     //   y += params[4] * 0.;
431     //   y += params[5] * 0.;
432   }  
433   else
434     y=0.;
435   
436   return y;
437   
438 }
439
440 //_______________________________________________________
441 Double_t AliEMCALPID::PolynomialMixed2(const Double_t x, const Double_t *params) const 
442 {
443   //
444   // Compute a polynomial for a given value of 'x'
445   // with the array of parameters passed as the second arg
446   //
447   
448   Double_t y;
449   if(x>0){
450     y  = params[0] / ( x * x);
451     y += params[1] / x;
452     y += params[2] ;
453     y += params[3] * x ;
454     y += params[4] * x * x ;
455     //   y += params[5] * 0.;
456   }  
457   else
458     y=0.;
459   //   cout << "y = " << y << endl;
460   return y;
461   
462 }
463
464 //_______________________________________________________
465 Double_t AliEMCALPID::PowerExp(const Double_t x, const Double_t *params) const 
466 {
467   //
468   // Compute a polynomial for a given value of 'x'
469   // with the array of parameters passed as the second arg
470   // par[0]*TMath::Power(x[0],par[1])
471   // par[0]*TMath::Exp((x[0]-par[1])*par[2]);
472   
473   Double_t y;
474   
475   y  = params[0] *TMath::Power( x,params[1]);
476   y += params[2] *TMath::Exp((x-params[3])*params[4]);
477   
478   return y;
479   
480 }
481
482
483 //_______________________________________________________
484 void AliEMCALPID::InitParameters()
485 {
486   // Initialize PID parameters, depending on the use or not of the reconstructor
487   // and the kind of event type if the reconstructor is not used.
488   //  fWeightHadronEnergy=0.;
489   //  fWeightPiZeroEnergy=0.;
490   //  fWeightGammaEnergy=0.;
491   
492   fPIDWeight[0] = -1;
493   fPIDWeight[1] = -1;
494   fPIDWeight[2] = -1;
495   
496   for(Int_t i=0; i<AliPID::kSPECIESN+1; i++)
497     fPIDFinal[i]= 0;
498   
499   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
500   
501   if(fReconstructor){
502     
503     if(!recParam) {
504       AliFatal("Reconstruction parameters for EMCAL not set!");
505     }
506     else {
507       
508       for(Int_t i=0; i<6; i++){
509         for(Int_t j=0; j<6; j++){
510           fGamma[i][j]       = recParam->GetGamma(i,j);
511           fGamma1to10[i][j]  = recParam->GetGamma1to10(i,j);
512           fHadron[i][j]      = recParam->GetHadron(i,j);
513           fHadron1to10[i][j] = recParam->GetHadron1to10(i,j);
514           fPiZero[i][j]      = recParam->GetPiZero(i,j);
515           
516           
517           //    AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
518           //                    i,j, fGamma[i][j],fPiZero[i][j],fHadron[i][j] ));
519           //    cout << "PID parameters (" << i << " ,"<<j<<") fGamma= "<<  fGamma[i][j]<<" fPi0 ="<<  fPiZero[i][j]<< endl;
520           
521         } // end loop j
522         fHadronEnergyProb[i] = recParam->GetHadronEnergyProb(i);
523         fPiZeroEnergyProb[i] = recParam->GetPiZeroEnergyProb(i);
524         fGammaEnergyProb[i]  = recParam->GetGammaEnergyProb(i);
525       } //end loop i
526       
527       
528     } // end if !recparam 
529     
530   } 
531   
532   else{
533     //   init the parameters here instead of from loading from recparam
534     //   default parameters are PbPb parameters.
535     SetHighFluxParam();
536     
537   }
538   
539 }
540
541
542 //_______________________________________________________
543 void AliEMCALPID::SetLowFluxParam()
544 {
545   
546   // as a first step, all array elements are initialized to 0.0
547   Int_t i, j;
548   
549   for (i = 0; i < 6; i++) {
550     for (j = 0; j < 6; j++) {
551       fGamma[i][j]      = fHadron[i][j] =  fPiZero[i][j] = 0.;
552       fGamma1to10[i][j] = fHadron1to10[i][j] = 0.;
553     }
554        fGammaEnergyProb[i]  =  fGammaEnergyProb[i];
555        fPiZeroEnergyProb[i] = fPiZeroEnergyProb[i];
556        fHadronEnergyProb[i] = fHadronEnergyProb[i];
557   }
558   
559   // New parametrisation for lambda0^2 (=x): f(x) = normLandau*TMath::Landau(x,mpvLandau,widthLandau)+normgaus*TMath::Gaus(x,meangaus,sigmagaus)
560   // See AliEMCALPid (index j) refers to the polynomial parameters of the fit of each parameter vs energy
561   // pp
562
563   // paramtype[0][j] = norm gauss
564   // paramtype[1][j] = mean gaus
565   // paramtype[2][j] = sigma gaus
566   // paramtype[3][j] = norm landau
567   // paramtype[4][j] = mpv landau
568   // paramtype[5][j] = sigma landau
569
570   fGamma[0][0] = -7.656908e-01; 
571   fGamma[0][1] =  2.352536e-01; 
572   fGamma[0][2] =  1.555996e-02;
573   fGamma[0][3] =  2.243525e-04;
574   fGamma[0][4] = -2.560087e-06;
575   
576   fGamma[1][0] =  6.500216e+00;
577   fGamma[1][1] = -2.564958e-01;
578   fGamma[1][2] =  1.967894e-01;
579   fGamma[1][3] = -3.982273e-04;
580   fGamma[1][4] =  2.797737e-06;
581
582   fGamma[2][0] =  2.416489e+00;
583   fGamma[2][1] = -1.601258e-01;
584   fGamma[2][2] =  3.126839e-02;
585   fGamma[2][3] =  3.387532e-04;
586   fGamma[2][4] = -4.089145e-06;
587
588   fGamma[3][0] =  0.;
589   fGamma[3][1] = -2.696008e+00;
590   fGamma[3][2] =  6.920305e-01;
591   fGamma[3][3] = -2.281122e-03;
592   fGamma[3][4] =  0.;
593
594   fGamma[4][0] =  2.281564e-01;
595   fGamma[4][1] = -7.575040e-02;
596   fGamma[4][2] =  3.813423e-01;
597   fGamma[4][3] = -1.243854e-04;
598   fGamma[4][4] =  1.232045e-06;
599
600   fGamma[5][0] = -3.290107e-01;
601   fGamma[5][1] =  3.707545e-02;
602   fGamma[5][2] =  2.917397e-03;
603   fGamma[5][3] =  4.695306e-05;
604   fGamma[5][4] = -3.572981e-07;
605
606   fHadron[0][0] = 9.482243e-01; 
607   fHadron[0][1] =  -2.780896e-01; 
608   fHadron[0][2] =  2.223507e-02;
609   fHadron[0][3] =  7.294263e-04; 
610   fHadron[0][4] =  -5.665872e-06;
611
612   fHadron[1][0] = 0.;
613   fHadron[1][1] = 0.;
614   fHadron[1][2] = 2.483298e-01;
615   fHadron[1][3] = 0.;
616   fHadron[1][4] = 0.;
617
618   fHadron[2][0] = -5.601199e+00; 
619   fHadron[2][1] =  2.097382e+00; 
620   fHadron[2][2] = -2.307965e-01;
621   fHadron[2][3] =  9.206871e-03;
622   fHadron[2][4] = -8.887548e-05;
623  
624   fHadron[3][0] =  6.543101e+00;
625   fHadron[3][1] =  -2.305203e+00;
626   fHadron[3][2] =  2.761673e-01; 
627   fHadron[3][3] = -5.465855e-03;
628   fHadron[3][4] =  2.784329e-05;
629  
630   fHadron[4][0] = -2.443530e+01;
631   fHadron[4][1] =  8.902578e+00 ;
632   fHadron[4][2] = -5.265901e-01;
633   fHadron[4][3] = 2.549111e-02;
634   fHadron[4][4] =  -2.196801e-04; 
635
636   fHadron[5][0] = 2.102007e-01;
637   fHadron[5][1] =  -3.844418e-02;
638   fHadron[5][2] =  1.234682e-01;
639   fHadron[5][3] = -3.866733e-03;
640   fHadron[5][4] = 3.362719e-05 ;
641
642   fPiZero[0][0] =  5.072157e-01;
643   fPiZero[0][1] = -5.352747e-01;
644   fPiZero[0][2] =  8.499259e-02;
645   fPiZero[0][3] = -3.687401e-03;
646   fPiZero[0][4] =  5.482280e-05;
647
648   fPiZero[1][0] =  4.590137e+02; 
649   fPiZero[1][1] = -7.079341e+01;
650   fPiZero[1][2] =  4.990735e+00;
651   fPiZero[1][3] = -1.241302e-01;
652   fPiZero[1][4] =  1.065772e-03;
653
654   fPiZero[2][0] =  1.376415e+02;
655   fPiZero[2][1] = -3.031577e+01;
656   fPiZero[2][2] =  2.474338e+00;
657   fPiZero[2][3] = -6.903410e-02;
658   fPiZero[2][4] =  6.244089e-04;
659
660   fPiZero[3][0] = 0.;
661   fPiZero[3][1] =  1.145983e+00;
662   fPiZero[3][2] = -2.476052e-01;
663   fPiZero[3][3] =  1.367373e-02;
664   fPiZero[3][4] = 0.;
665
666   fPiZero[4][0] = -2.097586e+02;
667   fPiZero[4][1] =  6.300800e+01;
668   fPiZero[4][2] = -4.038906e+00;
669   fPiZero[4][3] =  1.088543e-01;
670   fPiZero[4][4] = -9.362485e-04;
671
672   fPiZero[5][0] = -1.671477e+01; 
673   fPiZero[5][1] =  2.995415e+00;
674   fPiZero[5][2] = -6.040360e-02;
675   fPiZero[5][3] = -6.137459e-04;
676   fPiZero[5][4] =  1.847328e-05;
677   
678   fHadronEnergyProb[0] = 4.767543e-02;
679   fHadronEnergyProb[1] = -1.537523e+00;
680   fHadronEnergyProb[2] = 2.956727e-01;
681   fHadronEnergyProb[3] = -3.051022e+01;
682   fHadronEnergyProb[4] =-6.036931e-02;
683
684   Int_t ii= 0;
685   Int_t jj= 3;
686   AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
687                         ii,jj, fGamma[ii][jj],fPiZero[ii][jj],fHadron[ii][jj] ));
688   //cout << " LowFlux Parameters fGamma [2][2] = " << fGamma[2][2] << endl;
689   //cout << " LowFlux Parameters fHadron [2][2] = " << fHadron[2][2] << endl;
690    
691   // end for proton-proton  
692
693 }
694
695 //_______________________________________________________
696 void AliEMCALPID::SetHighFluxParam()
697 {
698   
699   // as a first step, all array elements are initialized to 0.0
700   Int_t i, j;
701   for (i = 0; i < 6; i++) {
702     for (j = 0; j < 6; j++) {
703       fGamma[i][j]      = fHadron[i][j] = fPiZero[i][j] = 0.;
704       fGamma1to10[i][j] = fHadron1to10[i][j] = 0.;
705     }
706     fGammaEnergyProb[i]  = 0.;
707     fPiZeroEnergyProb[i] = 0.;
708     fHadronEnergyProb[i] = 0.;
709   }
710   
711   // Pb Pb  this goes with inverted landau + gaussian for gammas, landau+gaussian for Pi0 and hadrons
712   
713   fGamma[0][0] = -7.656908e-01; 
714   fGamma[0][1] =  2.352536e-01; 
715   fGamma[0][2] =  1.555996e-02;
716   fGamma[0][3] =  2.243525e-04;
717   fGamma[0][4] = -2.560087e-06;
718   
719   fGamma[1][0] =  6.500216e+00;
720   fGamma[1][1] = -2.564958e-01;
721   fGamma[1][2] =  1.967894e-01;
722   fGamma[1][3] = -3.982273e-04;
723   fGamma[1][4] =  2.797737e-06;
724
725   fGamma[2][0] =  2.416489e+00;
726   fGamma[2][1] = -1.601258e-01;
727   fGamma[2][2] =  3.126839e-02;
728   fGamma[2][3] =  3.387532e-04;
729   fGamma[2][4] = -4.089145e-06;
730  
731   fGamma[3][0] =  0.;
732   fGamma[3][1] = -2.696008e+00;
733   fGamma[3][2] =  6.920305e-01;
734   fGamma[3][3] = -2.281122e-03;
735   fGamma[3][4] =  0.;
736
737   fGamma[4][0] =  2.281564e-01;
738   fGamma[4][1] = -7.575040e-02;
739   fGamma[4][2] =  3.813423e-01;
740   fGamma[4][3] = -1.243854e-04;
741   fGamma[4][4] =  1.232045e-06;
742
743   fGamma[5][0] = -3.290107e-01;
744   fGamma[5][1] =  3.707545e-02;
745   fGamma[5][2] =  2.917397e-03;
746   fGamma[5][3] =  4.695306e-05;
747   fGamma[5][4] = -3.572981e-07;
748    
749   fHadron[0][0] =   1.519112e-01;
750   fHadron[0][1] = -8.267603e-02;
751   fHadron[0][2] =  1.914574e-02;
752   fHadron[0][3] = -2.677921e-04;
753   fHadron[0][4] =  5.447939e-06;
754
755   fHadron[1][0] = 0.;
756   fHadron[1][1] = -7.549870e-02; 
757   fHadron[1][2] = 3.930087e-01;
758   fHadron[1][3] = -2.368500e-03; 
759   fHadron[1][4] = 0.;
760
761   fHadron[2][0] = 0.;
762   fHadron[2][1] =  -2.463152e-02;
763   fHadron[2][2] = 1.349257e-01;
764   fHadron[2][3] = -1.089440e-03;
765   fHadron[2][4] = 0.;
766
767   fHadron[3][0] = 0.;
768   fHadron[3][1] = 5.101560e-01;
769   fHadron[3][2] = 1.458679e-01;
770   fHadron[3][3] = 4.903068e-04;
771   fHadron[3][4] = 0.;
772
773   fHadron[4][0] = 0.;
774   fHadron[4][1] = -6.693943e-03; 
775   fHadron[4][2] =  2.444753e-01;
776   fHadron[4][3] = -5.553749e-05;
777   fHadron[4][4] = 0.;
778
779   fHadron[5][0] = -4.414030e-01;
780   fHadron[5][1] = 2.292277e-01;
781   fHadron[5][2] = -2.433737e-02;
782   fHadron[5][3] =  1.758422e-03;
783   fHadron[5][4] = -3.001493e-05;
784   
785   fPiZero[0][0] =  5.072157e-01;
786   fPiZero[0][1] = -5.352747e-01;
787   fPiZero[0][2] =  8.499259e-02;
788   fPiZero[0][3] = -3.687401e-03;
789   fPiZero[0][4] =  5.482280e-05;
790   
791   fPiZero[1][0] =  4.590137e+02; 
792   fPiZero[1][1] = -7.079341e+01;
793   fPiZero[1][2] =  4.990735e+00;
794   fPiZero[1][3] = -1.241302e-01;
795   fPiZero[1][4] =  1.065772e-03;
796   
797   fPiZero[2][0] =  1.376415e+02;
798   fPiZero[2][1] = -3.031577e+01;
799   fPiZero[2][2] =  2.474338e+00;
800   fPiZero[2][3] = -6.903410e-02;
801   fPiZero[2][4] =  6.244089e-04;
802
803   fPiZero[3][0] = 0.;
804   fPiZero[3][1] =  1.145983e+00;
805   fPiZero[3][2] = -2.476052e-01;
806   fPiZero[3][3] =  1.367373e-02;
807   fPiZero[3][4] = 0.;
808
809   fPiZero[4][0] = -2.097586e+02;
810   fPiZero[4][1] =  6.300800e+01;
811   fPiZero[4][2] = -4.038906e+00;
812   fPiZero[4][3] =  1.088543e-01;
813   fPiZero[4][4] = -9.362485e-04;
814
815   fPiZero[5][0] = -1.671477e+01; 
816   fPiZero[5][1] =  2.995415e+00;
817   fPiZero[5][2] = -6.040360e-02;
818   fPiZero[5][3] = -6.137459e-04;
819   fPiZero[5][4] =  1.847328e-05;
820
821   // those are the High Flux PbPb ones
822   fHadronEnergyProb[0] = 0.;
823   fHadronEnergyProb[1] = 0.;
824   fHadronEnergyProb[2] =  6.188452e-02;
825   fHadronEnergyProb[3] =  2.030230e+00;
826   fHadronEnergyProb[4] = -6.402242e-02;
827
828  Int_t ii= 0;
829  Int_t jj= 3;
830  AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
831                         ii,jj, fGamma[ii][jj],fPiZero[ii][jj],fHadron[ii][jj] ));
832   //cout << " HighFlux Parameters fGamma [2][2] = " << fGamma[2][2] << endl;
833   //cout << " HighFlux Parameters fHadron [2][2] = " << fHadron[2][2] << endl;
834    
835 }