PID parameters tuned for low and high flux environments (Marie Germain)
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Jul 2009 15:24:41 +0000 (15:24 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 14 Jul 2009 15:24:41 +0000 (15:24 +0000)
EMCAL/AliEMCALPID.cxx
EMCAL/AliEMCALPID.h
EMCAL/AliEMCALRecParam.cxx
EMCAL/AliEMCALRecParam.h
EMCAL/macros/RecParamDB/AliEMCALSetRecParamCDB.C
OCDB/EMCAL/Calib/RecoParam/Run0_999999999_v0_s0.root

index a983825..5fad2a0 100644 (file)
  **************************************************************************/
 
 /* $Id$ */
-/* History of cvs commits:
- *
- * $Log$
- * Revision 1.16  2007/11/23 13:39:05  gustavo
- * Track matching and PID parameters added to AliEMCALRecParam
- *
- * Revision 1.15  2007/10/09 08:46:10  hristov
- * The data members fEMCALClusterCluster and fPHOSCluster are removed from AliESDCaloCluster, the fClusterType is used to select PHOS or EMCAL clusters. Changes, needed to use correctly the new AliESDCaloCluster. (Christian)
- *
- * Revision 1.14  2007/07/26 16:54:53  morsch
- * Changes in AliESDEvent fwd declarartions.
- *
- * Revision 1.13  2007/07/11 13:43:29  hristov
- * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
- *
- * Revision 1.12  2007/06/11 20:43:06  hristov
- * Changes required by the updated AliESDCaloCluster (Gustavo)
- *
- * Revision 1.11  2007/03/30 13:50:34  gustavo
- * PID for particles with E < 5 GeV was not done, temporal solution found (Guenole)
- *
- * Revision 1.10  2007/03/09 14:34:11  gustavo
- * Correct probability calculation, added missing initialization of data members
- *
- * Revision 1.9  2007/02/20 20:17:43  hristov
- * Corrected array size, removed warnings (icc)
- *
- * Revision 1.8  2006/12/19 08:49:35  gustavo
- * New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled when calling AliEMCALPID in AliEMCALReconstructor::FillESD()
- *
- *
- */
-//    to compute PID for all the clusters in ESDs.root file
-//     the ESDs.root have to be in the same directory as the class
-//
-//
-//
-//
+
+//   Compute PID weights for all the clusters that are in AliESDs.root file
+//   the AliESDs.root have to be in the same directory as the class
 //
-//     AliEMCALPID::CalculPID(Energy,Lambda0)
-//       Calcul PID for all clusters in AliESDs.root file
+//   and do:    
+//   AliEMCALPID *pid = new AliEMCALPID(kFALSE); // this calls the constructor which avoids the call to recparam 
+//   pid->SetReconstructor(kFALSE);
+//   pid->SetPrintInfo(kTRUE);
+//   pid->SetHighFluxParam(); //   pid->SetLowFluxParam(); 
+//   
+//   then in cluster loop do
+//   pid->ComputePID(energy, lambda0);
+//       
+//        Compute PID Weight for all clusters in AliESDs.root file
 //       keep this function for the moment for a simple verification, could be removed
 //
+//   pid->GetPIDFinal(idx) gives the probabilities
 //
-//
-//   AliEMCALPID::CalculPID(Energy,Lambda0)
-//    calcul PID Weght for a cluster with Energy, Lambda0 .
-//    Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
+//   Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
 //
 //
 //
 //                   Pi0  PID[1]
 //                Hadron  PID[2]
 //
-//  
-//
-//
-//
-// --- ROOT system ---
+// --- standard c ---
 
 // standard C++ includes
-#include <Riostream.h>
+//#include <Riostream.h>
 
 // ROOT includes
-#include "TTree.h"
-#include "TStyle.h"
-#include "TVector3.h"
-#include "TBranch.h"
-#include "TClonesArray.h"
-#include "TCanvas.h"
-#include "TLorentzVector.h"
+//#include "TTree.h"
+//#include "TVector3.h"
+//#include "TBranch.h"
+//#include "TClonesArray.h"
+//#include "TLorentzVector.h"
 #include "TMath.h"
-#include "TFile.h"
-#include "TH1.h"
-#include "TH2.h"
-#include "TParticle.h"
+//#include "TRefArray.h"
+#include "TArrayD.h"
 
 // STEER includes
-#include "AliLog.h"
+#include "AliESDEvent.h"
+//#include "AliLog.h"
 #include "AliEMCALPID.h"
 #include "AliESDCaloCluster.h"
-#include "AliEMCALRecParam.h"
+//#include "AliEMCALRecParam.h"
 #include "AliEMCALReconstructor.h"
 
   
 ClassImp(AliEMCALPID)
-
+  
 //______________________________________________
   AliEMCALPID::AliEMCALPID():
-    fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE)
+    fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.), fWeightHadronEnergy(1.), fWeightGammaEnergy(1.),fWeightPiZeroEnergy(1.),fReconstructor(kTRUE)
 {
   //
   // Constructor.
   // Initialize all constant values which have to be used
   // during PID algorithm execution
   //
-  fPIDWeight[0] = -1;
-  fPIDWeight[1] = -1;
-  fPIDWeight[2] = -1;
-
-  for(Int_t i=0; i<AliPID::kSPECIESN+1; i++)
-    fPIDFinal[i]= 0;
+  
+  InitParameters(); 
+  
+  
+}
 
-  const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
-  if(!recParam) {
-    AliFatal("Reconstruction parameters for EMCAL not set!");
-  }
-  else {
-    for(Int_t i=0; i<6; i++){
-      for(Int_t j=0; j<6; j++){
-       fGamma[i][j] = recParam->GetGamma(i,j);
-       fHadron[i][j] = recParam->GetHadron(i,j);
-       fPiZero5to10[i][j] = recParam->GetPiZero5to10(i,j);
-       fPiZero10to60[i][j] = recParam->GetPiZero10to60(i,j);
-       AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
-                       i,j, fGamma[i][j],fPiZero5to10[i][j],fHadron[i][j] ));
-      }
-    }
-    
-  }
+//______________________________________________
+AliEMCALPID::AliEMCALPID(Bool_t reconstructor):
+  fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.), fWeightHadronEnergy(1.), fWeightGammaEnergy(1.),fWeightPiZeroEnergy(1.),fReconstructor(reconstructor)
+{
+  //
+  // Constructor.
+  // Initialize all constant values which have to be used
+  // during PID algorithm execution called when used in standalone mode 
+  //
+  
+  InitParameters(); 
   
 }
 
 //______________________________________________
 void AliEMCALPID::RunPID(AliESDEvent *esd)
 {
-//
-// Make the PID for all the EMCAL clusters containedin the ESDs File
-// but just gamma/PiO/Hadron
-//
-       // trivial check against NULL object passed
-
+  //
+  // Make the PID for all the EMCAL clusters containedin the ESDs File
+  // but just gamma/PiO/Hadron
+  //
+  // trivial check against NULL object passed
+  
   if (esd == 0x0) {
     AliInfo("NULL ESD object passed !!" );
     return ;
   }
-
+  
   Int_t nClusters = esd->GetNumberOfCaloClusters();
   Int_t firstCluster = 0;
   Double_t energy, lambda0;
   for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
-
+    
     AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
-       if (!clust->IsEMCAL()) continue ; 
+    if (!clust->IsEMCAL()) continue ; 
+    
     energy = clust->E();
     lambda0 = clust->GetM02();
     // verify cluster type
     Int_t clusterType= clust->GetClusterType();
     if (clusterType == AliESDCaloCluster::kEMCALClusterv1 && lambda0 != 0  && energy < 1000) {
-
-
+      
+      //      if (lambda0 != 0  && energy < 1000) {
+      
       // reject clusters with lambda0 = 0
-
-
+      
+      
       ComputePID(energy, lambda0);
-
-
+      
+      
       if (fPrintInfo) {
        AliInfo("___________________________________________________");
        AliInfo(Form( "Particle Energy = %f",energy));
@@ -209,9 +167,10 @@ void AliEMCALPID::RunPID(AliESDEvent *esd)
        AliInfo(Form( " kUnknown  : %f", fPIDFinal[10] ));
        AliInfo("___________________________________________________");
       }
-
-      if(fReconstructor) // In case it is called during reconstruction.
-       clust->SetPid(fPIDFinal);
+      
+      if(fReconstructor){ // In case it is called during reconstruction.
+       //      cout << "#############On remplit l esd avec les PIDWeight##########" << endl;
+       clust->SetPid(fPIDFinal);}
     } // end if (clusterType...)
   } // end for (iCluster...)
 }
@@ -223,60 +182,102 @@ void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
 // This is the main command, which uses the distributions computed and parametrised, 
 // and gives the PID by the bayesian method.
 //
-
-if (energy<5){energy =6;}
-
-
+//   cout << "ENERGY  " <<energy <<" lambda0 "<< lambda0<<  endl;
+  
+  Double_t weightGammaEnergy  = DistEnergy(energy, 1);
+  Double_t weightPiZeroEnergy = DistEnergy(energy, 2);
+  Double_t weightHadronEnergy = DistEnergy(energy, 3);
+  
+  //Double_t weightHadronEnergy = 1.;
+  
+  Double_t energyhadron=energy;
+  if(energyhadron<1.)energyhadron=1.; // no energy dependance of  parametrisation for hadrons below 1 GeV
+  if (energy<2){energy =2;} // no energy dependance of parametrisation for gamma and pi0 below 2 GeV
+  
+  if (energy>55){
+    energy =55.;
+    energyhadron=55.;
+  } // same parametrisation for gamma and hadrons above 55 GeV 
+  //   for the pi0 above 55GeV the 2 gammas supperposed no way to distinguish from real gamma  PIDWeight[1]=0
+  
   TArrayD paramDistribGamma  = DistLambda0(energy, 1);
   TArrayD paramDistribPiZero = DistLambda0(energy, 2);
-  TArrayD paramDistribHadron = DistLambda0(energy, 3);
+  TArrayD paramDistribHadron = DistLambda0(energyhadron, 3);
   
   Bool_t norm = kFALSE;
   
+  
   fProbGamma   = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0];
-  fProbGamma  += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3];
-  fProbPiZero  = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
-  fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
+  fProbGamma  += TMath::Landau(((1-paramDistribGamma[4])-lambda0),paramDistribGamma[4],paramDistribGamma[5],norm)* paramDistribGamma[3];
+  if(fProbGamma<0.)fProbGamma=0.;
+  
+  fProbGamma = fProbGamma*weightGammaEnergy;
+  
+  if(energy>10. || energy < 55.){
+    fProbPiZero  = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
+    fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
+    if(fProbPiZero<0. || energy<5.)fProbPiZero=0.;
+    fProbPiZero = fProbPiZero*weightPiZeroEnergy;
+  }
+  else {
+    fProbPiZero = 0.;
+  }
+  
   fProbHadron  = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
   fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
+  if(fProbHadron<0.)fProbHadron=0.;
+  fProbHadron = fProbHadron*weightHadronEnergy; // to take into account the probability for a hadron to have a given reconstructed energy 
   
   // compute PID Weight
-  fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
-  fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
-  fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
+  if( (fProbGamma + fProbPiZero + fProbHadron)>0.){
+    fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
+    fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
+    fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
+  }
+  else{   
+// cases where  energy and lambda0 large,  probably du to 2 clusters folded the clusters PID not assigned to hadron nor Pi0 nor gammas
+    fPIDWeight[0] = 0.;
+    fPIDWeight[1] = 0.;
+    fPIDWeight[2] = 0.;
+  }
+  
+  
+  // cout << " PID[0] "<<  fPIDWeight[0] <<  " PID[1] "<<  fPIDWeight[1] <<  " PID[2] "<<  fPIDWeight[2] << endl;
   
   SetPID(fPIDWeight[0], 0);
   SetPID(fPIDWeight[1], 1);
   SetPID(fPIDWeight[2], 2);
   
-  // sortie ecran pid Weight only for control (= in english ???)
+  // print  pid Weight only for control (= in english ???)
   if (fPrintInfo) {
     AliInfo(Form( "Energy in loop = %f", energy) );
     AliInfo(Form( "Lambda0 in loop = %f", lambda0) );
     AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) );
-    // AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) );
     AliInfo(Form( "fProbaPiZero = %f", fProbPiZero ));
     AliInfo(Form( "fProbaHadron = %f", fProbHadron) );
     AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f",  fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
-    AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] ));
     AliInfo("********************************************************" );
   }
   
-  fPIDFinal[0]  = fPIDWeight[0]/2;
+  fPIDFinal[0]  = fPIDWeight[0]/2; // photon
   fPIDFinal[1]  = fPIDWeight[2]/8;
   fPIDFinal[2]  = fPIDWeight[2]/8;
   fPIDFinal[3]  = fPIDWeight[2]/8;
   fPIDFinal[4]  = fPIDWeight[2]/8;
-  fPIDFinal[5]  = fPIDWeight[0]/2;
-  fPIDFinal[6]  = fPIDWeight[1]  ;
+  fPIDFinal[5]  = fPIDWeight[0]/2; // electron
+  fPIDFinal[6]  = fPIDWeight[1]  ; // Pi0
   fPIDFinal[7]  = fPIDWeight[2]/8;
   fPIDFinal[8]  = fPIDWeight[2]/8;
   fPIDFinal[9]  = fPIDWeight[2]/8;
   fPIDFinal[10] = fPIDWeight[2]/8;
+
 }
 
+
+
+
 //________________________________________________________
-TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
+TArrayD AliEMCALPID::DistLambda0(const Double_t energy, const Int_t type) 
 {
   //
   // Compute the values of the parametrised distributions using the data initialised before.
@@ -286,43 +287,35 @@ TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
   TArrayD  distributionParam(6);
   
   switch (type) {
+    
   case 1:
-    constGauss  = Polynomial(energy, fGamma[0]);
-    meanGauss   = Polynomial(energy, fGamma[1]);
-    sigmaGauss  = Polynomial(energy, fGamma[2]);
-    constLandau = Polynomial(energy, fGamma[3]);
-    mpvLandau   = Polynomial(energy, fGamma[4]);
-    sigmaLandau = Polynomial(energy, fGamma[5]);
     
+    constGauss  = PolynomialMixed2(energy, fGamma[0]);
+    meanGauss   = PolynomialMixed2(energy, fGamma[1]);
+    sigmaGauss  = PolynomialMixed2(energy, fGamma[2]);
+    constLandau = PolynomialMixed2(energy, fGamma[3]);
+    mpvLandau   = PolynomialMixed2(energy, fGamma[4]);
+    sigmaLandau = PolynomialMixed2(energy, fGamma[5]);
+   break;
 
-    break;
   case 2:
-    if (energy < 10) {
-      constGauss  = Polynomial(energy, fPiZero5to10[0]);
-      meanGauss   = Polynomial(energy, fPiZero5to10[1]);
-      sigmaGauss  = Polynomial(energy, fPiZero5to10[2]);
-      constLandau = Polynomial(energy, fPiZero5to10[3]);
-      mpvLandau   = Polynomial(energy, fPiZero5to10[4]);
-      sigmaLandau = Polynomial(energy, fPiZero5to10[5]);
 
-    }
-    else {
-      constGauss  = Polynomial(energy, fPiZero10to60[0]);
-      meanGauss   = Polynomial(energy, fPiZero10to60[1]);
-      sigmaGauss  = Polynomial(energy, fPiZero10to60[2]);
-      constLandau = Polynomial(energy, fPiZero10to60[3]);
-      mpvLandau   = Polynomial(energy, fPiZero10to60[4]);
-      sigmaLandau = Polynomial(energy, fPiZero10to60[5]);
-    }
+    constGauss  = PolynomialMixed2(energy, fPiZero[0]);
+    meanGauss   = PolynomialMixed2(energy, fPiZero[1]);
+    sigmaGauss  = PolynomialMixed2(energy, fPiZero[2]);
+    constLandau = PolynomialMixed2(energy, fPiZero[3]);
+    mpvLandau   = PolynomialMixed2(energy, fPiZero[4]);
+    sigmaLandau = PolynomialMixed2(energy, fPiZero[5]);
+    
     break;
   case 3:
-    constGauss  = Polynomial(energy, fHadron[0]);
-    meanGauss   = Polynomial(energy, fHadron[1]);
-    sigmaGauss  = Polynomial(energy, fHadron[2]);
-    constLandau = Polynomial(energy, fHadron[3]);
-    mpvLandau   = Polynomial(energy, fHadron[4]);
-    sigmaLandau = Polynomial(energy, fHadron[5]);
+    
+    constGauss  = PolynomialMixed2(energy, fHadron[0]);
+    meanGauss   = PolynomialMixed2(energy, fHadron[1]);
+    sigmaGauss  = PolynomialMixed2(energy, fHadron[2]);
+    constLandau = PolynomialMixed2(energy, fHadron[3]);
+    mpvLandau   = PolynomialMixed2(energy, fHadron[4]);
+    sigmaLandau = PolynomialMixed2(energy, fHadron[5]);
 
     break;
   }
@@ -337,8 +330,38 @@ TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
   return distributionParam;
 }
 
+//________________________________________________________
+Double_t AliEMCALPID::DistEnergy(const Double_t energy, const Int_t type) 
+{
+  //
+  // Compute the values of the weigh for a given energy the parametrised distribution using the data initialised before.
+  //
+  Double_t constante = 0.;
+  Double_t  energyParam;
+  
+  switch (type) {
+    
+  case 1:  
+    constante  = 1.;    
+    break;
+  case 2:
+      constante  = 1.;
+    break;
+  case 3:
+    constante  = PowerExp(energy, fHadronEnergyProb);
+    break;
+  }
+  
+  energyParam = constante;
+  
+  // //   cout << "Weight   " << constante << " for energy  "<< energy<< " GeV "<<  endl;
+  
+  return energyParam;
+}
+
+
 //_______________________________________________________
-Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
+Double_t AliEMCALPID::Polynomial(const Double_t x, const Double_t *params) const
 {
   //
   // Compute a polynomial for a given value of 'x'
@@ -355,3 +378,458 @@ Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
   
   return y;
 }
+//_______________________________________________________
+Double_t AliEMCALPID::Polynomial0(const Double_t *params) const 
+{
+  //
+  // Compute a polynomial for a given value of 'x'
+  // with the array of parameters passed as the second arg
+  //
+  
+  Double_t y;
+  y  = params[0];
+  return y;
+}
+
+//_______________________________________________________
+Double_t AliEMCALPID::Polynomialinv(const Double_t x, const Double_t *params) const
+{
+  //
+  // Compute a polynomial for a given value of 'x'
+  // with the array of parameters passed as the second arg
+  //
+  
+  Double_t y;
+  if(x>0){
+  y  = params[0];
+  y += params[1] / x;
+  y += params[2] / (x * x);
+  y += params[3] / (x * x * x);
+  y += params[4] / (x * x * x * x);
+  y += params[5] / (x * x * x * x * x);
+  }  
+  else
+    y=0.;
+  return y;
+  
+}
+//_______________________________________________________
+Double_t AliEMCALPID::PolynomialMixed1(const Double_t x, const Double_t *params) const 
+{
+  //
+  // Compute a polynomial for a given value of 'x'
+  // with the array of parameters passed as the second arg
+  //
+  
+  Double_t y;
+  if(x>0){
+    y  = params[0] / x;
+    y += params[1] ;
+    y += params[2] * x ;
+    //   y += params[3] * 0.;
+    //   y += params[4] * 0.;
+    //   y += params[5] * 0.;
+  }  
+  else
+    y=0.;
+  
+  return y;
+  
+}
+
+//_______________________________________________________
+Double_t AliEMCALPID::PolynomialMixed2(const Double_t x, const Double_t *params) const 
+{
+  //
+  // Compute a polynomial for a given value of 'x'
+  // with the array of parameters passed as the second arg
+  //
+  
+  Double_t y;
+  if(x>0){
+    y  = params[0] / ( x * x);
+    y += params[1] / x;
+    y += params[2] ;
+    y += params[3] * x ;
+    y += params[4] * x * x ;
+    //   y += params[5] * 0.;
+  }  
+  else
+    y=0.;
+  //   cout << "y = " << y << endl;
+  return y;
+  
+}
+
+//_______________________________________________________
+Double_t AliEMCALPID::PowerExp(const Double_t x, const Double_t *params) const 
+{
+  //
+  // Compute a polynomial for a given value of 'x'
+  // with the array of parameters passed as the second arg
+  // par[0]*TMath::Power(x[0],par[1])
+  // par[0]*TMath::Exp((x[0]-par[1])*par[2]);
+  
+  Double_t y;
+  
+  y  = params[0] *TMath::Power( x,params[1]);
+  y += params[2] *TMath::Exp((x-params[3])*params[4]);
+  
+  return y;
+  
+}
+
+
+//_______________________________________________________
+void AliEMCALPID::InitParameters()
+{
+  // Initialize PID parameters, depending on the use or not of the reconstructor
+  // and the kind of event type if the reconstructor is not used.
+  //  fWeightHadronEnergy=0.;
+  //  fWeightPiZeroEnergy=0.;
+  //  fWeightGammaEnergy=0.;
+  
+  fPIDWeight[0] = -1;
+  fPIDWeight[1] = -1;
+  fPIDWeight[2] = -1;
+  
+  for(Int_t i=0; i<AliPID::kSPECIESN+1; i++)
+    fPIDFinal[i]= 0;
+  
+  const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
+  
+  if(fReconstructor){
+    
+    if(!recParam) {
+      AliFatal("Reconstruction parameters for EMCAL not set!");
+    }
+    else {
+      
+      for(Int_t i=0; i<6; i++){
+       for(Int_t j=0; j<6; j++){
+         fGamma[i][j]       = recParam->GetGamma(i,j);
+         fGamma1to10[i][j]  = recParam->GetGamma1to10(i,j);
+         fHadron[i][j]      = recParam->GetHadron(i,j);
+         fHadron1to10[i][j] = recParam->GetHadron1to10(i,j);
+         fPiZero[i][j]      = recParam->GetPiZero(i,j);
+         
+         
+         //    AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
+         //                    i,j, fGamma[i][j],fPiZero[i][j],fHadron[i][j] ));
+         //    cout << "PID parameters (" << i << " ,"<<j<<") fGamma= "<<  fGamma[i][j]<<" fPi0 ="<<  fPiZero[i][j]<< endl;
+         
+       } // end loop j
+       fHadronEnergyProb[i] = recParam->GetHadronEnergyProb(i);
+       fPiZeroEnergyProb[i] = recParam->GetPiZeroEnergyProb(i);
+       fGammaEnergyProb[i]  = recParam->GetGammaEnergyProb(i);
+      } //end loop i
+      
+      
+    } // end if !recparam 
+    
+  } 
+  
+  else{
+    //   init the parameters here instead of from loading from recparam
+    //   default parameters are PbPb parameters.
+    SetHighFluxParam();
+    
+  }
+  
+}
+
+
+//_______________________________________________________
+void AliEMCALPID::SetLowFluxParam()
+{
+  
+  // as a first step, all array elements are initialized to 0.0
+  Int_t i, j;
+  
+  for (i = 0; i < 6; i++) {
+    for (j = 0; j < 6; j++) {
+      fGamma[i][j]      = fHadron[i][j] =  fPiZero[i][j] = 0.;
+      fGamma1to10[i][j] = fHadron1to10[i][j] = 0.;
+    }
+       fGammaEnergyProb[i]  =  fGammaEnergyProb[i];
+       fPiZeroEnergyProb[i] = fPiZeroEnergyProb[i];
+       fHadronEnergyProb[i] = fHadronEnergyProb[i];
+  }
+  
+  // New parametrisation for lambda0^2 (=x): f(x) = normLandau*TMath::Landau(x,mpvLandau,widthLandau)+normgaus*TMath::Gaus(x,meangaus,sigmagaus)
+  // See AliEMCALPid (index j) refers to the polynomial parameters of the fit of each parameter vs energy
+  // pp
+
+  // paramtype[0][j] = norm gauss
+  // paramtype[1][j] = mean gaus
+  // paramtype[2][j] = sigma gaus
+  // paramtype[3][j] = norm landau
+  // paramtype[4][j] = mpv landau
+  // paramtype[5][j] = sigma landau
+
+  fGamma[0][0] = -7.656908e-01; 
+  fGamma[0][1] =  2.352536e-01; 
+  fGamma[0][2] =  1.555996e-02;
+  fGamma[0][3] =  2.243525e-04;
+  fGamma[0][4] = -2.560087e-06;
+  
+  fGamma[1][0] =  6.500216e+00;
+  fGamma[1][1] = -2.564958e-01;
+  fGamma[1][2] =  1.967894e-01;
+  fGamma[1][3] = -3.982273e-04;
+  fGamma[1][4] =  2.797737e-06;
+
+  fGamma[2][0] =  2.416489e+00;
+  fGamma[2][1] = -1.601258e-01;
+  fGamma[2][2] =  3.126839e-02;
+  fGamma[2][3] =  3.387532e-04;
+  fGamma[2][4] = -4.089145e-06;
+
+  fGamma[3][0] =  0.;
+  fGamma[3][1] = -2.696008e+00;
+  fGamma[3][2] =  6.920305e-01;
+  fGamma[3][3] = -2.281122e-03;
+  fGamma[3][4] =  0.;
+
+  fGamma[4][0] =  2.281564e-01;
+  fGamma[4][1] = -7.575040e-02;
+  fGamma[4][2] =  3.813423e-01;
+  fGamma[4][3] = -1.243854e-04;
+  fGamma[4][4] =  1.232045e-06;
+
+  fGamma[5][0] = -3.290107e-01;
+  fGamma[5][1] =  3.707545e-02;
+  fGamma[5][2] =  2.917397e-03;
+  fGamma[5][3] =  4.695306e-05;
+  fGamma[5][4] = -3.572981e-07;
+
+  fHadron[0][0] = 9.482243e-01; 
+  fHadron[0][1] =  -2.780896e-01; 
+  fHadron[0][2] =  2.223507e-02;
+  fHadron[0][3] =  7.294263e-04; 
+  fHadron[0][4] =  -5.665872e-06;
+
+  fHadron[1][0] = 0.;
+  fHadron[1][1] = 0.;
+  fHadron[1][2] = 2.483298e-01;
+  fHadron[1][3] = 0.;
+  fHadron[1][4] = 0.;
+
+  fHadron[2][0] = -5.601199e+00; 
+  fHadron[2][1] =  2.097382e+00; 
+  fHadron[2][2] = -2.307965e-01;
+  fHadron[2][3] =  9.206871e-03;
+  fHadron[2][4] = -8.887548e-05;
+  fHadron[3][0] =  6.543101e+00;
+  fHadron[3][1] =  -2.305203e+00;
+  fHadron[3][2] =  2.761673e-01; 
+  fHadron[3][3] = -5.465855e-03;
+  fHadron[3][4] =  2.784329e-05;
+  fHadron[4][0] = -2.443530e+01;
+  fHadron[4][1] =  8.902578e+00 ;
+  fHadron[4][2] = -5.265901e-01;
+  fHadron[4][3] = 2.549111e-02;
+  fHadron[4][4] =  -2.196801e-04; 
+
+  fHadron[5][0] = 2.102007e-01;
+  fHadron[5][1] =  -3.844418e-02;
+  fHadron[5][2] =  1.234682e-01;
+  fHadron[5][3] = -3.866733e-03;
+  fHadron[5][4] = 3.362719e-05 ;
+
+  fPiZero[0][0] =  5.072157e-01;
+  fPiZero[0][1] = -5.352747e-01;
+  fPiZero[0][2] =  8.499259e-02;
+  fPiZero[0][3] = -3.687401e-03;
+  fPiZero[0][4] =  5.482280e-05;
+
+  fPiZero[1][0] =  4.590137e+02; 
+  fPiZero[1][1] = -7.079341e+01;
+  fPiZero[1][2] =  4.990735e+00;
+  fPiZero[1][3] = -1.241302e-01;
+  fPiZero[1][4] =  1.065772e-03;
+
+  fPiZero[2][0] =  1.376415e+02;
+  fPiZero[2][1] = -3.031577e+01;
+  fPiZero[2][2] =  2.474338e+00;
+  fPiZero[2][3] = -6.903410e-02;
+  fPiZero[2][4] =  6.244089e-04;
+
+  fPiZero[3][0] = 0.;
+  fPiZero[3][1] =  1.145983e+00;
+  fPiZero[3][2] = -2.476052e-01;
+  fPiZero[3][3] =  1.367373e-02;
+  fPiZero[3][4] = 0.;
+
+  fPiZero[4][0] = -2.097586e+02;
+  fPiZero[4][1] =  6.300800e+01;
+  fPiZero[4][2] = -4.038906e+00;
+  fPiZero[4][3] =  1.088543e-01;
+  fPiZero[4][4] = -9.362485e-04;
+
+  fPiZero[5][0] = -1.671477e+01; 
+  fPiZero[5][1] =  2.995415e+00;
+  fPiZero[5][2] = -6.040360e-02;
+  fPiZero[5][3] = -6.137459e-04;
+  fPiZero[5][4] =  1.847328e-05;
+  
+  fHadronEnergyProb[0] = 4.767543e-02;
+  fHadronEnergyProb[1] = -1.537523e+00;
+  fHadronEnergyProb[2] = 2.956727e-01;
+  fHadronEnergyProb[3] = -3.051022e+01;
+  fHadronEnergyProb[4] =-6.036931e-02;
+
+  Int_t ii= 0;
+  Int_t jj= 3;
+  AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
+                       ii,jj, fGamma[ii][jj],fPiZero[ii][jj],fHadron[ii][jj] ));
+  //cout << " LowFlux Parameters fGamma [2][2] = " << fGamma[2][2] << endl;
+  //cout << " LowFlux Parameters fHadron [2][2] = " << fHadron[2][2] << endl;
+   
+  // end for proton-proton  
+
+}
+
+//_______________________________________________________
+void AliEMCALPID::SetHighFluxParam()
+{
+  
+  // as a first step, all array elements are initialized to 0.0
+  Int_t i, j;
+  for (i = 0; i < 6; i++) {
+    for (j = 0; j < 6; j++) {
+      fGamma[i][j]      = fHadron[i][j] = fPiZero[i][j] = 0.;
+      fGamma1to10[i][j] = fHadron1to10[i][j] = 0.;
+    }
+    fGammaEnergyProb[i]  = 0.;
+    fPiZeroEnergyProb[i] = 0.;
+    fHadronEnergyProb[i] = 0.;
+  }
+  
+  // Pb Pb  this goes with inverted landau + gaussian for gammas, landau+gaussian for Pi0 and hadrons
+  
+  fGamma[0][0] = -7.656908e-01; 
+  fGamma[0][1] =  2.352536e-01; 
+  fGamma[0][2] =  1.555996e-02;
+  fGamma[0][3] =  2.243525e-04;
+  fGamma[0][4] = -2.560087e-06;
+  
+  fGamma[1][0] =  6.500216e+00;
+  fGamma[1][1] = -2.564958e-01;
+  fGamma[1][2] =  1.967894e-01;
+  fGamma[1][3] = -3.982273e-04;
+  fGamma[1][4] =  2.797737e-06;
+
+  fGamma[2][0] =  2.416489e+00;
+  fGamma[2][1] = -1.601258e-01;
+  fGamma[2][2] =  3.126839e-02;
+  fGamma[2][3] =  3.387532e-04;
+  fGamma[2][4] = -4.089145e-06;
+  fGamma[3][0] =  0.;
+  fGamma[3][1] = -2.696008e+00;
+  fGamma[3][2] =  6.920305e-01;
+  fGamma[3][3] = -2.281122e-03;
+  fGamma[3][4] =  0.;
+
+  fGamma[4][0] =  2.281564e-01;
+  fGamma[4][1] = -7.575040e-02;
+  fGamma[4][2] =  3.813423e-01;
+  fGamma[4][3] = -1.243854e-04;
+  fGamma[4][4] =  1.232045e-06;
+
+  fGamma[5][0] = -3.290107e-01;
+  fGamma[5][1] =  3.707545e-02;
+  fGamma[5][2] =  2.917397e-03;
+  fGamma[5][3] =  4.695306e-05;
+  fGamma[5][4] = -3.572981e-07;
+   
+  fHadron[0][0] =   1.519112e-01;
+  fHadron[0][1] = -8.267603e-02;
+  fHadron[0][2] =  1.914574e-02;
+  fHadron[0][3] = -2.677921e-04;
+  fHadron[0][4] =  5.447939e-06;
+
+  fHadron[1][0] = 0.;
+  fHadron[1][1] = -7.549870e-02; 
+  fHadron[1][2] = 3.930087e-01;
+  fHadron[1][3] = -2.368500e-03; 
+  fHadron[1][4] = 0.;
+
+  fHadron[2][0] = 0.;
+  fHadron[2][1] =  -2.463152e-02;
+  fHadron[2][2] = 1.349257e-01;
+  fHadron[2][3] = -1.089440e-03;
+  fHadron[2][4] = 0.;
+
+  fHadron[3][0] = 0.;
+  fHadron[3][1] = 5.101560e-01;
+  fHadron[3][2] = 1.458679e-01;
+  fHadron[3][3] = 4.903068e-04;
+  fHadron[3][4] = 0.;
+
+  fHadron[4][0] = 0.;
+  fHadron[4][1] = -6.693943e-03; 
+  fHadron[4][2] =  2.444753e-01;
+  fHadron[4][3] = -5.553749e-05;
+  fHadron[4][4] = 0.;
+
+  fHadron[5][0] = -4.414030e-01;
+  fHadron[5][1] = 2.292277e-01;
+  fHadron[5][2] = -2.433737e-02;
+  fHadron[5][3] =  1.758422e-03;
+  fHadron[5][4] = -3.001493e-05;
+  
+  fPiZero[0][0] =  5.072157e-01;
+  fPiZero[0][1] = -5.352747e-01;
+  fPiZero[0][2] =  8.499259e-02;
+  fPiZero[0][3] = -3.687401e-03;
+  fPiZero[0][4] =  5.482280e-05;
+  
+  fPiZero[1][0] =  4.590137e+02; 
+  fPiZero[1][1] = -7.079341e+01;
+  fPiZero[1][2] =  4.990735e+00;
+  fPiZero[1][3] = -1.241302e-01;
+  fPiZero[1][4] =  1.065772e-03;
+  
+  fPiZero[2][0] =  1.376415e+02;
+  fPiZero[2][1] = -3.031577e+01;
+  fPiZero[2][2] =  2.474338e+00;
+  fPiZero[2][3] = -6.903410e-02;
+  fPiZero[2][4] =  6.244089e-04;
+
+  fPiZero[3][0] = 0.;
+  fPiZero[3][1] =  1.145983e+00;
+  fPiZero[3][2] = -2.476052e-01;
+  fPiZero[3][3] =  1.367373e-02;
+  fPiZero[3][4] = 0.;
+
+  fPiZero[4][0] = -2.097586e+02;
+  fPiZero[4][1] =  6.300800e+01;
+  fPiZero[4][2] = -4.038906e+00;
+  fPiZero[4][3] =  1.088543e-01;
+  fPiZero[4][4] = -9.362485e-04;
+
+  fPiZero[5][0] = -1.671477e+01; 
+  fPiZero[5][1] =  2.995415e+00;
+  fPiZero[5][2] = -6.040360e-02;
+  fPiZero[5][3] = -6.137459e-04;
+  fPiZero[5][4] =  1.847328e-05;
+
+  // those are the High Flux PbPb ones
+  fHadronEnergyProb[0] = 0.;
+  fHadronEnergyProb[1] = 0.;
+  fHadronEnergyProb[2] =  6.188452e-02;
+  fHadronEnergyProb[3] =  2.030230e+00;
+  fHadronEnergyProb[4] = -6.402242e-02;
+
+ Int_t ii= 0;
+ Int_t jj= 3;
+ AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
+                       ii,jj, fGamma[ii][jj],fPiZero[ii][jj],fHadron[ii][jj] ));
+  //cout << " HighFlux Parameters fGamma [2][2] = " << fGamma[2][2] << endl;
+  //cout << " HighFlux Parameters fHadron [2][2] = " << fHadron[2][2] << endl;
+   
+}
index b4fcf06..d649371 100644 (file)
@@ -1,29 +1,19 @@
-#ifndef AliEMCALPID_H
-#define AliEMCALPID_H
+#ifndef ALIEMCALPID_H
+#define ALIEMCALPID_H
 
 /* $Id$ */
