**************************************************************************/
/* $Id$ */
-/* History of cvs commits:
- *
- * $Log$
- * 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
+
+// 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
+//
+// 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::kSPECIESCN] is the standard PID for :
//
// kElectron : fPIDFinal[0]
// kMuon : fPIDFinal[1]
-// kPion : fPIDFinal[2]
-// kKaon : fPIDFinal[3]
+// kPion : fPIDFinal[2]
+// kKaon : fPIDFinal[3]
// kProton : fPIDFinal[4]
// kPhoton : fPIDFinal[5]
-// kPi0 : fPIDFinal[6]
+// kPi0 : fPIDFinal[6]
// kNeutron : fPIDFinal[7]
// kKaon0 : fPIDFinal[8]
// kEleCon : fPIDFinal[9]
// Pi0 PID[1]
// Hadron PID[2]
//
-//
-//
-//
-//
-// --- ROOT system ---
+// Author: Genole Bourdaud 2007 (SUBATECH)
+// Marie Germain 07/2009 (SUBATECH), new parametrization for low and high flux environment
+// Gustavo Conesa 08/2009 (LNF), divide class in AliEMCALPID and AliEMCALPIDUtils, PIDUtils belong to library EMCALUtils
+// --- standard c ---
// standard C++ includes
-#include <Riostream.h>
-// #include <cstdlib>
-// #include <cmath>
+//#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 "AliRun.h"
-// #include "AliRunLoader.h"
-// #include "AliHeader.h"
-// #include "AliLoader.h"
-// #include "AliStack.h"
-// #include "AliESDtrack.h"
-// #include "AliESDEvent.h"
-#include "AliLog.h"
+// STEER includes
+#include "AliESDEvent.h"
#include "AliEMCALPID.h"
#include "AliESDCaloCluster.h"
+#include "AliEMCALReconstructor.h"
+
ClassImp(AliEMCALPID)
-
+
//______________________________________________
- AliEMCALPID::AliEMCALPID():
- fPrintInfo(kFALSE), fProbGamma(0.),fProbPiZero(0.),fProbHadron(0.),fReconstructor(kFALSE)
+ AliEMCALPID::AliEMCALPID()
+ : AliEMCALPIDUtils(), fReconstructor(kTRUE)
{
//
// Constructor.
// 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;
+ InitParameters();
- 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;
+}
+
+//______________________________________________
+AliEMCALPID::AliEMCALPID(Bool_t reconstructor)
+: AliEMCALPIDUtils(), fReconstructor(reconstructor)
+{
+ //
+ // Constructor.
+ // Initialize all constant values which have to be used
+ // during PID algorithm execution called when used in standalone mode
+ //
- fPiZero10to60[5][0] = 0.002942;
- fPiZero10to60[5][1] = -3.976e-05;
+ InitParameters();
- fPIDWeight[0] = -1;
- fPIDWeight[1] = -1;
- fPIDWeight[2] = -1;
- fReconstructor = kFALSE;
}
//______________________________________________
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->GetNumberOfEMCALClusters();
- Int_t firstCluster = esd->GetFirstEMCALCluster();
- Double_t energy, lambda0;
+
+ Int_t nClusters = esd->GetNumberOfCaloClusters();
+ Int_t firstCluster = 0;
+ Double_t energy=0., lambda0=0.;
for (Int_t iCluster = firstCluster; iCluster < (nClusters + firstCluster); iCluster++) {
AliESDCaloCluster *clust = esd->GetCaloCluster(iCluster);
+ 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));
AliInfo(Form( " kUnknown : %f", fPIDFinal[10] ));
AliInfo("___________________________________________________");
}
-
- if(fReconstructor) // In case it is called during reconstruction.
- clust->SetPid(fPIDFinal);
- } // end if (clusterType...)
+
+ if(fReconstructor){ // In case it is called during reconstruction.
+ // cout << "############# Fill ESDs with PIDWeight ##########" << endl;
+ clust->SetPID(fPIDFinal);}
+ } // end if (lambda0 != 0 && energy < 1000)
} // 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);
+//_______________________________________________________
+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.;
- SetPID(fPIDWeight[0], 0);
- SetPID(fPIDWeight[1], 1);
- SetPID(fPIDWeight[2], 2);
+ fPIDWeight[0] = -1;
+ fPIDWeight[1] = -1;
+ fPIDWeight[2] = -1;
- // 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("********************************************************" );
- }
+ for(Int_t i=0; i<AliPID::kSPECIESCN+1; i++)
+ fPIDFinal[i]= 0;
- 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);
+ const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
- 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]);
+ if(fReconstructor){
+
+ if(!recParam) {
+ AliFatal("Reconstruction parameters for EMCAL not set!");
}
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;
- }
+
+ 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
+
+ }
- distributionParam[0] = constGauss;
- distributionParam[1] = meanGauss;
- distributionParam[2] = sigmaGauss;
- distributionParam[3] = constLandau;
- distributionParam[4] = mpvLandau;
- distributionParam[5] = sigmaLandau;
+ else{
+ // init the parameters here instead of from loading from recparam
+ // default parameters are PbPb parameters.
+ SetHighFluxParam();
+
+ }
- 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;
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff