X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALPID.cxx;h=71cd6631c4ab9f0204673267d40304e550855536;hb=54d34aac48bb923502ca39356d982e5f84398da5;hp=60be817782690760f6fd187a3eef0cdee310a585;hpb=46363b1507c2550f9c5e71c0fbf0bb06533cf991;p=u%2Fmrichter%2FAliRoot.git

diff --git a/EMCAL/AliEMCALPID.cxx b/EMCAL/AliEMCALPID.cxx
index 60be8177826..71cd6631c4a 100644
--- a/EMCAL/AliEMCALPID.cxx
+++ b/EMCAL/AliEMCALPID.cxx
@@ -14,64 +14,33 @@
  **************************************************************************/
 
 /* $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
+
+//   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::kSPECIESN]  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]
@@ -83,109 +52,87 @@
 //	              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 <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 "AliESDEvent.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)
+  AliEMCALPID::AliEMCALPID()
+	: AliEMCALPIDUtils(), 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)
+: AliEMCALPIDUtils(), 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->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));
@@ -208,149 +155,70 @@ void AliEMCALPID::RunPID(AliESDEvent *esd)
 	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::kSPECIESN+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]);
+  if(fReconstructor){
     
-
-    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(!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;
-}