-/* History of cvs commits:
- *
- * $Log$
- * Revision 1.13  2007/07/11 13:43:29  hristov
- * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
- *
- * Revision 1.12  2007/02/20 20:17:43  hristov
- * Corrected array size, removed warnings (icc)
- *
- * Revision 1.11  2006/12/19 08:49:35  gustavo
- * New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled when calling AliEMCALPID in AliEMCALReconstructor::FillESD()
- *
- *
- */
 
 ///////////////////////////////////////////////////////////////////////////////
 // Class AliEMCALPID
+// Compute PID weights for all the clusters
 ///////////////////////////////////////////////////////////////////////////////
 
+//Root includes
 #include "TTask.h"
-#include "TArrayD.h"
-#include "AliESDEvent.h"
+//#include "TArrayD.h"
+class TArrayD ;
+//AliRoot includes
+class AliESDEvent ;
 #include "AliPID.h" 
 
 class AliEMCALPID : public TTask {
@@ -31,42 +21,64 @@ class AliEMCALPID : public TTask {
 public:
   
   AliEMCALPID();
+  AliEMCALPID(Bool_t reconstructor);
   virtual ~AliEMCALPID() { }
   
   void     RunPID(AliESDEvent *esd);
   void     ComputePID(Double_t energy, Double_t lambda0); // give the PID of a cluster
-  TArrayD  DistLambda0(Double_t energy, Int_t nature); // compute lambda0 distributions
+
+  void     InitParameters();
+  //void     InitParameters(Bool_t reconstructor);
+  void     SetLowFluxParam();
+  void     SetHighFluxParam();
+
+  TArrayD  DistLambda0(const Double_t energy, const Int_t nature) ; // compute lambda0 distributions
   
+  Double_t DistEnergy(const Double_t energy, const Int_t nature) ;
+
   Double_t GetPID(Int_t idx) const {if (idx>=0&&idx<3) return fPID[idx]; else return 0.;}
   Double_t GetPIDFinal(Int_t idx) const {if (idx>=0&&idx<AliPID::kSPECIESN) return fPIDFinal[idx]; else return 0.;}
   Double_t GetPIDWeight(Int_t idx) const {if (idx>=0&&idx<3) return fPIDWeight[idx]; else return 0.;}
   
-  void     SetPID(Double_t val, Int_t idx) {if (idx>=0&&idx<3) fPID[idx] = val;}
-  void     SetPIDFinal(Double_t val, Int_t idx) {if (idx>=0&&idx<AliPID::kSPECIESN) fPIDFinal[idx] = val;}
-  void     SetPIDWeight(Double_t val, Int_t idx) {if (idx>=0&&idx<3) fPIDWeight[idx] = val;}
-  void     SetPrintInfo(Bool_t yesno) {fPrintInfo = yesno;}
-   void     SetReconstructor(Bool_t yesno) {fReconstructor = yesno;}
+  void    SetPID(Double_t val, Int_t idx) {if (idx>=0&&idx<3) fPID[idx] = val;}
+  void    SetPIDFinal(Double_t val, Int_t idx) {if (idx>=0&&idx<AliPID::kSPECIESN) fPIDFinal[idx] = val;}
+  void    SetPIDWeight(Double_t val, Int_t idx) {if (idx>=0&&idx<3) fPIDWeight[idx] = val;}
+  void    SetPrintInfo(Bool_t yesno) {fPrintInfo = yesno;}
+  void    SetReconstructor(Bool_t yesno) {fReconstructor = yesno;}
+       
  private:
   
-  Double_t Polynomial(Double_t x, Double_t *params);
-  
+  Double_t Polynomial(const Double_t x, const Double_t *params) const ;
+  Double_t Polynomialinv(const Double_t x, const Double_t *params) const ;
+  Double_t PolynomialMixed1(const Double_t x, const Double_t *params) const ;
+  Double_t PolynomialMixed2(const Double_t x, const Double_t *params) const ;
+  Double_t Polynomial0(const Double_t *params) const ;
+  Double_t PowerExp(const Double_t x, const Double_t *params) const ;
+       
   Bool_t   fPrintInfo;          // flag to decide if details about PID must be printed
   
-  Double_t fGamma[6][6];        // Parameter to Compute PID
-  Double_t fHadron[6][6];                // Parameter to Compute PID
-  Double_t fPiZero5to10[6][6];  // Parameter to Compute PID
-  Double_t fPiZero10to60[6][6]; // Parameter to Compute PID
-  
+  Double_t fGamma[6][6];            // Parameter to Compute PID for photons
+  Double_t fGamma1to10[6][6];       // Parameter to Compute PID not used
+  Double_t fHadron[6][6];              // Parameter to Compute PID for hadrons, 1 to 10 GeV
+  Double_t fHadron1to10[6][6];     // Parameter to Compute PID for hadrons, 1 to 10 GeV
+  Double_t fPiZero[6][6];           // Parameter to Compute PID for pi0
+  Double_t fHadronEnergyProb[6];       // Parameter to Compute PID for energy ponderation for hadrons           
+  Double_t fPiZeroEnergyProb[6];       // Parameter to Compute PID for energy ponderation for Pi0       
+  Double_t fGammaEnergyProb[6];        // Parameter to Compute PID for energy ponderation for gamma     
+   
   Float_t fPID[3];
   
-  Float_t fPIDFinal[AliPID::kSPECIESN+1];  // final PID format
-  Float_t fPIDWeight[3];                 // order: gamma, pi0, hadrons,
-  Double_t fProbGamma;                 // probility to be a Gamma
-  Double_t fProbPiZero;                        // probility to be a PiO
-  Double_t fProbHadron;                        // probility to be a Hadron
-  Bool_t    fReconstructor;               //Fill esdcalocluster when called from EMCALReconstructor
+  Float_t fPIDFinal[AliPID::kSPECIESN+1]; // final PID format
+  Float_t fPIDWeight[3];                  // order: gamma, pi0, hadrons,
+  Double_t fProbGamma;                   // probility to be a Gamma
+  Double_t fProbPiZero;                          // probility to be a PiO
+  Double_t fProbHadron;                          // probility to be a Hadron
+  Double_t fWeightHadronEnergy;                  // Weight for a  a Hadron to have a given energy  (parametr from a flat distrib from 0 to 100)
+  Double_t fWeightGammaEnergy;           // Weight for a  Gamma to have a given energy  (for the moment =1.)
+  Double_t fWeightPiZeroEnergy;                  // Weight for a Pi0 Hadron to have a given energy (for the moment =1.)
+  Bool_t   fReconstructor;                // Fill esdcalocluster when called from EMCALReconstructor
   
-  ClassDef(AliEMCALPID, 0)
+  ClassDef(AliEMCALPID, 4)
 };
 
 #endif // ALIEMCALPID_H
index c7b2560..bdecde8 100644 (file)
@@ -1,4 +1,4 @@
- /**************************************************************************
+/**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
  * Author: The ALICE Off-line Project.                                    *
  **************************************************************************/
 
 /* $Id$ */
-// --- AliRoot header files ---
-#include "TObjArray.h"
-#include "AliCDBManager.h"
-#include "AliCDBEntry.h"
-#include "AliEMCALRecParam.h"
-#include "AliLog.h"
-
-ClassImp(AliEMCALRecParam)
-
-TObjArray* AliEMCALRecParam::fgkMaps =0; //ALTRO mappings 
 
 //-----------------------------------------------------------------------------
 // Container of EMCAL reconstruction parameters
 // The purpose of this object is to store it to OCDB
-// and retrieve it in AliEMCALClusterizerv1
+// and retrieve it in the corresponding reconstruction class:
+// AliEMCALClusterizer, AliEMCALPID, AliEMCALTracker ...
+//
 // Author: Yuri Kharlov
 //-----------------------------------------------------------------------------
 
+// --- Root header files
+//#include "TObjArray.h"
+
+// --- AliRoot header files ---
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
+#include "AliEMCALRecParam.h"
+
+ClassImp(AliEMCALRecParam)
+  
+  TObjArray* AliEMCALRecParam::fgkMaps =0; //ALTRO mappings 
+
 AliEMCALRecParam::AliEMCALRecParam() :
   AliDetectorRecoParam(),
   fClusteringThreshold(0.5),
@@ -55,116 +59,169 @@ AliEMCALRecParam::AliEMCALRecParam() :
   fNPedSamples(5) //raw signal
 {
   // default reco values
-
-  //PID parameters (Guenole)
-
+  
+  // PID parameters for Pb Pb from Lambda0 distributions fitted by  
+  // a landau inverted + Gaussian for Gammas
+  // and a Landau +Gaussian for Pi0 and hadrons 
+  // New parametrisation for 
+  // lambda0^2 (=x): f(x) = normLandau*TMath::Landau(((1-mpvlandau)-x),mpvLandau,widthLandau)+normgaus*TMath::Gaus(x,meangaus,sigmagaus) for gammas
+  // lambda0^2 (=x): f(x) = normLandau*TMath::Landau(x,mpvLandau,widthLandau)+normgaus*TMath::Gaus(x,meangaus,sigmagaus) for pi0 & hadrons
+  
+  // See AliEMCALPid 
+  // (index i) refers to each parameters of the f(lambda0^2) 
+  // i=0: normGaus
+  // i=1: meanGaus
+  // i=2: sigmaGaus
+  // i=3: normLandau
+  // i=4: mpvLandau
+  // i=5: sigmaLanda
+  // (index j) refers to the polynomial parameters of the fit of each parameter vs energy
+  // Pb Pb
+  
   // as a first step, all array elements are initialized to 0.0
   Int_t i, j;
   for (i = 0; i < 6; i++) {
-    for (j = 0; j < 6; j++) {
-      fGamma[i][j] = fHadron[i][j] = fPiZero5to10[i][j] = fPiZero10to60[i][j] = 0.;
+    for (j = 0; j < 6; j++) {      
+      fGamma[i][j] =  fPiZero[i][j] = fHadron[i][j] = 0.; 
+      fGamma1to10[i][j] =  fHadron1to10[i][j]= 0.;
     }
+    fGammaEnergyProb[i]=0.; // not yet implemented
+    fHadronEnergyProb[i]=0.; 
+    fPiZeroEnergyProb[i]=0.; // not yet implemented
+    
+    
   }
-
-  // then, only the ones which must be not zero are initialized
-  // while the others will remain to the value 0.0
-
-  fGamma[0][0] =  0.038022;
-  fGamma[0][1] = -0.0001883;
-  fGamma[0][2] =  5.449e-06;
-
-  fGamma[1][0] =  0.207313;
-  fGamma[1][1] = -0.000978;
-  fGamma[1][2] =  0.00001634;
-
-  fGamma[2][0] =  0.043364;
-  fGamma[2][1] = -0.0002048;
-  fGamma[2][2] =  8.661e-06;
-  fGamma[2][3] = -1.353e-07;
-
-  fGamma[3][0] =  0.265004;
-  fGamma[3][1] =  0.061298;
-  fGamma[3][2] = -0.003203;
-  fGamma[3][3] =  4.73e-05;
-
-  fGamma[4][0] =  0.243579;
-  fGamma[4][1] = -1.614e-05;
-
-  fGamma[5][0] =  0.002942;
-  fGamma[5][1] = -3.976e-05;
-
-  fHadron[0][0] =  0.011945 / 3.;
-  fHadron[0][1] =  0.000386 / 3.;
-  fHadron[0][2] = -0.000014 / 3.;
-  fHadron[0][3] =  1.336e-07 / 3.;
-
-  fHadron[1][0] =  0.496544;
-  fHadron[1][1] = -0.003226;
-  fHadron[1][2] =  0.00001678;
-
-  fHadron[2][0] =  0.144838;
-  fHadron[2][1] = -0.002954;
-  fHadron[2][2] =  0.00008754;
-  fHadron[2][3] = -7.587e-07;
-
-  fHadron[3][0] =  1.264461 / 7.;
-  fHadron[3][1] =  0.002097 / 7.;
-
-  fHadron[4][0] =  0.261950;
-  fHadron[4][1] = -0.001078;
-  fHadron[4][2] =  0.00003237;
-  fHadron[4][3] = -3.241e-07;
-  fHadron[4][4] =  0.;
-  fHadron[4][5] =  0.;
-  fHadron[5][0] =  0.010317;
-  fHadron[5][1] =  0.;
-  fHadron[5][2] =  0.;
-  fHadron[5][3] =  0.;
-  fHadron[5][4] =  0.;
-  fHadron[5][5] =  0.;
-
-  fPiZero5to10[0][0] = 0.009138;
-  fPiZero5to10[0][1] = 0.0006377;
-
-  fPiZero5to10[1][0] = 0.08;
-
-  fPiZero5to10[2][0] = -0.061119;
-  fPiZero5to10[2][1] =  0.019013;
-
-  fPiZero5to10[3][0] =  0.2;
-
-  fPiZero5to10[4][0] =  0.252044;
-  fPiZero5to10[4][1] = -0.002315;
-
-  fPiZero5to10[5][0] =  0.002942;
-  fPiZero5to10[5][1] = -3.976e-05;
+  // Pb Pb 
   
-  fPiZero10to60[0][0] =  0.009138;
-  fPiZero10to60[0][1] =  0.0006377;
-
-  fPiZero10to60[1][0] =  1.272837;
-  fPiZero10to60[1][1] = -0.069708;
-  fPiZero10to60[1][2] =  0.001568;
-  fPiZero10to60[1][3] = -1.162e-05;
-
-  fPiZero10to60[2][0] =  0.139703;
-  fPiZero10to60[2][1] =  0.003687;
-  fPiZero10to60[2][2] = -0.000568;
-  fPiZero10to60[2][3] =  1.498e-05;
-  fPiZero10to60[2][4] = -1.174e-07;
-
-  fPiZero10to60[3][0] = -0.826367;
-  fPiZero10to60[3][1] =  0.096951;
-  fPiZero10to60[3][2] = -0.002215;
-  fPiZero10to60[3][3] =  2.523e-05;
+  fGamma[0][0] = -7.656908e-01; 
+  fGamma[0][1] =  2.352536e-01; 
+  fGamma[0][2] =  1.555996e-02;
+  fGamma[0][3] =  2.243525e-04;
+  fGamma[0][4] = -2.560087e-06;
+  
+  fGamma[1][0] =  6.500216e+00;
+  fGamma[1][1] = -2.564958e-01;
+  fGamma[1][2] =  1.967894e-01;
+  fGamma[1][3] = -3.982273e-04;
+  fGamma[1][4] =  2.797737e-06;
+  
+  fGamma[2][0] =  2.416489e+00;
+  fGamma[2][1] = -1.601258e-01;
+  fGamma[2][2] =  3.126839e-02;
+  fGamma[2][3] =  3.387532e-04;
+  fGamma[2][4] = -4.089145e-06;
+  
+  
+  fGamma[3][0] =  0.;
+  fGamma[3][1] = -2.696008e+00;
+  fGamma[3][2] =  6.920305e-01;
+  fGamma[3][3] = -2.281122e-03;
+  fGamma[3][4] =  0.;
+  
+  fGamma[4][0] =  2.281564e-01;
+  fGamma[4][1] = -7.575040e-02;
+  fGamma[4][2] =  3.813423e-01;
+  fGamma[4][3] = -1.243854e-04;
+  fGamma[4][4] =  1.232045e-06;
+  
+  fGamma[5][0] = -3.290107e-01;
+  fGamma[5][1] =  3.707545e-02;
+  fGamma[5][2] =  2.917397e-03;
+  fGamma[5][3] =  4.695306e-05;
+  fGamma[5][4] = -3.572981e-07;
+  
+  
+  fHadron[0][0] =   1.519112e-01;
+  fHadron[0][1] = -8.267603e-02;
+  fHadron[0][2] =  1.914574e-02;
+  fHadron[0][3] = -2.677921e-04;
+  fHadron[0][4] =  5.447939e-06;
+  
+  
+  fHadron[1][0] = 0.;
+  fHadron[1][1] = -7.549870e-02; 
+  fHadron[1][2] = 3.930087e-01;
+  fHadron[1][3] = -2.368500e-03; 
+  fHadron[1][4] = 0.;
+  
+  
+  fHadron[2][0] = 0.;
+  fHadron[2][1] =  -2.463152e-02;
+  fHadron[2][2] = 1.349257e-01;
+  fHadron[2][3] = -1.089440e-03;
+  fHadron[2][4] = 0.;
+  
+  
+  
+  fHadron[3][0] = 0.;
+  fHadron[3][1] = 5.101560e-01;
+  fHadron[3][2] = 1.458679e-01;
+  fHadron[3][3] = 4.903068e-04;
+  fHadron[3][4] = 0.;
+  
+  fHadron[4][0] = 0.;
+  fHadron[4][1] = -6.693943e-03; 
+  fHadron[4][2] =  2.444753e-01;
+  fHadron[4][3] = -5.553749e-05;
+  fHadron[4][4] = 0.;
+  
+  fHadron[5][0] = -4.414030e-01;
+  fHadron[5][1] = 2.292277e-01;
+  fHadron[5][2] = -2.433737e-02;
+  fHadron[5][3] =  1.758422e-03;
+  fHadron[5][4] = -3.001493e-05;
+  
+  
+  fPiZero[0][0] =  5.072157e-01;
+  fPiZero[0][1] = -5.352747e-01;
+  fPiZero[0][2] =  8.499259e-02;
+  fPiZero[0][3] = -3.687401e-03;
+  fPiZero[0][4] =  5.482280e-05;
+  
+  
+  fPiZero[1][0] =  4.590137e+02; 
+  fPiZero[1][1] = -7.079341e+01;
+  fPiZero[1][2] =  4.990735e+00;
+  fPiZero[1][3] = -1.241302e-01;
+  fPiZero[1][4] =  1.065772e-03;
+  
+  
+  fPiZero[2][0] =  1.376415e+02;
+  fPiZero[2][1] = -3.031577e+01;
+  fPiZero[2][2] =  2.474338e+00;
+  fPiZero[2][3] = -6.903410e-02;
+  fPiZero[2][4] =  6.244089e-04;
+  
+  fPiZero[3][0] = 0.;
+  fPiZero[3][1] =  1.145983e+00;
+  fPiZero[3][2] = -2.476052e-01;
+  fPiZero[3][3] =  1.367373e-02;
+  fPiZero[3][4] = 0.;
+  
+  fPiZero[4][0] = -2.097586e+02;
+  fPiZero[4][1] =  6.300800e+01;
+  fPiZero[4][2] = -4.038906e+00;
+  fPiZero[4][3] =  1.088543e-01;
+  fPiZero[4][4] = -9.362485e-04;
+  
+  fPiZero[5][0] = -1.671477e+01; 
+  fPiZero[5][1] =  2.995415e+00;
+  fPiZero[5][2] = -6.040360e-02;
+  fPiZero[5][3] = -6.137459e-04;
+  fPiZero[5][4] =  1.847328e-05;
+  
+  // High flux ones pp 
+  
+  fHadronEnergyProb[0]= 0.;
+  fHadronEnergyProb[1]= 0.;
+  fHadronEnergyProb[2]=  6.188452e-02;
+  fHadronEnergyProb[3]=  2.030230e+00;
+  fHadronEnergyProb[4]= -6.402242e-02;
+  
+  
+}
 
-  fPiZero10to60[4][0] =  0.249890;
-  fPiZero10to60[4][1] = -0.000063;
 
-  fPiZero10to60[5][0] =  0.002942;
-  fPiZero10to60[5][1] = -3.976e-05;
-
-}
 
 //-----------------------------------------------------------------------------
 AliEMCALRecParam::AliEMCALRecParam(const AliEMCALRecParam& rp) :
@@ -180,7 +237,7 @@ AliEMCALRecParam::AliEMCALRecParam(const AliEMCALRecParam& rp) :
   fTrkCutR(rp.fTrkCutR),
   fTrkCutAlphaMin(rp.fTrkCutAlphaMin), 
   fTrkCutAlphaMax(rp.fTrkCutAlphaMax), 
-  fTrkCutAngle(rp.fTrkCutAngle),
+  fTrkCutAngle(rp.fTrkCutAngle), 
   fTrkCutNITS(rp.fTrkCutNITS),
   fTrkCutNTPC(rp.fTrkCutNTPC), // track matching
   fHighLowGainFactor(rp.fHighLowGainFactor), 
@@ -190,25 +247,30 @@ AliEMCALRecParam::AliEMCALRecParam(const AliEMCALRecParam& rp) :
   fNPedSamples(rp.fNPedSamples) //raw signal
 {
   //copy constructor
-
+  
   //PID values
   Int_t i, j;
   for (i = 0; i < 6; i++) {
     for (j = 0; j < 6; j++) {
       fGamma[i][j] = rp.fGamma[i][j];
+      fGamma1to10[i][j] = rp.fGamma1to10[i][j];
       fHadron[i][j] = rp.fHadron[i][j];
-      fPiZero5to10[i][j] = rp.fPiZero5to10[i][j];
-      fPiZero10to60[i][j] = rp.fPiZero10to60[i][j];
+      fHadron1to10[i][j] = rp.fHadron1to10[i][j];
+      fPiZero[i][j] = rp.fPiZero[i][j];
     }
+    fGammaEnergyProb[i] = rp.fGammaEnergyProb[i];
+    fPiZeroEnergyProb[i] = rp.fPiZeroEnergyProb[i];
+    fHadronEnergyProb[i] = rp.fHadronEnergyProb[i];
+    
   }
-
+  
 }
 
 //-----------------------------------------------------------------------------
 AliEMCALRecParam& AliEMCALRecParam::operator = (const AliEMCALRecParam& rp)
 {
   //assignment operator
-
+  
   if(this != &rp) {
     fClusteringThreshold = rp.fClusteringThreshold;
     fW0 = rp.fW0;
@@ -221,7 +283,7 @@ AliEMCALRecParam& AliEMCALRecParam::operator = (const AliEMCALRecParam& rp)
     fTrkCutR = rp.fTrkCutR;
     fTrkCutAlphaMin = rp.fTrkCutAlphaMin;
     fTrkCutAlphaMax = rp.fTrkCutAlphaMax;
-    fTrkCutAngle = rp.fTrkCutAngle;
+    fTrkCutAngle = rp.fTrkCutAngle; 
     fTrkCutNITS = rp.fTrkCutNITS;
     fTrkCutNTPC = rp.fTrkCutNTPC; //track matching
     fHighLowGainFactor = rp.fHighLowGainFactor; 
@@ -229,22 +291,26 @@ AliEMCALRecParam& AliEMCALRecParam::operator = (const AliEMCALRecParam& rp)
     fTau = rp.fTau;
     fNoiseThreshold = rp.fNoiseThreshold;
     fNPedSamples = rp.fNPedSamples; //raw signal
-
+    
     //PID values
     Int_t i, j;
     for (i = 0; i < 6; i++) {
       for (j = 0; j < 6; j++) {
        fGamma[i][j] = rp.fGamma[i][j];
+       fGamma1to10[i][j] = rp.fGamma1to10[i][j];
        fHadron[i][j] = rp.fHadron[i][j];
-       fPiZero5to10[i][j] = rp.fPiZero5to10[i][j];
-       fPiZero10to60[i][j] = rp.fPiZero10to60[i][j];
+       fHadron1to10[i][j] = rp.fHadron1to10[i][j];
+       fPiZero[i][j] = rp.fPiZero[i][j];
       }
+      fGammaEnergyProb[i] = rp.fGammaEnergyProb[i];
+      fPiZeroEnergyProb[i] = rp.fPiZeroEnergyProb[i];
+      fHadronEnergyProb[i] = rp.fHadronEnergyProb[i];
     }
-
+    
   }    
   
   return *this;
-
+  
 }
 
 //-----------------------------------------------------------------------------
@@ -255,39 +321,40 @@ AliEMCALRecParam* AliEMCALRecParam::GetDefaultParameters()
   params->SetName("Default - Pb+Pb");
   params->SetTitle("Default - Pb+Pb");
   return params;
-
+  
 }
 
 //-----------------------------------------------------------------------------
 AliEMCALRecParam* AliEMCALRecParam::GetCalibParam()
 {
-       //parameters for the reconstruction of calibration runs
-       AliEMCALRecParam* params = new AliEMCALRecParam();
-       params->SetClusteringThreshold(0.2); // 200 MeV                                             
-       params->SetMinECut(0.01);  //10 MeV       
-       params->SetName("Calibration - LED");
-       params->SetTitle("Calibration - LED");
-       params->SetEventSpecie(AliRecoParam::kCalib);
-       
-       return params;
-       
+  //parameters for the reconstruction of calibration runs
+  AliEMCALRecParam* params = GetLowFluxParam();
+  //params->SetClusteringThreshold(0.2); // 200 MeV                                             
+  //params->SetMinECut(0.01);  //10 MeV       
+  params->SetName("Calibration - LED");
+  params->SetTitle("Calibration - LED");
+  params->SetEventSpecie(AliRecoParam::kCalib);
+  
+  return params;
+  
 }
 
 //-----------------------------------------------------------------------------
 AliEMCALRecParam* AliEMCALRecParam::GetCosmicParam()
 {
-       //parameters for the reconstruction of cosmic runs
-       AliEMCALRecParam* params = new AliEMCALRecParam();
-       params->SetClusteringThreshold(0.2); // 200 MeV                                             
-       params->SetMinECut(0.01);  //10 MeV       
-       params->SetName("Cosmic");
-       params->SetTitle("Cosmic");
-       params->SetEventSpecie(AliRecoParam::kCosmic);
-       
-       return params;
-       
+  //parameters for the reconstruction of cosmic runs
+  AliEMCALRecParam* params = GetLowFluxParam();
+  //params->SetClusteringThreshold(0.2); // 200 MeV                                             
+  //params->SetMinECut(0.01);  //10 MeV       
+  params->SetName("Cosmic");
+  params->SetTitle("Cosmic");
+  params->SetEventSpecie(AliRecoParam::kCosmic);
+  
+  return params;
+  
 }
 
