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 }