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