+
 //-----------------------------------------------------------------------------
 AliEMCALRecParam* AliEMCALRecParam::GetLowFluxParam()
 {
@@ -298,7 +365,154 @@ AliEMCALRecParam* AliEMCALRecParam::GetLowFluxParam()
   params->SetName("Low Flux - p+p");
   params->SetTitle("Low Flux - p+p");
   params->SetEventSpecie(AliRecoParam::kLowMult);
-
+  
+  
+  //PID parameters for pp  implemented 
+  // as a first step, all array elements are initialized to 0.0
+  Int_t i, j;
+  for (i = 0; i < 6; i++) {
+    for (j = 0; j < 6; j++) {
+      params->SetGamma(i,j,0.);
+      params->SetGamma1to10(i,j,0.);
+      params->SetHadron(i,j,0.);
+      params->SetHadron1to10(i,j,0.);
+      params->SetPiZero(i,j,0.);
+      
+    }
+    params->SetGammaEnergyProb(i,0.); // not yet implemented
+    params->SetHadronEnergyProb(i,0.);
+    params->SetPiZeroEnergyProb(i,0.); // not yet implemented
+  }
+  
+  
+  params->SetGamma(0,0,-7.656908e-01);
+  params->SetGamma(0,1,2.352536e-01);
+  params->SetGamma(0,2,1.555996e-02);
+  params->SetGamma(0,3,2.243525e-04);
+  params->SetGamma(0,4,-2.560087e-06);
+
+  params->SetGamma(1,0,6.500216e+00);
+  params->SetGamma(1,1,-2.564958e-01);
+  params->SetGamma(1,2,1.967894e-01);
+  params->SetGamma(1,3,-3.982273e-04);
+  params->SetGamma(1,4,2.797737e-06);
+  
+  params->SetGamma(2,0,2.416489e+00);
+  params->SetGamma(2,1,-1.601258e-01);
+  params->SetGamma(2,2,3.126839e-02);
+  params->SetGamma(2,3,3.387532e-04);
+  params->SetGamma(2,4,-4.089145e-06);
+  params->SetGamma(3,0,0.);
+  params->SetGamma(3,1,-2.696008e+00); 
+  params->SetGamma(3,2, 6.920305e-01);
+  params->SetGamma(3,3,-2.281122e-03);
+  params->SetGamma(3,4,0.);
+
+  params->SetGamma(4,0,2.281564e-01); 
+  params->SetGamma(4,1,-7.575040e-02);
+  params->SetGamma(4,2,3.813423e-01);
+  params->SetGamma(4,3,-1.243854e-04);
+  params->SetGamma(4,4,1.232045e-06);
+
+  params->SetGamma(5,0,-3.290107e-01);
+  params->SetGamma(5,1,3.707545e-02);
+  params->SetGamma(5,2,2.917397e-03);
+  params->SetGamma(5,3,4.695306e-05);
+  params->SetGamma(5,4,-3.572981e-07);
+
+  params->SetHadron(0,0,9.482243e-01); 
+  params->SetHadron(0,1,-2.780896e-01);
+  params->SetHadron(0,2, 2.223507e-02);
+  params->SetHadron(0,3,7.294263e-04);
+  params->SetHadron(0,4,-5.665872e-06); 
+
+  params->SetHadron(1,0,0.);
+  params->SetHadron(1,1,0.);
+  params->SetHadron(1,2,2.483298e-01);
+  params->SetHadron(1,3,0.);
+  params->SetHadron(1,4,0.);
+
+  params->SetHadron(2,0,-5.601199e+00);
+  params->SetHadron(2,1,2.097382e+00);
+  params->SetHadron(2,2,-2.307965e-01);
+  params->SetHadron(2,3,9.206871e-03);
+  params->SetHadron(2,4,-8.887548e-05);
+
+  params->SetHadron(3,0,6.543101e+00);
+  params->SetHadron(3,1,-2.305203e+00);
+  params->SetHadron(3,2,2.761673e-01);
+  params->SetHadron(3,3,-5.465855e-03);
+  params->SetHadron(3,4,2.784329e-05);
+
+  params->SetHadron(4,0,-2.443530e+01);
+  params->SetHadron(4,1,8.902578e+00);
+  params->SetHadron(4,2,-5.265901e-01);
+  params->SetHadron(4,3,2.549111e-02);
+  params->SetHadron(4,4,-2.196801e-04);
+
+  params->SetHadron(5,0,2.102007e-01);
+  params->SetHadron(5,1,-3.844418e-02);
+  params->SetHadron(5,2,1.234682e-01);
+  params->SetHadron(5,3,-3.866733e-03);
+  params->SetHadron(5,4,3.362719e-05);
+
+  params->SetPiZero(0,0,5.07215e-01);
+  params->SetPiZero(0,1,-5.35274e-01);
+  params->SetPiZero(0,2,8.49925e-02);
+  params->SetPiZero(0,3,-3.68740e-03);
+  params->SetPiZero(0,4,5.48228e-05);
+
+  params->SetPiZero(1,0,4.590137e+02);
+  params->SetPiZero(1,1,-7.079341e+01);
+  params->SetPiZero(1,2,4.990735e+00);
+  params->SetPiZero(1,3,-1.241302e-01);
+  params->SetPiZero(1,4,1.065772e-03);
+
+  params->SetPiZero(2,0,1.376415e+02); 
+  params->SetPiZero(2,1,-3.031577e+01);
+  params->SetPiZero(2,2,2.474338e+00);
+  params->SetPiZero(2,3,-6.903410e-02); 
+  params->SetPiZero(2,4,6.244089e-04);
+
+  params->SetPiZero(3,0,0.);
+  params->SetPiZero(3,1,1.145983e+00);
+  params->SetPiZero(3,2,-2.476052e-01);
+  params->SetPiZero(3,3,1.367373e-02);
+  params->SetPiZero(3,4,0.);
+
+  params->SetPiZero(4,0,-2.09758e+02);
+  params->SetPiZero(4,1,6.30080e+01);
+  params->SetPiZero(4,2,-4.03890e+00);
+  params->SetPiZero(4,3,1.08854e-01);
+  params->SetPiZero(4,4,-9.36248e-04);
+
+  params->SetPiZero(5,0,-1.671477e+01);
+  params->SetPiZero(5,1,2.995415e+00);
+  params->SetPiZero(5,2,-6.040360e-02);
+  params->SetPiZero(5,3,-6.137459e-04);
+  params->SetPiZero(5,4,1.847328e-05);
+
+//     params->SetHadronEnergyProb(0,0.);
+//     params->SetHadronEnergyProb(1,0.);
+//     params->SetHadronEnergyProb(2,1.);
+//     params->SetHadronEnergyProb(3,0.);
+//     params->SetHadronEnergyProb(4,0.);
+  params->SetHadronEnergyProb(0, 4.767543e-02);
+  params->SetHadronEnergyProb(1,-1.537523e+00);
+  params->SetHadronEnergyProb(2,2.956727e-01);
+  params->SetHadronEnergyProb(3,-3.051022e+01);
+  params->SetHadronEnergyProb(4,-6.036931e-02);
+
+//   Int_t ii= 0;
+//   Int_t jj= 3;
+//     AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
+//                     ii,jj, params->GetGamma(ii,jj), params->GetPiZero(ii,jj), params->GetHadron(ii,jj)));
+//     cout << " Low Flux Parameters fGamma [2][2] = " << params->GetGamma(2,2) << endl;
+//     cout << " Low Flux Parameters fHadron [2][2] = " << params->GetHadron(2,2) << endl;
+   
   return params;
 
 }
