EMCAL Response in OADB
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 19 Oct 2011 07:34:04 +0000 (07:34 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 19 Oct 2011 07:34:04 +0000 (07:34 +0000)
M. Weber

STEER/STEERBase/AliEMCALPIDResponse.cxx
STEER/STEERBase/AliEMCALPIDResponse.h
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliPIDResponse.h

index c11cb8d..1f8af9d 100644 (file)
-/**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- *                                                                        *\r
- * Author: The ALICE Off-line Project.                                    *\r
- * Contributors are mentioned in the code where appropriate.              *\r
- *                                                                        *\r
- * Permission to use, copy, modify and distribute this software and its   *\r
- * documentation strictly for non-commercial purposes is hereby granted   *\r
- * without fee, provided that the above copyright notice appears in all   *\r
- * copies and that both the copyright notice and this permission notice   *\r
- * appear in the supporting documentation. The authors make no claims     *\r
- * about the suitability of this software for any purpose. It is          *\r
- * provided "as is" without express or implied warranty.                  *\r
- **************************************************************************/\r
-//////////////////////////////////////////////////////////////////////////\r
-//                                                                      //\r
-// AliEMCALPIDResponse                                                  //\r
-//                                                                      //\r
-// EMCAL class to perfom PID                                            //\r
-// This is a prototype and still under development                      //\r
-//                                                                      //\r
-// ---------------------------------------------------------------------//\r
-// GetNumberOfSigmas():                                                 //\r
-//                                                                      //\r
-// Electrons:  Number of Sigmas for E/p value                           //\r
-//             Parametrization of LHC11a (after recalibration)          //   \r
-//                                                                      //\r
-// NON electrons:                                                       //\r
-//             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //\r
-//             --> return +/- 999                                       //\r
-//             Otherwise                                                //\r
-//             --> return nsigma (parametrization of LHC10e)            //   \r
-//                                                                      //\r
-// ---------------------------------------------------------------------//\r
-// ComputeEMCALProbability():                                           //\r
-//                                                                      //\r
-// Electrons:  Probability from Gaussian distribution                    //\r
-//                                                                      //\r
-// NON electrons:                                                       //\r
-//             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //\r
-//             --> probability to find particles below or above thr.    //\r
-//             Otherwise                                                //\r
-//             -->  Probability from Gaussian distribution              //  \r
-//                  (proper normalization to each other?)               //\r
-//                                                                      //\r
-//////////////////////////////////////////////////////////////////////////\r
-\r
-#include <TF1.h>\r
-#include <TMath.h>\r
-\r
-#include "AliEMCALPIDResponse.h"       //class header\r
-\r
-#include "AliLog.h"   \r
-\r
-ClassImp(AliEMCALPIDResponse)\r
-\r
-const Float_t AliEMCALPIDResponse::fLowEoP  = 0.5;   // lower E/p threshold for NON electrons\r
-const Float_t AliEMCALPIDResponse::fHighEoP = 1.5;   // upper E/p threshold for NON electrons\r
-\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-AliEMCALPIDResponse::AliEMCALPIDResponse():\r
-  TObject(),\r
-  fNorm(NULL)\r
-{\r
-  //\r
-  //  The default constructor\r
-  //\r
-\r
-  for(Int_t i = 0; i < fNptBins; i++){\r
-\r
-    fPtCutMin[i] = 0.0;\r
-\r
-    for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++){\r
-\r
-      fMeanEoP[j][i]  = 0.0;\r
-      fSigmaEoP[j][i] = 0.0;\r
-      fProbLow[j][i]  = 0.0;\r
-      fProbHigh[j][i] = 0.0;\r
-\r
-    }\r
-  } \r
-  fPtCutMin[fNptBins] = 0.0;\r
-\r
-  fNorm = new TF1("fNorm","gaus",-20,20); \r
-\r
-  SetPtBoundary();\r
-  SetParametrizations();\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):\r
-  TObject(other),\r
-  fNorm(other.fNorm)\r
-{\r
-  //\r
-  //  The copy constructor\r
-  //\r
-    for(Int_t i = 0; i < fNptBins; i++)\r
-    {\r
-       fPtCutMin[i] = 0.0;\r
-       for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++)\r
-       {\r
-           fMeanEoP[j][i]  = 0.0;\r
-           fSigmaEoP[j][i] = 0.0;\r
-           fProbLow[j][i]  = 0.0;\r
-           fProbHigh[j][i] = 0.0;\r
-       }\r
-    }\r
-  \r
-    fPtCutMin[fNptBins] = 0.0;\r
-    SetPtBoundary();\r
-    SetParametrizations();\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)\r
-{\r
-  //\r
-  //  The assignment operator\r
-  //\r
-\r
-  if(this == &other) return *this;\r
-  \r
-  // Make copy\r
-  TObject::operator=(other);\r
-  fNorm = other.fNorm;\r
-\r
-  for(Int_t i = 0; i < fNptBins; i++)\r
-    {\r
-       fPtCutMin[i] = 0.0;\r
-       for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++)\r
-       {\r
-           fMeanEoP[j][i]  = 0.0;\r
-           fSigmaEoP[j][i] = 0.0;\r
-           fProbLow[j][i]  = 0.0;\r
-           fProbHigh[j][i] = 0.0;\r
-       }\r
-    }\r
-  \r
-  fPtCutMin[fNptBins] = 0.0;\r
-  SetPtBoundary();\r
-  SetParametrizations();\r
-\r
-  return *this;\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-AliEMCALPIDResponse::~AliEMCALPIDResponse() {\r
-\r
-  delete fNorm;\r
-\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-void AliEMCALPIDResponse::SetPtBoundary(){\r
-  //\r
-  // Set boundaries for momentum bins\r
-  //\r
-  fPtCutMin[0] = 1.5;\r
-  fPtCutMin[1] = 2.5;\r
-  fPtCutMin[2] = 3.5;\r
-  fPtCutMin[3] = 4.5;\r
-  fPtCutMin[4] = 5.5;\r
-  fPtCutMin[5] = 6.5;\r
-\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-void AliEMCALPIDResponse::SetParametrizations(){\r
-\r
-// This are the preliminary parametrizations (hard coded)\r
-// For electrons from LHC11a (Deepa Thomas)\r
-// For NON-electrons from LHC10e (TOF/TPC analysis)\r
-  \r
-  // Gaussian mean\r
-  Float_t mean[4][6]    = {\r
-    { 0.932, 0.997, 0.998, 1.001, 1.011, 1.011 },                              // electrons\r
-    { 0.227804, 0.34839, 0.404077, -0.107795, -4.14584, 0.5 },             // NON electrons\r
-    { -2.10377, 0.0582898, 0.0582898, 0.0582898, 0.0582898, 0.0582898  },   // protons\r
-    { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5}                                              // anti-protons\r
-  };        \r
-  \r
-  // Gaussian sigma\r
-  Float_t sigma[4][6]= {\r
-    { 0.0866, 0.0693, 0.0664, 0.0583, 0.0488, 0.0515  },                      // electrons\r
-    { 0.310831, 0.267586, 0.404077, 0.381968, 1.46183, 0.314687 },          // NON electrons\r
-    { 0.603209, 0.255332, 0.255332, 0.255332, 0.255332, 0.255332},          // protons\r
-    { 0.516837, 0.351516,0.351516,0.351516,0.351516,0.351516 }              // anti - protons\r
-  };\r
-\r
-  // lower probability\r
-  Float_t probL[3][6] = {\r
-    { 0.928689, 0.938455, 0.940448, 0.948496, 0.955981, 0.951923 },         // NON electrons\r
-    { 0.974518, 0.978088, 0.975089, 0.975089, 0.975089,0.975089},           // protons\r
-    { 0.824037, 0.861149, 0.898734, 0.898734, 0.898734, 0.898734},          // anti - protons\r
-  };\r
-\r
-  // upper probability\r
-  Float_t probH[3][6] = {\r
-    { 0.00030227, 4.04106e-05, 0.000147406, 0., 0.000956938, 0.00106838 },  // NON electrons\r
-    { 0.000157945, 0., 0., 0., 0., 0. },                                    // protons\r
-    { 0.00343237, 0., 0., 0., 0., 0.}                                       // anti - protons\r
-  };\r
-  \r
-  \r
-  // set parametrizations\r
-  Int_t spec = 0;\r
-  for (Int_t species = 0; species < 2*AliPID::kSPECIES; species++) {          // first negative particles and then positive\r
-    for (Int_t pt = 0; pt < fNptBins; pt++){\r
-\r
-      switch(species){\r
-      case 0:      // electrons\r
-       spec = 0;\r
-       break;\r
-      case 4:      // anti - protons\r
-       spec = 3;\r
-       break;\r
-      case 5:      // positrons\r
-       spec = 0;\r
-       break;\r
-      case 9:      // protons\r
-       spec = 2;\r
-       break;\r
-      default:     // NON electrons\r
-       spec = 1;\r
-       break;\r
-      }\r
-    \r
-\r
-      fMeanEoP[species][pt]  = mean[spec][pt];    \r
-      fSigmaEoP[species][pt] = sigma[spec][pt];  \r
-      if( spec == 0) {     // electrons have NO lower and upper probability thresholds --> set to 0\r
-       fProbLow[species][pt]  = 0.;\r
-       fProbHigh[species][pt] = 0.;    \r
-      }\r
-      else{\r
-       fProbLow[species][pt]  = probL[spec-1][pt];\r
-       fProbHigh[species][pt] = probH[spec-1][pt];     \r
-      } \r
-    \r
-    }//loop pt bins\r
-  }//loop species\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-Int_t AliEMCALPIDResponse::GetPtBin(Float_t pt) const {\r
-  //\r
-  // Returns the momentum bin index\r
-  //\r
-\r
-  Int_t i = -1;\r
-  while(pt > fPtCutMin[i+1] && i+1 < fNptBins) i++;\r
-\r
-  return i;\r
-}\r
-\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-Double_t AliEMCALPIDResponse::GetExpectedSignal( Float_t pt, AliPID::EParticleType n, Int_t charge) const  {\r
-  //\r
-  // Calculates the expected PID signal as the function of \r
-  // the information stored in the track, for the specified particle type \r
-  //  \r
-\r
-  Double_t signal = 0.;\r
-\r
-  // Check the charge\r
-  if( charge != -1 && charge != 1){\r
-    \r
-    return signal;\r
-  }\r
-  \r
-  // Get the pt bin\r
-  Int_t ptBin = GetPtBin(pt);\r
-  \r
-  // Get the species (first negative , then positive)\r
-  Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
-\r
-  // Get the signal\r
-  if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
-    signal = fMeanEoP[species][ptBin];\r
-  }\r
-\r
-  return signal;\r
-\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-Double_t AliEMCALPIDResponse::GetExpectedSigma( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {\r
-  //\r
-  // Calculates the expected sigma of the PID signal as the function of \r
-  // the information stored in the track, for the specified particle type \r
-  //  \r
-  //\r
-  \r
-  Double_t sigma = 999.;\r
-\r
-  // Check the charge\r
-  if( charge != -1 && charge != 1){\r
-    \r
-    return sigma;\r
-  }  \r
-\r
-  // Get the pt bin\r
-  Int_t ptBin = GetPtBin(pt);\r
-  \r
-  // Get the species (first negative , then positive)\r
-  Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
-\r
-  // Get the sigma\r
-  if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
-    sigma = fSigmaEoP[species][ptBin];\r
-  }\r
-\r
-  return sigma;\r
\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {\r
-  //\r
-  // Calculates the expected sigma of the PID signal as the function of \r
-  // the information stored in the track, for the specified particle type \r
-  //  \r
-  //\r
-  \r
-  Double_t norm = 1.;\r
-\r
-  // Check the charge\r
-  if( charge != -1 && charge != 1){\r
-    \r
-    return norm;\r
-  }\r
-\r
-  // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )\r
-  fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*GetExpectedSigma(pt,n,charge)*GetExpectedSigma(pt,n,charge)),GetExpectedSignal(pt,n,charge),GetExpectedSigma(pt,n,charge));\r
-  norm = 1./fNorm->Integral(fLowEoP,fHighEoP)*(1-GetLowProb(pt,n,charge)-GetHighProb(pt,n,charge));\r
-\r
-  return norm;\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {\r
-      \r
-  Double_t mean  = GetExpectedSignal(pt,n,charge);\r
-  Double_t sigma = GetExpectedSigma(pt,n,charge);\r
-\r
-  // if electron\r
-  if(n == AliPID::kElectron){\r
-    return (eop - mean) / sigma;\r
-  }\r
-\r
-  // if NON electron\r
-  else{\r
-    if ( eop < fLowEoP )\r
-      return -999.;    // not parametrized \r
-    else if ( eop > fHighEoP )\r
-      return 999.;     // not parametrized \r
-    else{\r
-      return (eop - mean) / sigma; \r
-    }\r
-  }\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-Double_t AliEMCALPIDResponse::GetLowProb( Float_t pt, AliPID::EParticleType n,  Int_t charge) const {\r
-  //\r
-  //\r
-  \r
-  Double_t prob = 0.;\r
-\r
-  // Check the charge\r
-  if( charge != -1 && charge != 1){\r
-    \r
-    return prob;\r
-  }\r
-  \r
-  // Get the pt bin\r
-  Int_t ptBin = GetPtBin(pt);\r
-  \r
-  // Get the species (first negative , then positive)\r
-  Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
-\r
-  // Get the probability\r
-  if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
-    prob = fProbLow[species][ptBin];\r
-  }\r
-\r
-  return prob;\r
\r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-Double_t AliEMCALPIDResponse::GetHighProb( Float_t pt, AliPID::EParticleType n,  Int_t charge) const {\r
-  //\r
-  //\r
-  \r
-  Double_t prob = 0.;\r
-\r
-  // Check the charge\r
-  if( charge != -1 && charge != 1){\r
-    \r
-    return prob;\r
-  }\r
-  \r
-  // Get the pt bin\r
-  Int_t ptBin = GetPtBin(pt);\r
-  \r
-  // Get the species (first negative , then positive)\r
-  Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;\r
-\r
-  // Get the probability\r
-  if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){\r
-    prob = fProbHigh[species][ptBin];\r
-  }\r
-\r
-  return prob;\r
-  \r
-}\r
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\r
-Double_t AliEMCALPIDResponse::ComputeEMCALProbability(Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {\r
-  //\r
-  //\r
-  Double_t fRange  = 5.0;   // hardcoded \r
-  Double_t nsigma  = 0.0;\r
-  Double_t sigma   = 0.0;\r
-  Double_t proba   = 999.;\r
-  \r
-\r
-  // Check the charge\r
-  if( charge != -1 && charge != 1){\r
-    \r
-    return proba;\r
-  }\r
-\r
\r
-  // default value (will be returned, if pt below threshold)\r
-  for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
-    pEMCAL[species] = 999.;\r
-  }\r
-\r
-  if( GetPtBin(pt) > -1 ){\r
-\r
-    // set E/p range\r
-    if(eop < 0.05) eop = 0.05;\r
-    if(eop > 2.00) eop = 2.00;\r
-\r
-    for (Int_t species = 0; species < AliPID::kSPECIES; species++) {\r
-\r
-      AliPID::EParticleType type = AliPID::EParticleType(species);\r
-\r
-      // get nsigma value for each particle type at this E/p value\r
-      nsigma = GetNumberOfSigmas(pt,eop,type,charge);\r
-      sigma  = GetExpectedSigma(pt,type,charge);\r
-\r
-      // electrons (standard Gaussian calculation of probabilities)\r
-      if(type == AliPID::kElectron){\r
-       if (TMath::Abs(nsigma) > fRange) {\r
-         pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
-       }\r
-       else{\r
-         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
-       }\r
-      }\r
-      //NON electrons\r
-      else{\r
-       // E/p < 0.5  -->  return probability below E/p = 0.5\r
-       if ( nsigma == -999){\r
-         pEMCAL[species] = GetLowProb(pt,type,charge);\r
-       }\r
-       // E/p > 1.5  -->  return probability above E/p = 1.5\r
-       else if ( nsigma == 999){\r
-         pEMCAL[species] = GetHighProb(pt,type,charge);\r
-       }\r
-       // in parametrized region --> calculate probability for corresponding Gauss curve\r
-       else{\r
-         pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);\r
-       \r
-         // normalize to total probability == 1\r
-         pEMCAL[species]*=GetExpectedNorm(pt,type,charge);\r
-       }\r
-      }\r
-    }\r
-\r
-    // return the electron probability\r
-    proba = pEMCAL[AliPID::kElectron];  \r
-\r
-  }\r
-\r
-  return proba;\r
-\r
-}\r
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+//////////////////////////////////////////////////////////////////////////
+//                                                                      //
+// AliEMCALPIDResponse                                                  //
+//                                                                      //
+// EMCAL class to perfom PID                                            //
+// This is a prototype and still under development                      //
+//                                                                      //
+// ---------------------------------------------------------------------//
+// GetNumberOfSigmas():                                                 //
+//                                                                      //
+// Electrons:  Number of Sigmas for E/p value                           //
+//             Parametrization of LHC11a (after recalibration)          //
+//                                                                      //
+// NON electrons:                                                       //
+//             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //
+//             --> return +/- 99                                        //
+//             Otherwise                                                //
+//             --> return nsigma (parametrization of LHC10e)            //
+//                                                                      //
+// NO Parametrization (outside pT range): --> return -999               //
+//                                                                      //
+// ---------------------------------------------------------------------//
+// ComputeEMCALProbability():                                           //
+//                                                                      //
+// Electrons:  Probability from Gaussian distribution                   //
+//                                                                      //
+// NON electrons:                                                       //
+//             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //
+//             --> probability to find particles below or above thr.    //
+//             Otherwise                                                //
+//             -->  Probability from Gaussian distribution              //
+//                  (proper normalization to each other?)               //
+//                                                                      //
+// NO Parametrization (outside pT range): --> return 999                //
+//////////////////////////////////////////////////////////////////////////
+
+#include <TF1.h>
+#include <TMath.h>
+
+#include "AliEMCALPIDResponse.h"       //class header
+
+#include "AliLog.h"   
+
+ClassImp(AliEMCALPIDResponse)
+
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+AliEMCALPIDResponse::AliEMCALPIDResponse():
+  TObject(),
+  fNorm(NULL),
+  fkPIDParams(NULL)
+{
+  //
+  //  The default constructor
+  //
+
+
+  fNorm = new TF1("fNorm","gaus",-20,20); 
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
+  TObject(other),
+  fNorm(other.fNorm),
+  fkPIDParams(other.fkPIDParams)
+{
+  //
+  //  The copy constructor
+  //
+
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
+{
+  //
+  //  The assignment operator
+  //
+
+  if(this == &other) return *this;
+  
+  // Make copy
+  TObject::operator=(other);
+  fNorm = other.fNorm;
+  fkPIDParams = other.fkPIDParams;
+
+  return *this;
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+AliEMCALPIDResponse::~AliEMCALPIDResponse() {
+
+  delete fNorm;
+
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {
+  //
+  // Calculates the expected sigma of the PID signal as the function of 
+  // the information stored in the track, for the specified particle type 
+  //  
+  //
+  
+  Double_t norm = 1.;
+
+  // Check the charge
+  if( charge != -1 && charge != 1){
+    return norm;
+  }
+
+  // Get the parameters for this particle type and pt
+  const TVectorD *params = GetParams(n, pt);
+
+  // IF not in momentum range, NULL is returned --> return default value
+  if(!params) return norm;
+
+  Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
+  Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
+  Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
+  Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
+  Double_t probLow  = (*params)[6];   // probability to be below eopMin
+  Double_t probHigh = (*params)[7];   // probability to be above eopMax
+
+  // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
+  fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*sigma*sigma),mean,sigma);
+  norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);
+
+  return norm;
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {
+      
+  Double_t nsigma = -999.;
+
+  // Check the charge
+  if( charge != -1 && charge != 1){
+    return nsigma;
+  }
+
+  // Get the parameters for this particle type and pt
+  const TVectorD *params = GetParams(n, pt);
+
+  // IF not in momentum range, NULL is returned --> return default value
+  if(!params) return nsigma;
+
+  Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
+  Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
+  Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
+  Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
+
+  // if electron
+  if(n == AliPID::kElectron){
+    if(sigma != 0) nsigma = (eop - mean) / sigma;
+  }
+
+  // if NON electron
+  else{
+    if ( eop < eopMin )
+      nsigma = -99;    // not parametrized (smaller than eopMin)
+    else if ( eop > eopMax )
+      nsigma = 99.;     // not parametrized (bigger than eopMax)
+    else{
+      if(sigma != 0) nsigma = (eop - mean) / sigma; 
+    }
+  }
+
+  return nsigma;
+
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliEMCALPIDResponse::ComputeEMCALProbability(Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
+  //
+  //
+  Double_t fRange  = 5.0;   // hardcoded (???)
+  Double_t nsigma  = 0.0;
+  Double_t proba   = 999.;
+  
+
+  // Check the charge
+  if( charge != -1 && charge != 1){
+    return proba;
+  }
+
+  // default value (will be returned, if pt below threshold)
+  for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
+    pEMCAL[species] = 999.;
+  }
+
+  // set E/p range
+  if(eop < 0.05) eop = 0.05;
+  if(eop > 2.00) eop = 2.00;
+  
+  for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
+    
+    AliPID::EParticleType type = AliPID::EParticleType(species);
+
+    // Get the parameters for this particle type and pt
+    const TVectorD *params = GetParams(species, pt);
+    
+    // IF not in momentum range, NULL is returned --> return default value
+    if(!params) return proba;
+
+    Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
+    Double_t probLow  = (*params)[6];   // probability to be below eopMin
+    Double_t probHigh = (*params)[7];   // probability to be above eopMax
+
+    // get nsigma value for each particle type at this E/p value
+    nsigma = GetNumberOfSigmas(pt,eop,type,charge);
+
+    // electrons (standard Gaussian calculation of probabilities)
+    if(type == AliPID::kElectron){
+      if (TMath::Abs(nsigma) > fRange) {
+       pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
+      }
+      else{
+       pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
+      }
+    }
+    //NON electrons
+    else{
+      // E/p < eopMin  -->  return probability below E/p = eopMin
+      if ( nsigma == -99){
+       pEMCAL[species] = probLow;
+      }
+      // E/p > eopMax  -->  return probability above E/p = eopMax
+      else if ( nsigma == 99){
+       pEMCAL[species] = probHigh;
+      }
+      // in parametrized region --> calculate probability for corresponding Gauss curve
+      else{
+       pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
+       
+       // normalize to total probability == 1
+       pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
+      }
+    }
+    
+    // return the electron probability
+    proba = pEMCAL[AliPID::kElectron];  
+
+  }
+
+  return proba;
+
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt) const {
+  //
+  // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
+  //
+  // 0 = momMin
+  // 1 = momMax
+  // 2 = mean of Gaus
+  // 3 = sigma of Gaus
+  // 4 = eopLow   
+  // 5 = eopHig   
+  // 6 = probLow  (not used for electrons)
+  // 7 = probHigh (not used for electrons)
+  //
+
+  if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
+
+  TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
+  if(!particlePar) return NULL;
+  
+  TIter parIter(particlePar);
+  const TVectorD *parameters = NULL;
+  Double_t momMin = 0.;
+  Double_t momMax = 0.;
+
+  while((parameters = static_cast<const TVectorD *>(parIter()))){
+
+    momMin = (*parameters)[0];
+    momMax = (*parameters)[1];
+
+    if( fPt > momMin && fPt < momMax ) return parameters;
+
+  }  
+  AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));
+
+  return parameters;
+}
index 5a8af3e..9425c11 100644 (file)
@@ -1,72 +1,57 @@
-#ifndef AliEMCALPIDResponse_h\r
-#define AliEMCALPIDResponse_h\r
-\r
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice                               */\r
-\r
-//////////////////////////////////////////////////////////////////////////\r
-//                                                                      //\r
-// AliEMCALPIDResponse                                                  //\r
-//                                                                      //\r
-// EMCAL class to perfom PID                                            //\r
-// This is a prototype and still under development                      //\r
-//                                                                      //\r
-//                                                                      //\r
-//////////////////////////////////////////////////////////////////////////\r
-\r
-#include "AliPID.h"\r
-class TF1;\r
-\r
-class AliEMCALPIDResponse: public TObject \r
-{\r
-public : \r
-    AliEMCALPIDResponse();    //ctor\r
-    AliEMCALPIDResponse( const AliEMCALPIDResponse& other);                //copy ructor\r
-    AliEMCALPIDResponse &operator=( const AliEMCALPIDResponse& other);     //assignment operator\r
-\r
-    virtual ~AliEMCALPIDResponse();     //dtor\r
-  \r
-\r
-    // Getters\r
-    Int_t     GetPtBin(Float_t pt) const;\r
-    Int_t     GetNPtBins() const {return fNptBins;};\r
-\r
-    Float_t   GetLowEoP() const {return fLowEoP;};\r
-    Float_t   GetHighEoP() const {return fHighEoP;};\r
-\r
-    Double_t  GetExpectedSignal( Float_t pt, AliPID::EParticleType n,  Int_t charge) const;\r
-    Double_t  GetExpectedSigma ( Float_t pt, AliPID::EParticleType n,  Int_t charge) const;\r
-    Double_t  GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const;\r
-    Double_t  GetExpectedNorm  ( Float_t pt, AliPID::EParticleType n,  Int_t charge) const;\r
-    Double_t  GetLowProb       ( Float_t pt, AliPID::EParticleType n,  Int_t charge) const;\r
-    Double_t  GetHighProb      ( Float_t pt, AliPID::EParticleType n,  Int_t charge) const;\r
-\r
-    //Setters\r
-    void   SetPtBoundary();\r
-    void   SetParametrizations();\r
-\r
-    // EMCAL probability -> should go to another place?\r
-    Double_t ComputeEMCALProbability( Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const;\r
-\r
-protected:\r
-  \r
-private:\r
-\r
-  TF1 *fNorm;                            // Gauss function for normalizing NON electron probabilities \r
-\r
-  static const Int_t fNptBins   = 6;     // number of momentum bins\r
-  static const Float_t fLowEoP;   // lower E/p threshold for NON electrons\r
-  static const Float_t fHighEoP;   // upper E/p threshold for NON electrons\r
-\r
-  Float_t fPtCutMin[fNptBins+1];                       // min values for pt bins\r
-  Float_t fMeanEoP[2*AliPID::kSPECIES][fNptBins];      // mean value of E/p distribution (charge dependent)\r
-  Float_t fSigmaEoP[2*AliPID::kSPECIES][fNptBins];     // mean value of E/p distribution (charge dependent)\r
-  Float_t fProbLow[2*AliPID::kSPECIES][fNptBins];      // probability below E/p threshold for NON electrons (charge dependent)\r
-  Float_t fProbHigh[2*AliPID::kSPECIES][fNptBins];     // probability above E/p threshold for NON electrons (charge dependent)\r
-\r
-\r
-  ClassDef(AliEMCALPIDResponse, 1)\r
-};\r
-\r
-#endif // #ifdef AliEMCALPIDResponse_cxx\r
-\r
+#ifndef AliEMCALPIDResponse_h
+#define AliEMCALPIDResponse_h
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//////////////////////////////////////////////////////////////////////////
+//                                                                      //
+// AliEMCALPIDResponse                                                  //
+//                                                                      //
+// EMCAL class to perfom PID                                            //
+// This is a prototype and still under development                      //
+//                                                                      //
+// Author: Michael Weber (m.weber@cern.ch)                              //
+//////////////////////////////////////////////////////////////////////////
+
+#include "AliPID.h"
+#include <TVectorD.h>
+
+class TF1;
+
+class AliEMCALPIDResponse: public TObject 
+{
+public : 
+    AliEMCALPIDResponse();    //ctor
+    AliEMCALPIDResponse( const AliEMCALPIDResponse& other);                //copy ructor
+    AliEMCALPIDResponse &operator=( const AliEMCALPIDResponse& other);     //assignment operator
+
+    virtual ~AliEMCALPIDResponse();     //dtor
+  
+
+    // Getters
+    Double_t  GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const;
+    Double_t  GetExpectedNorm  ( Float_t pt, AliPID::EParticleType n,  Int_t charge) const;
+  
+    //Setters
+    void   SetPIDParams(const TObjArray * params) { fkPIDParams = params; }
+    
+
+    // EMCAL probability -> should go to another place?
+    Double_t ComputeEMCALProbability( Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const;
+
+protected:
+  
+private:
+
+  TF1 *fNorm;                            // Gauss function for normalizing NON electron probabilities 
+
+  const TObjArray *fkPIDParams;               // PID Params
+
+  const TVectorD* GetParams(Int_t nParticle, Float_t fPt) const; 
+
+  ClassDef(AliEMCALPIDResponse, 1)
+};
+
+#endif // #ifdef AliEMCALPIDResponse_cxx
+
index 1d2e7e4..f1ebda2 100644 (file)
@@ -35,6 +35,7 @@
 #include <AliLog.h>
 #include <AliPID.h>
 #include <AliOADBContainer.h>
+#include <AliTRDPIDParams.h>
 #include <AliTRDPIDReference.h>
 
 #include "AliPIDResponse.h"
@@ -67,6 +68,7 @@ fTRDPIDParams(0x0),
 fTRDPIDReference(0x0),
 fTOFTimeZeroType(kBest_T0),
 fTOFres(100.),
+fEMCALPIDParams(0x0),
 fCurrentEvent(0x0)
 {
   //
@@ -117,6 +119,7 @@ fTRDPIDParams(0x0),
 fTRDPIDReference(0x0),
 fTOFTimeZeroType(AliPIDResponse::kBest_T0),
 fTOFres(100.),
+fEMCALPIDParams(0x0),
 fCurrentEvent(0x0)
 {
   //
@@ -156,6 +159,7 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fResolutionCorrection=0x0;
     fTRDPIDParams=0x0;
     fTRDPIDReference=0x0;
+    fEMCALPIDParams=0x0;
     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
     fTOFTimeZeroType=AliPIDResponse::kBest_T0;
     fTOFres=100.;
@@ -537,6 +541,9 @@ void AliPIDResponse::ExecNewRun()
 
   SetTRDPidResponseMaster(); 
   InitializeTRDResponse();
+
+  SetEMCALPidResponseMaster(); 
+  InitializeEMCALResponse();
   
   fTOFResponse.SetTimeResolution(fTOFres);
 }
@@ -744,7 +751,7 @@ void AliPIDResponse::SetTRDPidResponseMaster()
     AliError("Failed initializing PID Params from OADB");
   } else {
     AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
-    fTRDPIDParams = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
+    fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
     if(!fTRDPIDParams){
       AliError(Form("TRD Params not found in run %d", fRun));
     }
@@ -798,3 +805,49 @@ Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t
   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
 }
 
+//______________________________________________________________________________
+void AliPIDResponse::SetEMCALPidResponseMaster()
+{
+  //
+  // Load the EMCAL pid response functions from the OADB
+  //
+  TObjArray* fEMCALPIDParamsRun      = NULL;
+  TObjArray* fEMCALPIDParamsPass     = NULL;
+
+  if(fEMCALPIDParams) return;
+  AliOADBContainer contParams("contParams"); 
+
+  Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
+  if(statusPars){
+    AliError("Failed initializing PID Params from OADB");
+  } 
+  else {
+    AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
+
+    fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
+    if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
+    if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
+
+    if(!fEMCALPIDParams){
+      AliError(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
+      AliInfo("Will take the standard LHC11a pass2 instead ...");
+
+      fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(144871));
+      if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",2)));
+      if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
+
+      if(!fEMCALPIDParams){
+       AliError(Form("DEFAULT EMCAL Params (LHC11a pass2) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));       
+      }
+    }
+  }
+}
+
+//______________________________________________________________________________
+void AliPIDResponse::InitializeEMCALResponse(){
+  //
+  // Set PID Params to the EMCAL PID response
+  // 
+  fEMCALResponse.SetPIDParams(fEMCALPIDParams);
+
+}
index 17a265e..18fad85 100644 (file)
@@ -127,6 +127,8 @@ private:
   Int_t   fTOFTimeZeroType;            //! default start time type for tof (ESD)
   Float_t fTOFres;                     //! TOF resolution
 
+  TObjArray *fEMCALPIDParams;             //! EMCAL PID Params
+
   AliVEvent *fCurrentEvent;            //! event currently being processed
   
   void ExecNewRun();
@@ -148,7 +150,11 @@ private:
   void InitializeTRDResponse();
 
   //TOF
-  
+
+  //EMCAL
+  void SetEMCALPidResponseMaster();
+  void InitializeEMCALResponse();
+
   //
   void SetRecoInfo();