[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.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
//
//
//
//
//
//	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>

// 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 "AliLog.h"
#include "AliEMCALPID.h"
#include "AliESDCaloCluster.h"
#include "AliEMCALRecParam.h"
#include "AliEMCALReconstructor.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
  //
 
  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(!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] ));
      }
    }
    
  }
  
}

//______________________________________________
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

  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::kEMCALClusterv1 && 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$
46363b15	20	* Revision 1.16 2007/11/23 13:39:05 gustavo
	21	* Track matching and PID parameters added to AliEMCALRecParam
	22	*
8ba062b1	23	* Revision 1.15 2007/10/09 08:46:10 hristov
	24	* 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)
	25	*
8ada0ffe	26	* Revision 1.14 2007/07/26 16:54:53 morsch
	27	* Changes in AliESDEvent fwd declarartions.
	28	*
89ffc0b0	29	* Revision 1.13 2007/07/11 13:43:29 hristov
	30	* New class AliESDEvent, backward compatibility with the old AliESD (Christian)
	31	*
af885e0f	32	* Revision 1.12 2007/06/11 20:43:06 hristov
	33	* Changes required by the updated AliESDCaloCluster (Gustavo)
	34	*
7592dfc4	35	* Revision 1.11 2007/03/30 13:50:34 gustavo
	36	* PID for particles with E < 5 GeV was not done, temporal solution found (Guenole)
	37	*
4ce73a00	38	* Revision 1.10 2007/03/09 14:34:11 gustavo
	39	* Correct probability calculation, added missing initialization of data members
	40	*
4cfd8ef9	41	* Revision 1.9 2007/02/20 20:17:43 hristov
	42	* Corrected array size, removed warnings (icc)
	43	*
d69ab345	44	* Revision 1.8 2006/12/19 08:49:35 gustavo
	45	* New PID class for EMCAL, bayesian analysis done with ESD data, PID information filled when calling AliEMCALPID in AliEMCALReconstructor::FillESD()
	46	*
dc293ae9	47	*
	48	*/
	49	// to compute PID for all the clusters in ESDs.root file
	50	// the ESDs.root have to be in the same directory as the class
	51	//
	52	//
	53	//
	54	//
	55	//
	56	// AliEMCALPID::CalculPID(Energy,Lambda0)
	57	// Calcul PID for all clusters in AliESDs.root file
	58	// keep this function for the moment for a simple verification, could be removed
	59	//
	60	//
	61	//
	62	// AliEMCALPID::CalculPID(Energy,Lambda0)
	63	// calcul PID Weght for a cluster with Energy, Lambda0 .
	64	// Double_t PIDFinal[AliPID::kSPECIESN] is the standard PID for :
	65	//
	66	//
	67	//
	68	// kElectron : fPIDFinal[0]
	69	// kMuon : fPIDFinal[1]
	70	// kPion : fPIDFinal[2]
	71	// kKaon : fPIDFinal[3]
	72	// kProton : fPIDFinal[4]
	73	// kPhoton : fPIDFinal[5]
	74	// kPi0 : fPIDFinal[6]
	75	// kNeutron : fPIDFinal[7]
	76	// kKaon0 : fPIDFinal[8]
	77	// kEleCon : fPIDFinal[9]
	78	// kUnknown : fPIDFinal[10]
	79	//
	80	//
	81	// PID[3] is a simple PID for
	82	// Electron & Photon PID[0]
	83	// Pi0 PID[1]
	84	// Hadron PID[2]
	85	//
	86	//
	87	//
	88	//
	89	//
	90	// --- ROOT system ---
	91
	92	// standard C++ includes
	93	#include <Riostream.h>
dc293ae9	94
	95	// ROOT includes
	96	#include "TTree.h"
	97	#include "TStyle.h"
	98	#include "TVector3.h"
	99	#include "TBranch.h"
	100	#include "TClonesArray.h"
	101	#include "TCanvas.h"
	102	#include "TLorentzVector.h"
	103	#include "TMath.h"
	104	#include "TFile.h"
	105	#include "TH1.h"
	106	#include "TH2.h"
	107	#include "TParticle.h"
	108
46363b15	109	// STEER includes
dc293ae9	110	#include "AliLog.h"
dc293ae9	111	#include "AliEMCALPID.h"
89ffc0b0	112	#include "AliESDCaloCluster.h"
8ba062b1	113	#include "AliEMCALRecParam.h"
	114	#include "AliEMCALReconstructor.h"
	115
dc293ae9	116
	117	ClassImp(AliEMCALPID)
	118
4cfd8ef9	119	//______________________________________________
	120	AliEMCALPID::AliEMCALPID():
	121	fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE)
dc293ae9	122	{
4cfd8ef9	123	//
	124	// Constructor.
	125	// Initialize all constant values which have to be used
	126	// during PID algorithm execution
	127	//
46363b15	128
4cfd8ef9	129	fPIDWeight[0] = -1;
	130	fPIDWeight[1] = -1;
	131	fPIDWeight[2] = -1;
8ba062b1	132
46363b15	133	for(Int_t i=0; i<AliPID::kSPECIESN+1; i++)
	134	fPIDFinal[i]= 0;
	135
	136	const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
	137	if(!recParam) {
8ba062b1	138	AliFatal("Reconstruction parameters for EMCAL not set!");
	139	}
	140	else {
46363b15	141	for(Int_t i=0; i<6; i++){
	142	for(Int_t j=0; j<6; j++){
	143	fGamma[i][j] = recParam->GetGamma(i,j);
	144	fHadron[i][j] = recParam->GetHadron(i,j);
	145	fPiZero5to10[i][j] = recParam->GetPiZero5to10(i,j);
	146	fPiZero10to60[i][j] = recParam->GetPiZero10to60(i,j);
	147	AliDebug(1,Form("PID parameters (%d, %d): fGamma=%.3f, fPi=%.3f, fHadron=%.3f",
	148	i,j, fGamma[i][j],fPiZero5to10[i][j],fHadron[i][j] ));
	149	}
	150	}
	151
8ba062b1	152	}
46363b15	153
dc293ae9	154	}
4cfd8ef9	155
4cfd8ef9	156	//______________________________________________
af885e0f	157	void AliEMCALPID::RunPID(AliESDEvent *esd)
dc293ae9	158	{
	159	//
	160	// Make the PID for all the EMCAL clusters containedin the ESDs File
	161	// but just gamma/PiO/Hadron
	162	//
	163	// trivial check against NULL object passed
8ba062b1	164
dc293ae9	165	if (esd == 0x0) {
4ce73a00	166	AliInfo("NULL ESD object passed !!" );
dc293ae9	167	return ;
dc293ae9	168	}
4cfd8ef9	169
dc293ae9	170	Int_t nClusters = esd->GetNumberOfEMCALClusters();
	171	Int_t firstCluster = esd->GetFirstEMCALCluster();
	172	Double_t energy, lambda0;
	173	for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
8ba062b1	174
dc293ae9	175	AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
7592dfc4	176	energy = clust->E();
dc293ae9	177	lambda0 = clust->GetM02();
	178	// verify cluster type
	179	Int_t clusterType= clust->GetClusterType();
8ada0ffe	180	if (clusterType == AliESDCaloCluster::kEMCALClusterv1 && lambda0 != 0 && energy < 1000) {
4ce73a00	181
4ce73a00	182
dc293ae9	183	// reject clusters with lambda0 = 0
4ce73a00	184
4ce73a00	185
dc293ae9	186	ComputePID(energy, lambda0);
4ce73a00	187
4ce73a00	188
dc293ae9	189	if (fPrintInfo) {
	190	AliInfo("___________________________________________________");
	191	AliInfo(Form( "Particle Energy = %f",energy));
	192	AliInfo(Form( "Particle Lambda0 of the particle = %f", lambda0) );
	193	AliInfo("PIDWeight of the particle :" );
	194	AliInfo(Form( " GAMMA : %f",fPID[0] ));
	195	AliInfo(Form( " PiZero : %f",fPID[1] ));
	196	AliInfo(Form( " HADRON : %f", fPID[2] ));
	197	AliInfo("_________________________________________");
	198	AliInfo(Form( " kElectron : %f", fPIDFinal[0]) );
	199	AliInfo(Form( " kMuon : %f", fPIDFinal[1] ));
	200	AliInfo(Form( " kPion : %f", fPIDFinal[2] ));
	201	AliInfo(Form( " kKaon : %f", fPIDFinal[3] ));
	202	AliInfo(Form( " kProton : %f", fPIDFinal[4] ));
	203	AliInfo(Form( " kPhoton : %f", fPIDFinal[5] ));
	204	AliInfo(Form( " kPi0 : %f", fPIDFinal[6] ));
	205	AliInfo(Form( " kNeutron : %f", fPIDFinal[7] ));
	206	AliInfo(Form( " kKaon0 : %f", fPIDFinal[8] ));
	207	AliInfo(Form( " kEleCon : %f", fPIDFinal[9] ));
	208	AliInfo(Form( " kUnknown : %f", fPIDFinal[10] ));
	209	AliInfo("___________________________________________________");
	210	}
4cfd8ef9	211
dc293ae9	212	if(fReconstructor) // In case it is called during reconstruction.
	213	clust->SetPid(fPIDFinal);
	214	} // end if (clusterType...)
	215	} // end for (iCluster...)
	216	}
4cfd8ef9	217
4cfd8ef9	218	//__________________________________________________________
dc293ae9	219	void AliEMCALPID::ComputePID(Double_t energy, Double_t lambda0)
	220	{
	221	//
	222	// This is the main command, which uses the distributions computed and parametrised,
	223	// and gives the PID by the bayesian method.
	224	//
4ce73a00	225
	226	if (energy<5){energy =6;}
	227
	228
4cfd8ef9	229	TArrayD paramDistribGamma = DistLambda0(energy, 1);
	230	TArrayD paramDistribPiZero = DistLambda0(energy, 2);
	231	TArrayD paramDistribHadron = DistLambda0(energy, 3);
	232
	233	Bool_t norm = kFALSE;
	234
	235	fProbGamma = TMath::Gaus(lambda0, paramDistribGamma[1], paramDistribGamma[2], norm) * paramDistribGamma[0];
	236	fProbGamma += TMath::Landau(lambda0, paramDistribGamma[4], paramDistribGamma[5], norm) * paramDistribGamma[3];
	237	fProbPiZero = TMath::Gaus(lambda0, paramDistribPiZero[1], paramDistribPiZero[2], norm) * paramDistribPiZero[0];
	238	fProbPiZero += TMath::Landau(lambda0, paramDistribPiZero[4], paramDistribPiZero[5], norm) * paramDistribPiZero[3];
	239	fProbHadron = TMath::Gaus(lambda0, paramDistribHadron[1], paramDistribHadron[2], norm) * paramDistribHadron[0];
	240	fProbHadron += TMath::Landau(lambda0, paramDistribHadron[4], paramDistribHadron[5], norm) * paramDistribHadron[3];
	241
	242	// compute PID Weight
	243	fPIDWeight[0] = fProbGamma / (fProbGamma + fProbPiZero + fProbHadron);
	244	fPIDWeight[1] = fProbPiZero / (fProbGamma+fProbPiZero+fProbHadron);
	245	fPIDWeight[2] = fProbHadron / (fProbGamma+fProbPiZero+fProbHadron);
	246
	247	SetPID(fPIDWeight[0], 0);
	248	SetPID(fPIDWeight[1], 1);
	249	SetPID(fPIDWeight[2], 2);
	250
	251	// sortie ecran pid Weight only for control (= in english ???)
	252	if (fPrintInfo) {
	253	AliInfo(Form( "Energy in loop = %f", energy) );
	254	AliInfo(Form( "Lambda0 in loop = %f", lambda0) );
	255	AliInfo(Form( "fProbGamma in loop = %f", fProbGamma) );
	256	// AliInfo(Form( "fParametresDistribGamma[2] = %f", fParamDistribGamma[2]) );
	257	AliInfo(Form( "fProbaPiZero = %f", fProbPiZero ));
	258	AliInfo(Form( "fProbaHadron = %f", fProbHadron) );
	259	AliInfo(Form( "PIDWeight in loop = %f \|\|\| %f \|\|\| %f", fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
	260	AliInfo(Form( "fGamma[2][2] = %f", fGamma[2][2] ));
	261	AliInfo("********************************************************" );
	262	}
	263
	264	fPIDFinal[0] = fPIDWeight[0]/2;
	265	fPIDFinal[1] = fPIDWeight[2]/8;
	266	fPIDFinal[2] = fPIDWeight[2]/8;
	267	fPIDFinal[3] = fPIDWeight[2]/8;
	268	fPIDFinal[4] = fPIDWeight[2]/8;
	269	fPIDFinal[5] = fPIDWeight[0]/2;
	270	fPIDFinal[6] = fPIDWeight[1] ;
	271	fPIDFinal[7] = fPIDWeight[2]/8;
	272	fPIDFinal[8] = fPIDWeight[2]/8;
	273	fPIDFinal[9] = fPIDWeight[2]/8;
	274	fPIDFinal[10] = fPIDWeight[2]/8;
dc293ae9	275	}
4cfd8ef9	276
4cfd8ef9	277	//________________________________________________________
dc293ae9	278	TArrayD AliEMCALPID::DistLambda0(Double_t energy, Int_t type)
dc293ae9	279	{
4cfd8ef9	280	//
	281	// Compute the values of the parametrised distributions using the data initialised before.
	282	//
	283	Double_t constGauss = 0., meanGauss = 0., sigmaGauss = 0.;
	284	Double_t constLandau=0., mpvLandau=0., sigmaLandau=0.;
	285	TArrayD distributionParam(6);
	286
	287	switch (type) {
	288	case 1:
	289	constGauss = Polynomial(energy, fGamma[0]);
	290	meanGauss = Polynomial(energy, fGamma[1]);
	291	sigmaGauss = Polynomial(energy, fGamma[2]);
	292	constLandau = Polynomial(energy, fGamma[3]);
	293	mpvLandau = Polynomial(energy, fGamma[4]);
	294	sigmaLandau = Polynomial(energy, fGamma[5]);
46363b15	295
46363b15	296
4cfd8ef9	297	break;
	298	case 2:
	299	if (energy < 10) {
	300	constGauss = Polynomial(energy, fPiZero5to10[0]);
	301	meanGauss = Polynomial(energy, fPiZero5to10[1]);
	302	sigmaGauss = Polynomial(energy, fPiZero5to10[2]);
	303	constLandau = Polynomial(energy, fPiZero5to10[3]);
	304	mpvLandau = Polynomial(energy, fPiZero5to10[4]);
	305	sigmaLandau = Polynomial(energy, fPiZero5to10[5]);
46363b15	306
4cfd8ef9	307	}
	308	else {
	309	constGauss = Polynomial(energy, fPiZero10to60[0]);
	310	meanGauss = Polynomial(energy, fPiZero10to60[1]);
	311	sigmaGauss = Polynomial(energy, fPiZero10to60[2]);
	312	constLandau = Polynomial(energy, fPiZero10to60[3]);
	313	mpvLandau = Polynomial(energy, fPiZero10to60[4]);
	314	sigmaLandau = Polynomial(energy, fPiZero10to60[5]);
46363b15	315
4cfd8ef9	316	}
	317	break;
	318	case 3:
	319	constGauss = Polynomial(energy, fHadron[0]);
	320	meanGauss = Polynomial(energy, fHadron[1]);
	321	sigmaGauss = Polynomial(energy, fHadron[2]);
	322	constLandau = Polynomial(energy, fHadron[3]);
	323	mpvLandau = Polynomial(energy, fHadron[4]);
	324	sigmaLandau = Polynomial(energy, fHadron[5]);
46363b15	325
4cfd8ef9	326	break;
	327	}
	328
	329	distributionParam[0] = constGauss;
	330	distributionParam[1] = meanGauss;
	331	distributionParam[2] = sigmaGauss;
	332	distributionParam[3] = constLandau;
	333	distributionParam[4] = mpvLandau;
	334	distributionParam[5] = sigmaLandau;
	335
	336	return distributionParam;
dc293ae9	337	}
4cfd8ef9	338
4cfd8ef9	339	//_______________________________________________________
dc293ae9	340	Double_t AliEMCALPID::Polynomial(Double_t x, Double_t *params)
dc293ae9	341	{
4cfd8ef9	342	//
	343	// Compute a polynomial for a given value of 'x'
	344	// with the array of parameters passed as the second arg
	345	//
46363b15	346
4cfd8ef9	347	Double_t y;
	348	y = params[0];
	349	y += params[1] * x;
	350	y += params[2] * x * x;
	351	y += params[3] * x * x * x;
	352	y += params[4] * x * x * x * x;
	353	y += params[5] * x * x * x * x * x;
	354
	355	return y;
dc293ae9	356	}