@@ -314,9 +528,9 @@ AliEMCALRecParam* AliEMCALRecParam::GetHighFluxParam()
   params->SetName("High Flux - Pb+Pb");
   params->SetTitle("High Flux - Pb+Pb");
   params->SetEventSpecie(AliRecoParam::kHighMult);
-
+  
   return params;
-
+  
 }
 
 //-----------------------------------------------------------------------------
@@ -325,9 +539,9 @@ void AliEMCALRecParam::Print(Option_t *) const
   // Print reconstruction parameters to stdout
   AliInfo(Form("Clusterization parameters :\n fClusteringThreshold=%.3f,\n fW0=%.3f,\n fMinECut=%.3f,\n fUnfold=%d,\n fLocMaxCut=%.3f \n",
               fClusteringThreshold,fW0,fMinECut,fUnfold,fLocMaxCut));
-
-  AliInfo(Form("Track-matching cuts :\n x %f, y %f, z %f, R %f \n alphaMin %f, alphaMax %f, Angle %f, NITS %f, NTPC %f\n", fTrkCutX, fTrkCutY, fTrkCutZ, fTrkCutR,fTrkCutAlphaMin,fTrkCutAlphaMax, fTrkCutAngle,fTrkCutNITS,fTrkCutNTPC));
-
+  
+  AliInfo(Form("Track-matching cuts :\n x %f, y %f, z %f, R %f \n alphaMin %f, alphaMax %f, Angle %f, NITS %f, NTPC %\n", fTrkCutX, fTrkCutY, fTrkCutZ, fTrkCutR,fTrkCutAlphaMin,fTrkCutAlphaMax, fTrkCutAngle,fTrkCutNITS,fTrkCutNTPC));
+  
   AliInfo(Form("PID parameters, Gamma :\n"));
   for(Int_t i = 0; i < 6; i++){
     for(Int_t j = 0; j < 6; j++){
@@ -335,9 +549,8 @@ void AliEMCALRecParam::Print(Option_t *) const
     }
     printf("\n");
   }
-
-  printf("\n");
-
+  
+  
   AliInfo(Form("PID parameters, Hadron :\n"));
   for(Int_t i = 0; i < 6; i++){
     for(Int_t j = 0; j < 6; j++){
@@ -345,32 +558,23 @@ void AliEMCALRecParam::Print(Option_t *) const
     }
     printf("\n");
   }
-
-  printf("\n");
-
-  AliInfo(Form("PID parameters, Pi0zero5to10 :\n"));
-  for(Int_t i = 0; i < 6; i++){
-    for(Int_t j = 0; j < 6; j++){
-      printf(" %f, ", fPiZero5to10[i][j]);
-    }
-    printf("\n");
-  }
-
+  
   printf("\n");
-
-  AliInfo(Form("PID parameters, Pi0zero10to60 :\n"));
+  
+  AliInfo(Form("PID parameters, Pi0zero :\n"));
   for(Int_t i = 0; i < 6; i++){
     for(Int_t j = 0; j < 6; j++){
-      printf(" %f, ", fPiZero10to60[i][j]);
+      printf(" %f, ", fPiZero[i][j]);
     }
     printf("\n");
   }
-
+  
   printf("\n");
-
+  
+  
   AliInfo(Form("Raw signal parameters: \n gain factor=%f, order=%d, tau=%f, noise threshold=%d, nped samples=%d \n",
               fHighLowGainFactor,fOrderParameter,fTau,fNoiseThreshold,fNPedSamples));
-
+  
 }
 
 //-----------------------------------------------------------------------------
@@ -378,20 +582,20 @@ const TObjArray* AliEMCALRecParam::GetMappings()
 {
   //Returns array of AliAltroMappings for RCU0..RCUX.
   //If not found, read it from OCDB.                                           
-
+  
   //Quick check as follows:                                                   
-  //  root [0] AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT/OCDB"
+  //  root [0] AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT"
   //  root [1] AliCDBManager::Instance()->SetRun(1);                             
   //  root [2] TObjArray* maps = AliEMCALRecParam::GetMappings();                
   //  root [3] maps->Print();                                                    
-
+  
   if(fgkMaps) return fgkMaps;
-
+  
   AliCDBEntry* entry = AliCDBManager::Instance()->Get("EMCAL/Calib/Mapping");
   if(entry)
     fgkMaps = (TObjArray*)entry->GetObject();
-
+  
   return fgkMaps;
-
+  
 }
 
index 4c5783b..69bb77d 100644 (file)
 
 class AliEMCALRecParam : public AliDetectorRecoParam
 {
-public:
+ public:
   
   AliEMCALRecParam() ;
   AliEMCALRecParam(const AliEMCALRecParam& recParam);
   AliEMCALRecParam& operator = (const AliEMCALRecParam& recParam);
   virtual ~AliEMCALRecParam() {}
+  
   //Clustering (Unfolding : Cynthia)
-  Float_t GetClusteringThreshold() const     {return fClusteringThreshold;}
-  Float_t GetW0                 () const     {return fW0                 ;}
-  Float_t GetMinECut            () const     {return fMinECut            ;}
-  Float_t GetLocMaxCut          () const     {return fLocMaxCut            ;}
-  Bool_t  GetUnfold             () const     {return fUnfold            ;}
-  void SetClusteringThreshold(Float_t thrsh)   {fClusteringThreshold = thrsh;}
-  void SetW0                 (Float_t w0)      {fW0 = w0                    ;}
-  void SetMinECut            (Float_t minEcut) {fMinECut = minEcut          ;}
-  void SetLocMaxCut          (Float_t locMaxCut) {fLocMaxCut = locMaxCut    ;}
-  void SetUnfold             (Bool_t unfold)     {fUnfold = unfold          ; if(fUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!");}
-
+  Float_t GetClusteringThreshold() const     {return fClusteringThreshold ;}
+  Float_t GetW0                 () const     {return fW0                  ;}
+  Float_t GetMinECut            () const     {return fMinECut             ;}
+  Float_t GetLocMaxCut          () const     {return fLocMaxCut           ;}
+  Bool_t  GetUnfold             () const     {return fUnfold              ;}
+  void SetClusteringThreshold(Float_t thrsh)     {fClusteringThreshold = thrsh;}
+  void SetW0                 (Float_t w0)        {fW0 = w0                ;}
+  void SetMinECut            (Float_t minEcut)   {fMinECut = minEcut      ;}
+  void SetLocMaxCut          (Float_t locMaxCut) {fLocMaxCut = locMaxCut  ;}
+  void SetUnfold             (Bool_t unfold)     {fUnfold = unfold ; if(fUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!");}
+  
   //PID (Guenole)
-  Double_t GetGamma(Int_t i, Int_t j) const    {return fGamma[i][j];} 
-  Double_t GetHadron(Int_t i, Int_t j) const    {return fHadron[i][j];}
-  Double_t GetPiZero5to10(Int_t i, Int_t j) const    {return fPiZero5to10[i][j];}
-  Double_t GetPiZero10to60(Int_t i, Int_t j) const    {return fPiZero10to60[i][j];}
-
-  void SetGamma(Int_t i, Int_t j,Double_t param )   {fGamma[i][j]=param;}
-  void SetHadron(Int_t i, Int_t j,Double_t param )   {fHadron[i][j]=param;}
-  void SetPiZero5to10(Int_t i, Int_t j,Double_t param)   {fPiZero5to10[i][j]=param;}
-  void SetPiZero10to60(Int_t i, Int_t j,Double_t param)   {fPiZero10to60[i][j]=param;}
-
+  Double_t GetGamma(Int_t i, Int_t j) const       {return fGamma[i][j];} 
+  Double_t GetGammaEnergyProb(Int_t i) const      {return fGammaEnergyProb[i];} 
+  Double_t GetGamma1to10(Int_t i, Int_t j) const  {return fGamma1to10[i][j];}   // not used
+  Double_t GetHadron(Int_t i, Int_t j) const      {return fHadron[i][j];}
+  Double_t GetHadron1to10(Int_t i, Int_t j) const {return fHadron1to10[i][j];}   // not used
+  Double_t GetHadronEnergyProb(Int_t i) const     {return fHadronEnergyProb[i];}
+  Double_t GetPiZero(Int_t i, Int_t j) const      {return fPiZero[i][j];}
+  Double_t GetPiZeroEnergyProb(Int_t i) const     {return fPiZeroEnergyProb[i];}
+  
+  void SetGamma(Int_t i, Int_t j,Double_t param )       {fGamma[i][j]=param;}
+  void SetGammaEnergyProb(Int_t i, Double_t param )     {fGammaEnergyProb[i]=param;}
+  void SetGamma1to10(Int_t i, Int_t j,Double_t param )  {fGamma1to10[i][j]=param;}
+  void SetHadron(Int_t i, Int_t j,Double_t param )      {fHadron[i][j]=param;}
+  void SetHadron1to10(Int_t i, Int_t j,Double_t param ) {fHadron1to10[i][j]=param;}
+  void SetHadronEnergyProb(Int_t i,Double_t param )     {fHadronEnergyProb[i]=param;}
+  void SetPiZero(Int_t i, Int_t j,Double_t param)       {fPiZero[i][j]=param;}
+  void SetPiZeroEnergyProb(Int_t i,Double_t param)      {fPiZeroEnergyProb[i]=param;}
+  
   //Track Matching (Alberto)
   /* track matching cut setters */
   void SetTrkCutX(Double_t value)        {fTrkCutX = value;}
@@ -61,8 +69,8 @@ public:
   void SetTrkCutAlphaMin(Double_t value) {fTrkCutAlphaMin = value;}
   void SetTrkCutAlphaMax(Double_t value) {fTrkCutAlphaMax = value;}
   void SetTrkCutAngle(Double_t value)    {fTrkCutAngle = value;}
-  void SetTrkCutNITS(Double_t value)        {fTrkCutNITS = value;}
-  void SetTrkCutNTPC(Double_t value)        {fTrkCutNTPC = value;}
+  void SetTrkCutNITS(Double_t value)     {fTrkCutNITS = value;}
+  void SetTrkCutNTPC(Double_t value)     {fTrkCutNTPC = value;}
   /* track matching cut getters */
   Double_t GetTrkCutX() const        {return fTrkCutX;}
   Double_t GetTrkCutY() const        {return fTrkCutY;}
@@ -71,9 +79,9 @@ public:
   Double_t GetTrkCutAlphaMin() const {return fTrkCutAlphaMin;}
   Double_t GetTrkCutAlphaMax() const {return fTrkCutAlphaMax;}
   Double_t GetTrkCutAngle() const    {return fTrkCutAngle;}
-  Double_t GetTrkCutNITS() const        {return fTrkCutNITS;}
-  Double_t GetTrkCutNTPC() const        {return fTrkCutNTPC;}
-
+  Double_t GetTrkCutNITS() const     {return fTrkCutNITS;}
+  Double_t GetTrkCutNTPC() const     {return fTrkCutNTPC;}
+  
   //Raw signal fitting (Jenn)
   /* raw signal setters */
   void SetHighLowGainFactor(Double_t value) {fHighLowGainFactor = value;}
@@ -87,31 +95,36 @@ public:
   Double_t GetTau()               const {return fTau;}
   Int_t    GetNoiseThreshold()    const {return fNoiseThreshold;}
   Int_t    GetNPedSamples()       const {return fNPedSamples;}
+  
   virtual void Print(Option_t * option="") const ;
-
+  
   static AliEMCALRecParam* GetDefaultParameters();
   static AliEMCALRecParam* GetLowFluxParam();
   static AliEMCALRecParam* GetHighFluxParam();
   static AliEMCALRecParam* GetCalibParam();
   static AliEMCALRecParam* GetCosmicParam();
-
+  
   static const  TObjArray* GetMappings();
-
-private:
+  
+ private:
   //Clustering
   Float_t fClusteringThreshold ; // minimum energy to seed a EC digit in a cluster
   Float_t fW0 ;                  // logarithmic weight for the cluster center of gravity calculation
   Float_t fMinECut;              // Minimum energy for a digit to be a member of a cluster
-  Bool_t fUnfold;               // flag to perform cluster unfolding
+  Bool_t fUnfold;                // flag to perform cluster unfolding
   Float_t fLocMaxCut;            // minimum energy difference to consider local maxima in a cluster
-
+  
   //PID (Guenole)
-  Double_t fGamma[6][6];        // Parameter to Compute PID      
-  Double_t fHadron[6][6];      // Parameter to Compute PID      
-  Double_t fPiZero5to10[6][6];  // Parameter to Compute PID     
-  Double_t fPiZero10to60[6][6]; // Parameter to Compute PID     
-
+  Double_t fGamma[6][6];         // Parameter to Compute PID for photons     
+  Double_t fGamma1to10[6][6];    // Parameter to Compute PID not used
+  Double_t fHadron[6][6];           // Parameter to Compute PID for hadrons     
+  Double_t fHadron1to10[6][6];          // Parameter to Compute PID for hadrons between 1 and 10 GeV    
+  Double_t fHadronEnergyProb[6]; // Parameter to Compute PID for energy ponderation for hadrons         
+  Double_t fPiZeroEnergyProb[6]; // Parameter to Compute PID for energy ponderation for Pi0     
+  Double_t fGammaEnergyProb[6];  // Parameter to Compute PID for energy ponderation for gamma           
+  Double_t fPiZero[6][6];        // Parameter to Compute PID for pi0    
+  
+  
   //Track-Matching (Alberto)
   Double_t  fTrkCutX;              // X-difference cut for track matching
   Double_t  fTrkCutY;              // Y-difference cut for track matching
@@ -120,21 +133,21 @@ private:
   Double_t  fTrkCutAlphaMin;       // cut on 'alpha' parameter for track matching (min)
   Double_t  fTrkCutAlphaMax;       // cut on 'alpha' parameter for track matching (min)
   Double_t  fTrkCutAngle;          // cut on relative angle between different track points for track matching
-  Double_t  fTrkCutNITS;              // Number of ITS hits for track matching
-  Double_t  fTrkCutNTPC;              // Number of TPC hits for track matching
+  Double_t  fTrkCutNITS;           // Number of ITS hits for track matching
+  Double_t  fTrkCutNTPC;           // Number of TPC hits for track matching
+  
   //Raw signal fitting parameters (Jenn)
   Double_t fHighLowGainFactor;     //gain factor to convert between high and low gain
   Int_t    fOrderParameter;        //order parameter for raw signal fit
   Double_t fTau;                   //decay constant for raw signal fit
   Int_t    fNoiseThreshold;        //threshold to consider signal or noise
   Int_t    fNPedSamples;           //number of time samples to use in pedestal calculation
-
+  
   static TObjArray* fgkMaps;       // ALTRO mappings for RCU0..RCUX
-
-  ClassDef(AliEMCALRecParam,6)   // Reconstruction parameters
-
-} ;
+  
+  ClassDef(AliEMCALRecParam,7)     // Reconstruction parameters
+    
+    } ;
 
 #endif //  ALIEMCALRECPARAM_H
 
index 3cf595b..d83d56e 100644 (file)
 
 void AliEMCALSetRecParamCDB(AliRecoParam::EventSpecie_t default = AliRecoParam::kDefault)
 {
-
+  
   // Create an object AliEMCALRecParam and store it to OCDB
-
+  
   //Activate CDB storage
   AliCDBManager* cdb = AliCDBManager::Instance();
   if(!cdb->IsDefaultStorageSet()) cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
-
-
+  
   // Create reconstruction parameter object and set parameter values
   TObjArray* recParamArray = new TObjArray();
-
+  
   {
     //default
-    AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetDefaultParameters();
-    
-    //Clusterization
-    recParamDB->SetClusteringThreshold(0.5);
-    recParamDB->SetW0(4.5);
-    recParamDB->SetMinECut(0.45);
-    recParamDB->SetUnfold(kFALSE);
-    recParamDB->SetLocMaxCut(0.03);
-    
-    //Track matching
-    recParamDB->SetTrkCutX(6.0);
-    recParamDB->SetTrkCutY(6.0);
-    recParamDB->SetTrkCutZ(6.0);
-    recParamDB->SetTrkCutR(10.0);
-    recParamDB->SetTrkCutAlphaMin(-50.0);
-    recParamDB->SetTrkCutAlphaMax( 50.0);
-    recParamDB->SetTrkCutNITS(3.0);
-    recParamDB->SetTrkCutNTPC(20.0);
-    recParamDB->SetTrkCutAngle(10000.0);      // i.e. exclude this for the moment
-    
-    //PID
-    
-    // as a first step, all array elements are initialized to 0.0
-    Int_t i, j;
-    for (i = 0; i < 6; i++) {
-      for (j = 0; j < 6; j++) {
-       
-       recParamDB->SetGamma(i,j,0.);
-       recParamDB->SetHadron(i,j,0.);
-       recParamDB->SetPiZero5to10(i,j, 0.);
-       recParamDB->SetPiZero10to60(i,j,0.);
-      }
-    } 
-    
-    recParamDB->SetGamma(0,0,  0.038022);
-    recParamDB->SetGamma(0,1, -0.0001883);
-    recParamDB->SetGamma(0,2,  5.449e-06);
-    
-    recParamDB->SetGamma(1,0,  0.207313);
-    recParamDB->SetGamma(1,1, -0.000978);
-    recParamDB->SetGamma(1,2,  0.00001634);
-    
-    recParamDB->SetGamma(2,0,  0.043364);
-    recParamDB->SetGamma(2,1, -0.0002048);
-    recParamDB->SetGamma(2,2,  8.661e-06);
-    recParamDB->SetGamma(2,3, -1.353e-07);
-    
-    recParamDB->SetGamma(3,0,  0.265004);
-    recParamDB->SetGamma(3,1,  0.061298);
-    recParamDB->SetGamma(3,2, -0.003203);
-    recParamDB->SetGamma(3,3,  4.73e-05);
-    
-    recParamDB->SetGamma(4,0,  0.243579);
-    recParamDB->SetGamma(4,1, -1.614e-05);
-    
-    recParamDB->SetGamma(5,0,  0.002942);
-    recParamDB->SetGamma(5,1, -3.976e-05);
-    
-    recParamDB->SetHadron(0,0,  0.011945 / 3.);
-    recParamDB->SetHadron(0,1,  0.000386 / 3.);
-    recParamDB->SetHadron(0,2, -0.000014 / 3.);
-    recParamDB->SetHadron(0,3,  1.336e-07 / 3.);
-    
-    recParamDB->SetHadron(1,0,  0.496544);
-    recParamDB->SetHadron(1,1, -0.003226);
-    recParamDB->SetHadron(1,2,  0.00001678);
-    
-    recParamDB->SetHadron(2,0,  0.144838);
-    recParamDB->SetHadron(2,1, -0.002954);
-    recParamDB->SetHadron(2,2,  0.00008754);
-    recParamDB->SetHadron(2,3, -7.587e-07);
-    
-    recParamDB->SetHadron(3,0,  1.264461 / 7.);
-    recParamDB->SetHadron(3,1,  0.002097 / 7.);
-    
-    recParamDB->SetHadron(4,0,  0.261950);
-    recParamDB->SetHadron(4,1, -0.001078);
-    recParamDB->SetHadron(4,2,  0.00003237);
-    recParamDB->SetHadron(4,3, -3.241e-07);
-    recParamDB->SetHadron(4,4,  0.);
-    recParamDB->SetHadron(4,5,  0.);
-    recParamDB->SetHadron(5,0,  0.010317);
-    recParamDB->SetHadron(5,1,  0.);
-    recParamDB->SetHadron(5,2,  0.);
-    recParamDB->SetHadron(5,3,  0.);
-    recParamDB->SetHadron(5,4,  0.);
-    recParamDB->SetHadron(5,5,  0.);
-    
-    recParamDB->SetPiZero5to10(0,0, 0.009138);
-    recParamDB->SetPiZero5to10(0,1, 0.0006377);
-    
-    recParamDB->SetPiZero5to10(1,0, 0.08);
-    
-    recParamDB->SetPiZero5to10(2,0, -0.061119);
-    recParamDB->SetPiZero5to10(2,1,  0.019013);
-    
-    recParamDB->SetPiZero5to10(3,0,  0.2);
-    
-    recParamDB->SetPiZero5to10(4,0,  0.252044);
-    recParamDB->SetPiZero5to10(4,1, -0.002315);
-    
-    recParamDB->SetPiZero5to10(5,0,  0.002942);
-    recParamDB->SetPiZero5to10(5,1, -3.976e-05);
-    
-    recParamDB->SetPiZero10to60(0,0,  0.009138);
-    recParamDB->SetPiZero10to60(0,1,  0.0006377);
-    
-    recParamDB->SetPiZero10to60(1,0,  1.272837);
-    recParamDB->SetPiZero10to60(1,1, -0.069708);
-    recParamDB->SetPiZero10to60(1,2,  0.001568);
-    recParamDB->SetPiZero10to60(1,3, -1.162e-05);
-    
-    recParamDB->SetPiZero10to60(2,0,  0.139703);
-    recParamDB->SetPiZero10to60(2,1,  0.003687);
-    recParamDB->SetPiZero10to60(2,2, -0.000568);
-    recParamDB->SetPiZero10to60(2,3,  1.498e-05);
-    recParamDB->SetPiZero10to60(2,4, -1.174e-07);
-    
-    recParamDB->SetPiZero10to60(3,0, -0.826367);
-    recParamDB->SetPiZero10to60(3,1,  0.096951);
-    recParamDB->SetPiZero10to60(3,2, -0.002215);
-    recParamDB->SetPiZero10to60(3,3,  2.523e-05);
-    
-    recParamDB->SetPiZero10to60(4,0,  0.249890);
-    recParamDB->SetPiZero10to60(4,1, -0.000063);
-    
-    recParamDB->SetPiZero10to60(5,0,  0.002942);
-    recParamDB->SetPiZero10to60(5,1, -3.976e-05);
-    
-    // raw signal fitting
-    recParamDB->SetHighLowGainFactor(16.);
-    recParamDB->SetOrderParameter(2);
-    recParamDB->SetTau(2.35);
-    recParamDB->SetNoiseThreshold(3);
-    recParamDB->SetNPedSamples(5);
-
-       recParamDB->SetName("Default - Pb+Pb");
-       recParamDB->SetTitle("Default - Pb+Pb");
-
+    //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetDefaultParameters();
+    AliEMCALRecParam *recParamDB = GetHighMultiplicityParameters();
+    recParamDB->SetName("Default - Pb+Pb");
+    recParamDB->SetTitle("Default - Pb+Pb");
     //Add to the recParamArray
     recParamDB->SetEventSpecie(AliRecoParam::kDefault);
     recParamArray->AddLast(recParamDB);
   }
-
+  
   //Add other options here, if desired, for
   //Cosmic, LowMult and HighMult type events
   //and add them to the array
-
-
+  
   {
     //For now, default is Pb+Pb, but let's add it again as
     //the "high mult" version too...
-    AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetDefaultParameters();
-       recParamDB->SetName("High Flux - Pb+Pb");
-       recParamDB->SetTitle("High Flux - Pb+Pb");
+    //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetHighFluxParam();
+    AliEMCALRecParam *recParamDB = GetHighMultiplicityParameters();
+    recParamDB->SetName("High Flux - Pb+Pb");
+    recParamDB->SetTitle("High Flux - Pb+Pb");
     recParamDB->SetEventSpecie(AliRecoParam::kHighMult);
     recParamArray->AddLast(recParamDB);
   }
-
+  
   {
     //Low multiplicity parameter modifications:
-    AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetDefaultParameters();
-
-    recParamDB->SetClusteringThreshold(0.2); // 200 MeV
-    recParamDB->SetMinECut(0.01);  //10 MeV
-       recParamDB->SetName("Low Flux - p+p");
-       recParamDB->SetTitle("Low Flux - p+p");
+    //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetLowFluxParam();
+    AliEMCALRecParam *recParamDB = GetLowMultiplicityParameters();
+    recParamDB->SetName("Low Flux - p+p");
+    recParamDB->SetTitle("Low Flux - p+p");
     recParamDB->SetEventSpecie(AliRecoParam::kLowMult);
     recParamArray->AddLast(recParamDB);
     
   }
-
+  
   {
-       //Cosmic parameter modifications (same as low multiplicity):
-       AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetDefaultParameters();
-               
-       recParamDB->SetClusteringThreshold(0.2); // 200 MeV
-       recParamDB->SetMinECut(0.01);  //10 MeV
-       recParamDB->SetName("Cosmic");
-       recParamDB->SetTitle("Cosmic");
-       recParamDB->SetEventSpecie(AliRecoParam::kCosmic);
-       recParamArray->AddLast(recParamDB);
-               
+    //Cosmic parameter modifications (same as low multiplicity):
+    //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetLowFluxParam();
+    AliEMCALRecParam *recParamDB = GetLowMultiplicityParameters();
+    recParamDB->SetName("Cosmic");
+    recParamDB->SetTitle("Cosmic");
+    recParamDB->SetEventSpecie(AliRecoParam::kCosmic);
+    recParamArray->AddLast(recParamDB);
+    
   }
-       
+  
   {
-       //Calib parameter modifications (same as low multiplicity):
-       AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetDefaultParameters();
-               
-       recParamDB->SetClusteringThreshold(0.2); // 200 MeV
-       recParamDB->SetMinECut(0.01);  //10 MeV
-       recParamDB->SetName("Calibration - LED");
-       recParamDB->SetTitle("Calibration - LED");
-       recParamDB->SetEventSpecie(AliRecoParam::kCalib);
-       recParamArray->AddLast(recParamDB);
-               
-       }
-
+    //Calib parameter modifications (same as low multiplicity):
+    //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetLowFluxParam();
+    AliEMCALRecParam *recParamDB = GetLowMultiplicityParameters();
+    recParamDB->SetName("Calibration - LED");
+    recParamDB->SetTitle("Calibration - LED");
+    recParamDB->SetEventSpecie(AliRecoParam::kCalib);
+    recParamArray->AddLast(recParamDB);
+    
+  }
+  
   //Set the default version in the array
   Bool_t defaultIsSet = kFALSE;
   for(Int_t i = 0; i < recParamArray->GetEntriesFast(); i++) {
@@ -239,12 +97,12 @@ void AliEMCALSetRecParamCDB(AliRecoParam::EventSpecie_t default = AliRecoParam::
       defaultIsSet = kTRUE;
     }
   }
-
+  
   if(!defaultIsSet) {
     AliError("The default reconstruction parameters are not set!  Exiting...");
     return;
   }
-
+  
   // Store calibration data into database  
   AliCDBMetaData *md = new AliCDBMetaData();
   md->SetResponsible("J. Klay");
@@ -254,6 +112,334 @@ void AliEMCALSetRecParamCDB(AliRecoParam::EventSpecie_t default = AliRecoParam::
   
   AliCDBId id("EMCAL/Calib/RecoParam",0,AliCDBRunRange::Infinity());
   cdb->GetDefaultStorage()->Put(recParamArray, id, md);
-
+  
   return;
 }
+
+//-----------------------------------------------------------------------------
+AliEMCALRecParam* GetHighMultiplicityParameters()
+{
+  //Set here the high flux/multiplicity ("Pb+Pb") parameters for the reconstruction
+  //Right now it should be the same settings as with
+  //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetHighFluxParam();
+  //or
+  //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetDefaultParameters();
+  
+  AliEMCALRecParam* params =  AliEMCALRecParam::GetDefaultParameters();
+  
+  //Clusterization
+  params->SetClusteringThreshold(0.5);
+  params->SetW0(4.5);
+  params->SetMinECut(0.45);
+  params->SetUnfold(kFALSE);
+  params->SetLocMaxCut(0.03);
+  
+  //Track matching
+  params->SetTrkCutX(6.0);
+  params->SetTrkCutY(6.0);
+  params->SetTrkCutZ(6.0);
+  params->SetTrkCutR(10.0);
+  params->SetTrkCutAlphaMin(-50.0);
+  params->SetTrkCutAlphaMax( 50.0);
+  params->SetTrkCutNITS(3.0);
+  params->SetTrkCutNTPC(20.0);
+  params->SetTrkCutAngle(10000.0);      // i.e. exclude this for the moment
+  
+  //PID
+       
+  // as a first step, all array elements are initialized to 0.0
+  Int_t i, j;
+  for (i = 0; i < 6; i++) {
+    params->SetGammaEnergyProb(i,0.);
+    params->SetHadronEnergyProb(i,0.); 
+    params->SetPiZeroEnergyProb(i,0.);
+    
+    for (j = 0; j < 6; j++) {
+      
+      params->SetGamma(i,j,0.);
+      params->SetHadron(i,j,0.);
+      params->SetPiZero(i,j, 0.);
+      params->SetGamma1to10(i,j,0.);
+      params->SetHadron1to10(i,j,0.);
+    }
+  } 
+  
+  params->SetGamma(0,0, -7.656908e-01); 
+  params->SetGamma(0,1,  2.352536e-01); 
+  params->SetGamma(0,2,  1.555996e-02);
+  params->SetGamma(0,3,  2.243525e-04);
+  params->SetGamma(0,4, -2.560087e-06);
+  
+  params->SetGamma(1,0,  6.500216e+00);
+  params->SetGamma(1,1, -2.564958e-01);
+  params->SetGamma(1,2,  1.967894e-01);
+  params->SetGamma(1,3, -3.982273e-04);
+  params->SetGamma(1,4,  2.797737e-06);
+  
+  params->SetGamma(2,0,  2.416489e+00);
+  params->SetGamma(2,1, -1.601258e-01);
+  params->SetGamma(2,2,  3.126839e-02);
+  params->SetGamma(2,3,  3.387532e-04);
+  params->SetGamma(2,4, -4.089145e-06);
+  
+  params->SetGamma(3,0,  0.);
+  params->SetGamma(3,1, -2.696008e+00);
+  params->SetGamma(3,2,  6.920305e-01);
+  params->SetGamma(3,3, -2.281122e-03);
+  params->SetGamma(3,4,  0.);
+  
+  params->SetGamma(4,0,  2.281564e-01);
+  params->SetGamma(4,1, -7.575040e-02);
+  params->SetGamma(4,2,  3.813423e-01);
+  params->SetGamma(4,3, -1.243854e-04);
+  params->SetGamma(4,4,  1.232045e-06);
+  
+  params->SetGamma(5,0, -3.290107e-01);
+  params->SetGamma(5,1,  3.707545e-02);
+  params->SetGamma(5,2,  2.917397e-03);
+  params->SetGamma(5,3,  4.695306e-05);
+  params->SetGamma(5,4, -3.572981e-07);
+  
+  params->SetHadron(0,0,   1.519112e-01);
+  params->SetHadron(0,1, -8.267603e-02);
+  params->SetHadron(0,2,  1.914574e-02);
+  params->SetHadron(0,3, -2.677921e-04);
+  params->SetHadron(0,4,  5.447939e-06);
+  
+  params->SetHadron(1,0,  0.);
+  params->SetHadron(1,1, -7.549870e-02); 
+  params->SetHadron(1,2,  3.930087e-01);
+  params->SetHadron(1,3, -2.368500e-03); 
+  params->SetHadron(1,4,  0.);
+  
+  params->SetHadron(2,0,  0.);
+  params->SetHadron(2,1, -2.463152e-02);
+  params->SetHadron(2,2,  1.349257e-01);
+  params->SetHadron(2,3, -1.089440e-03);
+  params->SetHadron(2,4,  0.);
+  
+  params->SetHadron(3,0, 0.);
+  params->SetHadron(3,1, 5.101560e-01);
+  params->SetHadron(3,2, 1.458679e-01);
+  params->SetHadron(3,3, 4.903068e-04);
+  params->SetHadron(3,4, 0.);
+  
+  params->SetHadron(4,0, 0.);
+  params->SetHadron(4,1, -6.693943e-03); 
+  params->SetHadron(4,2,  2.444753e-01);
+  params->SetHadron(4,3, -5.553749e-05);
+  params->SetHadron(4,4, 0.);
+  
+  params->SetHadron(5,0, -4.414030e-01);
+  params->SetHadron(5,1, 2.292277e-01);
+  params->SetHadron(5,2, -2.433737e-02);
+  params->SetHadron(5,3,  1.758422e-03);
+  params->SetHadron(5,4, -3.001493e-05);
+  
+  params->SetPiZero(0,0,  5.072157e-01);
+  params->SetPiZero(0,1, -5.352747e-01);
+  params->SetPiZero(0,2,  8.499259e-02);
+  params->SetPiZero(0,3, -3.687401e-03);
+  params->SetPiZero(0,4,  5.482280e-05);
+  
+  params->SetPiZero(1,0,  4.590137e+02); 
+  params->SetPiZero(1,1, -7.079341e+01);
+  params->SetPiZero(1,2,  4.990735e+00);
+  params->SetPiZero(1,3, -1.241302e-01);
+  params->SetPiZero(1,4,  1.065772e-03);
+  
+  params->SetPiZero(2,0,  1.376415e+02);
+  params->SetPiZero(2,1, -3.031577e+01);
+  params->SetPiZero(2,2,  2.474338e+00);
+  params->SetPiZero(2,3, -6.903410e-02);
+  params->SetPiZero(2,4,  6.244089e-04);
+  
+  params->SetPiZero(3,0, 0.);
+  params->SetPiZero(3,1,  1.145983e+00);
+  params->SetPiZero(3,2, -2.476052e-01);
+  params->SetPiZero(3,3,  1.367373e-02);
+  params->SetPiZero(3,4, 0.);
+  
+  params->SetPiZero(4,0, -2.097586e+02);
+  params->SetPiZero(4,1,  6.300800e+01);
+  params->SetPiZero(4,2, -4.038906e+00);
+  params->SetPiZero(4,3,  1.088543e-01);
+  params->SetPiZero(4,4, -9.362485e-04);
+  
+  params->SetPiZero(5,0, -1.671477e+01); 
+  params->SetPiZero(5,1,  2.995415e+00);
+  params->SetPiZero(5,2, -6.040360e-02);
+  params->SetPiZero(5,3, -6.137459e-04);
+  params->SetPiZero(5,4,  1.847328e-05);
+  
+  // High flux ones pp 
+  
+  params->SetHadronEnergyProb(0, 0.);
+  params->SetHadronEnergyProb(1, 0.);
+  params->SetHadronEnergyProb(2,  6.188452e-02);
+  params->SetHadronEnergyProb(3,  2.030230e+00);
+  params->SetHadronEnergyProb(4, -6.402242e-02);
+  
+  // raw signal fitting
+  params->SetHighLowGainFactor(16.);
+  params->SetOrderParameter(2);
+  params->SetTau(2.35);
+  params->SetNoiseThreshold(3);
+  params->SetNPedSamples(5);
+  
+  return params ;
+}      
+
+//-----------------------------------------------------------------------------
+AliEMCALRecParam* GetLowMultiplicityParameters()
+{
+  // Set here the low flux/multiplicity ("p+p") parameters for the reconstruction
+  //Right now it should be the same settings as with
+  //AliEMCALRecParam *recParamDB = AliEMCALRecParam::GetLowFluxParam();
+  
+  AliEMCALRecParam* params =  AliEMCALRecParam::GetDefaultParameters();
+  params->SetClusteringThreshold(0.2); // 200 MeV                                             
+  params->SetMinECut(0.01);  //10 MeV          
+  
+  //PID parameters for pp  implemented 
+  // as a first step, all array elements are initialized to 0.0
+  Int_t i, j;
+  for (i = 0; i < 6; i++) {
+    for (j = 0; j < 6; j++) {
+      params->SetGamma(i,j,0.);
+      params->SetGamma1to10(i,j,0.);
+      params->SetHadron(i,j,0.);
+      params->SetHadron1to10(i,j,0.);
+      params->SetPiZero(i,j,0.);
+      
+    }
+    params->SetGammaEnergyProb(i,0.); // not yet implemented
+    params->SetHadronEnergyProb(i,0.);
+    params->SetPiZeroEnergyProb(i,0.); // not yet implemented
+  }
+  
+  params->SetGamma(0,0, -7.656908e-01);
+  params->SetGamma(0,1,  2.352536e-01);
+  params->SetGamma(0,2,  1.555996e-02);
+  params->SetGamma(0,3,  2.243525e-04);
+  params->SetGamma(0,4, -2.560087e-06);
+  
+  params->SetGamma(1,0,  6.500216e+00);
+  params->SetGamma(1,1, -2.564958e-01);
+  params->SetGamma(1,2,  1.967894e-01);
+  params->SetGamma(1,3, -3.982273e-04);
+  params->SetGamma(1,4,  2.797737e-06);
+  
+  params->SetGamma(2,0,  2.416489e+00);
+  params->SetGamma(2,1, -1.601258e-01);
+  params->SetGamma(2,2,  3.126839e-02);
+  params->SetGamma(2,3,  3.387532e-04);
+  params->SetGamma(2,4, -4.089145e-06);
+  
+  params->SetGamma(3,0,0.);
+  params->SetGamma(3,1,-2.696008e+00); 
+  params->SetGamma(3,2, 6.920305e-01);
+  params->SetGamma(3,3,-2.281122e-03);
+  params->SetGamma(3,4,0.);
+  
+  params->SetGamma(4,0,  2.281564e-01); 
+  params->SetGamma(4,1, -7.575040e-02);
+  params->SetGamma(4,2,  3.813423e-01);
+  params->SetGamma(4,3, -1.243854e-04);
+  params->SetGamma(4,4,  1.232045e-06);
+  
+  params->SetGamma(5,0, -3.290107e-01);
+  params->SetGamma(5,1,  3.707545e-02);
+  params->SetGamma(5,2,  2.917397e-03);
+  params->SetGamma(5,3,  4.695306e-05);
+  params->SetGamma(5,4, -3.572981e-07);
+  
+  params->SetHadron(0,0,  9.482243e-01); 
+  params->SetHadron(0,1, -2.780896e-01);
+  params->SetHadron(0,2,  2.223507e-02);
+  params->SetHadron(0,3,  7.294263e-04);
+  params->SetHadron(0,4, -5.665872e-06); 
+  
+  params->SetHadron(1,0,  0.);
+  params->SetHadron(1,1,  0.);
+  params->SetHadron(1,2,  2.483298e-01);
+  params->SetHadron(1,3,  0.);
+  params->SetHadron(1,4,  0.);
+  
+  params->SetHadron(2,0, -5.601199e+00);
+  params->SetHadron(2,1,  2.097382e+00);
+  params->SetHadron(2,2, -2.307965e-01);
+  params->SetHadron(2,3,  9.206871e-03);
+  params->SetHadron(2,4, -8.887548e-05);
+  
+  params->SetHadron(3,0,  6.543101e+00);
+  params->SetHadron(3,1, -2.305203e+00);
+  params->SetHadron(3,2,  2.761673e-01);
+  params->SetHadron(3,3, -5.465855e-03);
+  params->SetHadron(3,4,  2.784329e-05);
+  
+  params->SetHadron(4,0, -2.443530e+01);
+  params->SetHadron(4,1,  8.902578e+00);
+  params->SetHadron(4,2, -5.265901e-01);
+  params->SetHadron(4,3,  2.549111e-02);
+  params->SetHadron(4,4, -2.196801e-04);
+  
+  params->SetHadron(5,0,  2.102007e-01);
+  params->SetHadron(5,1, -3.844418e-02);
+  params->SetHadron(5,2,  1.234682e-01);
+  params->SetHadron(5,3, -3.866733e-03);
+  params->SetHadron(5,4,  3.362719e-05);
+  
+  params->SetPiZero(0,0,  5.07215e-01);
+  params->SetPiZero(0,1, -5.35274e-01);
+  params->SetPiZero(0,2,  8.49925e-02);
+  params->SetPiZero(0,3, -3.68740e-03);
+  params->SetPiZero(0,4,  5.48228e-05);
+  
+  params->SetPiZero(1,0,  4.590137e+02);
+  params->SetPiZero(1,1, -7.079341e+01);
+  params->SetPiZero(1,2,  4.990735e+00);
+  params->SetPiZero(1,3, -1.241302e-01);
+  params->SetPiZero(1,4,  1.065772e-03);
+  
+  params->SetPiZero(2,0,  1.376415e+02); 
+  params->SetPiZero(2,1, -3.031577e+01);
+  params->SetPiZero(2,2,  2.474338e+00);
+  params->SetPiZero(2,3, -6.903410e-02); 
+  params->SetPiZero(2,4,  6.244089e-04);
+  
+  params->SetPiZero(3,0,  0.);
+  params->SetPiZero(3,1,  1.145983e+00);
+  params->SetPiZero(3,2, -2.476052e-01);
+  params->SetPiZero(3,3,  1.367373e-02);
+  params->SetPiZero(3,4,  0.);
+  
+  params->SetPiZero(4,0, -2.09758e+02);
+  params->SetPiZero(4,1,  6.30080e+01);
+  params->SetPiZero(4,2, -4.03890e+00);
+  params->SetPiZero(4,3,  1.08854e-01);
+  params->SetPiZero(4,4, -9.36248e-04);
+  
+  params->SetPiZero(5,0, -1.671477e+01);
+  params->SetPiZero(5,1,  2.995415e+00);
+  params->SetPiZero(5,2, -6.040360e-02);
+  params->SetPiZero(5,3,  -6.137459e-04);
+  params->SetPiZero(5,4, 1.847328e-05);
+  
+  //     params->SetHadronEnergyProb(0,0.);
+  //     params->SetHadronEnergyProb(1,0.);
+  //     params->SetHadronEnergyProb(2,1.);
+  //     params->SetHadronEnergyProb(3,0.);
+  //     params->SetHadronEnergyProb(4,0.);
+  
+  params->SetHadronEnergyProb(0,  4.767543e-02);
+  params->SetHadronEnergyProb(1, -1.537523e+00);
+  params->SetHadronEnergyProb(2,  2.956727e-01);
+  params->SetHadronEnergyProb(3, -3.051022e+01);
+  params->SetHadronEnergyProb(4, -6.036931e-02);
+  
+  return params;
+  
+}
+
+
index 51cde3d..cb58f48 100644 (file)
Binary files a/OCDB/EMCAL/Calib/RecoParam/Run0_999999999_v0_s0.root and b/OCDB/EMCAL/Calib/RecoParam/Run0_999999999_v0_s0.root differ