-/**************************************************************************\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
-AliEMCALPIDResponse::AliEMCALPIDResponse():\r
-fNorm(0)\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
- 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
-\r
-AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)\r
-{\r
- //\r
- // The assignment operator\r
- //\r
- fNorm = other.fNorm;\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
-\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 kFALSE //
+//////////////////////////////////////////////////////////////////////////
+
+#include <TF1.h>
+#include <TMath.h>
+
+#include "AliEMCALPIDResponse.h" //class header
+
+#include "AliLog.h"
+
+ClassImp(AliEMCALPIDResponse)
+
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+AliEMCALPIDResponse::AliEMCALPIDResponse():
+ TObject(),
+ fNorm(NULL),
+ fCurrCentrality(-1.),
+ fkPIDParams(NULL)
+{
+ //
+ // The default constructor
+ //
+
+
+ fNorm = new TF1("fNorm","gaus",-20,20);
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
+ TObject(other),
+ fNorm(other.fNorm),
+ fCurrCentrality(other.fCurrCentrality),
+ 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;
+ fCurrCentrality = other.fCurrCentrality;
+ 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, charge);
+
+ // 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, charge);
+
+ // 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;
+
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
+ //
+ //
+ Double_t fRange = 5.0; // hardcoded (???)
+ Double_t nsigma = 0.0;
+
+
+ // Check the charge
+ if( charge != -1 && charge != 1){
+ return kFALSE;
+ }
+
+
+ // default value (will be returned, if pt below threshold)
+ for (Int_t species = 0; species < nSpecies; species++) {
+ pEMCAL[species] = 1./nSpecies;
+ }
+
+ // set E/p range
+ if(eop < 0.05) eop = 0.05;
+ if(eop > 2.00) eop = 2.00;
+
+ for (Int_t species = 0; species < nSpecies; species++) {
+
+ AliPID::EParticleType type = AliPID::EParticleType(species);
+
+ // Get the parameters for this particle type and pt
+ const TVectorD *params = GetParams(species, pt, charge);
+
+ // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
+ if(!params) return kFALSE;
+
+ 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 kTRUE;
+
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) 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)
+ //
+ // for PbPb the parametrization is done centrality dependent (marked by TString "Centrality")
+ // so first the correct centrality bin has to be found
+
+ // **** Centrality bins (hard coded for the moment)
+ const Int_t nCent = 7;
+ Int_t centBins[nCent+1] = {0,10,20,30,40,50,70,90};
+
+ if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
+ if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
+
+ TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
+ if(!particlePar) return NULL;
+
+ const TVectorD *parameters = NULL;
+ Double_t momMin = 0.;
+ Double_t momMax = 0.;
+
+ // is the centrality dependent parametrization used
+ TString arrayName = particlePar->GetName();
+
+ // centrality dependent parametrization
+ if(arrayName.Contains("Centrality")){
+
+ for(Int_t iCent = 0; iCent < nCent; iCent++){
+
+ if( fCurrCentrality > centBins[iCent] && fCurrCentrality < centBins[iCent+1] ){
+
+ TObjArray * centPar = dynamic_cast<TObjArray *>(particlePar->At(iCent));
+ if(!centPar) return NULL;
+
+ TIter centIter(centPar);
+ parameters = NULL;
+ momMin = 0.;
+ momMax = 0.;
+
+ while((parameters = static_cast<const TVectorD *>(centIter()))){
+
+ momMin = (*parameters)[0];
+ momMax = (*parameters)[1];
+
+ if( fPt > momMin && fPt < momMax ) return parameters;
+
+ }
+ }
+ }
+ }
+
+ // NO centrality dependent parametrization
+ else{
+
+ TIter parIter(particlePar);
+ 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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff