[u/mrichter/AliRoot.git] / EMCAL / AliEMCALPID.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */
/* History of cvs commits:
 *
 * $Log$
 * 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
//
//
//
//
//
//	AliEMCALPID::CalculPID(Energy,Lambda0)
//	  Calcul PID for all clusters in AliESDs.root file
//	  keep this function for the moment for a simple verification, could be removed
//
//
//
//   AliEMCALPID::CalculPID(Energy,Lambda0)
//    calcul PID Weght for a cluster with Energy, Lambda0 .
//    Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
//
//
//
//	kElectron :  fPIDFinal[0]
//	kMuon     :  fPIDFinal[1]
//	kPion	  :  fPIDFinal[2]
//	kKaon	  :  fPIDFinal[3]
//	kProton   :  fPIDFinal[4]
//	kPhoton   :  fPIDFinal[5]
//	kPi0	  :  fPIDFinal[6]
//	kNeutron  :  fPIDFinal[7]
//	kKaon0	  :  fPIDFinal[8]
//	kEleCon   :  fPIDFinal[9]
//	kUnknown  :  fPIDFinal[10]
//
//
//    PID[3] is a simple PID for
//      Electron & Photon  PID[0]
//	              Pi0  PID[1]
//		   Hadron  PID[2]
//
//  
//
//
//
// --- ROOT system ---

// standard C++ includes
#include <Riostream.h>
// #include <cstdlib>
// #include <cmath>

// ROOT includes
#include "TTree.h"
#include "TStyle.h"
#include "TVector3.h"
#include "TBranch.h"
#include "TClonesArray.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TParticle.h"

// // STEER includes
// #include "AliRun.h"
// #include "AliRunLoader.h"
// #include "AliHeader.h"
// #include "AliLoader.h"
// #include "AliStack.h"
// #include "AliESDtrack.h"
// #include "AliESD.h"
#include "AliLog.h"
#include "AliEMCALPID.h"
  
ClassImp(AliEMCALPID)

//______________________________________________
  AliEMCALPID::AliEMCALPID():
    fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE)
{
  //
  // Constructor.
  // Initialize all constant values which have to be used
  // during PID algorithm execution
  //
  
  // set flag for printing to FALSE by default
  fPrintInfo = kFALSE;
  
  // 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.;
    }
  }
  
  // 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;
  
  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;
  
  fPiZero10to60[4][0] =  0.249890;
  fPiZero10to60[4][1] = -0.000063;
  
  fPiZero10to60[5][0] =  0.002942;
  fPiZero10to60[5][1] = -3.976e-05;
  
  fPIDWeight[0] = -1;
  fPIDWeight[1] = -1;
  fPIDWeight[2] = -1;
  fReconstructor = kFALSE;
}

//______________________________________________
void AliEMCALPID::RunPID(AliESD *esd)
{
//
// 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->GetNumberOfEMCALClusters();
  Int_t firstCluster = esd->GetFirstEMCALCluster();
  Double_t energy, lambda0;
  for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
    
    AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
    energy = clust->E();
    lambda0 = clust->GetM02();
    // verify cluster type
    Int_t clusterType= clust->GetClusterType();
    if (clusterType == AliESDCaloCluster::kClusterv1 && lambda0 != 0  && energy < 1000) {


      // reject clusters with lambda0 = 0


      ComputePID(energy, lambda0);


      if (fPrintInfo) {
	AliInfo("___________________________________________________");
	AliInfo(Form( "Particle Energy = %f",energy));
	AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
	AliInfo("PIDWeight of the particle :" );
	AliInfo(Form( " GAMMA  : %f",fPID[0] ));
	AliInfo(Form( " PiZero : %f",fPID[1] ));
	AliInfo(Form( " HADRON : %f", fPID[2] ));
	AliInfo("_________________________________________");
	AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
	AliInfo(Form( " kMuon     : %f", fPIDFinal[1] ));
	AliInfo(Form( " kPion	    : %f", fPIDFinal[2] ));
	AliInfo(Form( " kKaon	    : %f", fPIDFinal[3] ));
	AliInfo(Form( " kProton   : %f", fPIDFinal[4] ));
	AliInfo(Form( " kPhoton   : %f", fPIDFinal[5] ));
	AliInfo(Form( " kPi0	    : %f", fPIDFinal[6] ));
	AliInfo(Form( " kNeutron  : %f", fPIDFinal[7] ));
	AliInfo(Form( " kKaon0	: %f", fPIDFinal[8] ));
	AliInfo(Form( " kEleCon   : %f", fPIDFinal[9] ));
	AliInfo(Form( " kUnknown  : %f", fPIDFinal[10] ));
	AliInfo("___________________________________________________");
      }

      if(fReconstructor) // In case it is called during reconstruction.
	clust->SetPid(fPIDFinal);
    } // end if (clusterType...)
  } // end for (iCluster...)
}

//__________________________________________________________
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;}


  TArrayD paramDistribGamma  = DistLambda0(energy, 1);
  TArrayD paramDistribPiZero = DistLambda0(energy, 2);
  TArrayD paramDistribHadron = DistLambda0(energy, 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];
  fProbHadron  = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
  fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
  
  // compute PID Weight
  fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
  fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
  fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
  
  SetPID(fPIDWeight[0], 0);
  SetPID(fPIDWeight[1], 1);
  SetPID(fPIDWeight[2], 2);
  
  // sortie ecran 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[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[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)
{
  //
  // Compute the values of the parametrised distributions using the data initialised before.
  //
  Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
  Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
  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]);
    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]);
    }
    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]);
    break;
  }
  
  distributionParam[0] = constGauss;
  distributionParam[1] = meanGauss;
  distributionParam[2] = sigmaGauss;
  distributionParam[3] = constLandau;
  distributionParam[4] = mpvLandau;
  distributionParam[5] = sigmaLandau;
  
  return distributionParam;
}

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