1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\r
7 * Permission to use, copy, modify and distribute this software and its *
\r
8 * documentation strictly for non-commercial purposes is hereby granted *
\r
9 * without fee, provided that the above copyright notice appears in all *
\r
10 * copies and that both the copyright notice and this permission notice *
\r
11 * appear in the supporting documentation. The authors make no claims *
\r
12 * about the suitability of this software for any purpose. It is *
\r
13 * provided "as is" without express or implied warranty. *
\r
14 **************************************************************************/
\r
15 //////////////////////////////////////////////////////////////////////////
\r
17 // AliEMCALPIDResponse //
\r
19 // EMCAL class to perfom PID //
\r
20 // This is a prototype and still under development //
\r
22 // ---------------------------------------------------------------------//
\r
23 // GetNumberOfSigmas(): //
\r
25 // Electrons: Number of Sigmas for E/p value //
\r
26 // Parametrization of LHC11a (after recalibration) //
\r
28 // NON electrons: //
\r
29 // Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5) //
\r
30 // --> return +/- 999 //
\r
32 // --> return nsigma (parametrization of LHC10e) //
\r
34 // ---------------------------------------------------------------------//
\r
35 // ComputeEMCALProbability(): //
\r
37 // Electrons: Probability from Gaussian distribution //
\r
39 // NON electrons: //
\r
40 // Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5) //
\r
41 // --> probability to find particles below or above thr. //
\r
43 // --> Probability from Gaussian distribution //
\r
44 // (proper normalization to each other?) //
\r
46 //////////////////////////////////////////////////////////////////////////
\r
51 #include "AliEMCALPIDResponse.h" //class header
\r
53 #include "AliLog.h"
\r
55 ClassImp(AliEMCALPIDResponse)
\r
56 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
57 AliEMCALPIDResponse::AliEMCALPIDResponse():
\r
61 // The default constructor
\r
64 for(Int_t i = 0; i < fNptBins; i++){
\r
68 for(Int_t j = 0; j < 2*AliPID::kSPECIES; j++){
\r
70 fMeanEoP[j][i] = 0.0;
\r
71 fSigmaEoP[j][i] = 0.0;
\r
72 fProbLow[j][i] = 0.0;
\r
73 fProbHigh[j][i] = 0.0;
\r
77 fPtCutMin[fNptBins] = 0.0;
\r
79 fNorm = new TF1("fNorm","gaus",-20,20);
\r
82 SetParametrizations();
\r
84 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
85 void AliEMCALPIDResponse::SetPtBoundary(){
\r
87 // Set boundaries for momentum bins
\r
97 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
98 void AliEMCALPIDResponse::SetParametrizations(){
\r
100 // This are the preliminary parametrizations (hard coded)
\r
101 // For electrons from LHC11a (Deepa Thomas)
\r
102 // For NON-electrons from LHC10e (TOF/TPC analysis)
\r
105 Float_t mean[4][6] = {
\r
106 { 0.932, 0.997, 0.998, 1.001, 1.011, 1.011 }, // electrons
\r
107 { 0.227804, 0.34839, 0.404077, -0.107795, -4.14584, 0.5 }, // NON electrons
\r
108 { -2.10377, 0.0582898, 0.0582898, 0.0582898, 0.0582898, 0.0582898 }, // protons
\r
109 { 0.5, 0.5, 0.5, 0.5, 0.5, 0.5} // anti-protons
\r
113 Float_t sigma[4][6]= {
\r
114 { 0.0866, 0.0693, 0.0664, 0.0583, 0.0488, 0.0515 }, // electrons
\r
115 { 0.310831, 0.267586, 0.404077, 0.381968, 1.46183, 0.314687 }, // NON electrons
\r
116 { 0.603209, 0.255332, 0.255332, 0.255332, 0.255332, 0.255332}, // protons
\r
117 { 0.516837, 0.351516,0.351516,0.351516,0.351516,0.351516 } // anti - protons
\r
120 // lower probability
\r
121 Float_t probL[3][6] = {
\r
122 { 0.928689, 0.938455, 0.940448, 0.948496, 0.955981, 0.951923 }, // NON electrons
\r
123 { 0.974518, 0.978088, 0.975089, 0.975089, 0.975089,0.975089}, // protons
\r
124 { 0.824037, 0.861149, 0.898734, 0.898734, 0.898734, 0.898734}, // anti - protons
\r
127 // upper probability
\r
128 Float_t probH[3][6] = {
\r
129 { 0.00030227, 4.04106e-05, 0.000147406, 0., 0.000956938, 0.00106838 }, // NON electrons
\r
130 { 0.000157945, 0., 0., 0., 0., 0. }, // protons
\r
131 { 0.00343237, 0., 0., 0., 0., 0.} // anti - protons
\r
135 // set parametrizations
\r
137 for (Int_t species = 0; species < 2*AliPID::kSPECIES; species++) { // first negative particles and then positive
\r
138 for (Int_t pt = 0; pt < fNptBins; pt++){
\r
141 case 0: // electrons
\r
144 case 4: // anti - protons
\r
147 case 5: // positrons
\r
153 default: // NON electrons
\r
159 fMeanEoP[species][pt] = mean[spec][pt];
\r
160 fSigmaEoP[species][pt] = sigma[spec][pt];
\r
161 if( spec == 0) { // electrons have NO lower and upper probability thresholds --> set to 0
\r
162 fProbLow[species][pt] = 0.;
\r
163 fProbHigh[species][pt] = 0.;
\r
166 fProbLow[species][pt] = probL[spec-1][pt];
\r
167 fProbHigh[species][pt] = probH[spec-1][pt];
\r
173 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
174 Int_t AliEMCALPIDResponse::GetPtBin(Float_t pt) const {
\r
176 // Returns the momentum bin index
\r
180 while(pt > fPtCutMin[i+1] && i+1 < fNptBins) i++;
\r
185 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
186 Double_t AliEMCALPIDResponse::GetExpectedSignal( Float_t pt, AliPID::EParticleType n, Int_t charge) const {
\r
188 // Calculates the expected PID signal as the function of
\r
189 // the information stored in the track, for the specified particle type
\r
192 Double_t signal = 0.;
\r
194 // Check the charge
\r
195 if( charge != -1 && charge != 1){
\r
201 Int_t ptBin = GetPtBin(pt);
\r
203 // Get the species (first negative , then positive)
\r
204 Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;
\r
207 if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){
\r
208 signal = fMeanEoP[species][ptBin];
\r
214 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
215 Double_t AliEMCALPIDResponse::GetExpectedSigma( Float_t pt, AliPID::EParticleType n, Int_t charge) const {
\r
217 // Calculates the expected sigma of the PID signal as the function of
\r
218 // the information stored in the track, for the specified particle type
\r
222 Double_t sigma = 999.;
\r
224 // Check the charge
\r
225 if( charge != -1 && charge != 1){
\r
231 Int_t ptBin = GetPtBin(pt);
\r
233 // Get the species (first negative , then positive)
\r
234 Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;
\r
237 if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){
\r
238 sigma = fSigmaEoP[species][ptBin];
\r
244 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
245 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n, Int_t charge) const {
\r
247 // Calculates the expected sigma of the PID signal as the function of
\r
248 // the information stored in the track, for the specified particle type
\r
252 Double_t norm = 1.;
\r
254 // Check the charge
\r
255 if( charge != -1 && charge != 1){
\r
260 // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
\r
261 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
262 norm = 1./fNorm->Integral(fLowEoP,fHighEoP)*(1-GetLowProb(pt,n,charge)-GetHighProb(pt,n,charge));
\r
266 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
267 Double_t AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt, Float_t eop, AliPID::EParticleType n, Int_t charge) const {
\r
269 Double_t mean = GetExpectedSignal(pt,n,charge);
\r
270 Double_t sigma = GetExpectedSigma(pt,n,charge);
\r
273 if(n == AliPID::kElectron){
\r
274 return (eop - mean) / sigma;
\r
279 if ( eop < fLowEoP )
\r
280 return -999.; // not parametrized
\r
281 else if ( eop > fHighEoP )
\r
282 return 999.; // not parametrized
\r
284 return (eop - mean) / sigma;
\r
288 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
289 Double_t AliEMCALPIDResponse::GetLowProb( Float_t pt, AliPID::EParticleType n, Int_t charge) const {
\r
293 Double_t prob = 0.;
\r
295 // Check the charge
\r
296 if( charge != -1 && charge != 1){
\r
302 Int_t ptBin = GetPtBin(pt);
\r
304 // Get the species (first negative , then positive)
\r
305 Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;
\r
307 // Get the probability
\r
308 if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){
\r
309 prob = fProbLow[species][ptBin];
\r
315 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
316 Double_t AliEMCALPIDResponse::GetHighProb( Float_t pt, AliPID::EParticleType n, Int_t charge) const {
\r
320 Double_t prob = 0.;
\r
322 // Check the charge
\r
323 if( charge != -1 && charge != 1){
\r
329 Int_t ptBin = GetPtBin(pt);
\r
331 // Get the species (first negative , then positive)
\r
332 Int_t species = n + AliPID::kSPECIES * ( charge + 1 ) / 2;
\r
334 // Get the probability
\r
335 if(species > -1 && species < 2*AliPID::kSPECIES && ptBin > -1 ){
\r
336 prob = fProbHigh[species][ptBin];
\r
342 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\r
343 Double_t AliEMCALPIDResponse::ComputeEMCALProbability(Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
\r
346 Double_t fRange = 5.0; // hardcoded
\r
347 Double_t nsigma = 0.0;
\r
348 Double_t sigma = 0.0;
\r
349 Double_t proba = -1.;
\r
352 // Check the charge
\r
353 if( charge != -1 && charge != 1){
\r
359 // default value (will be returned, if pt below threshold)
\r
360 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
\r
361 pEMCAL[species] = -1;
\r
364 if( GetPtBin(pt) > -1 ){
\r
367 if(eop < 0.05) eop = 0.05;
\r
368 if(eop > 2.00) eop = 2.00;
\r
370 for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
\r
372 AliPID::EParticleType type = AliPID::EParticleType(species);
\r
374 // get nsigma value for each particle type at this E/p value
\r
375 nsigma = GetNumberOfSigmas(pt,eop,type,charge);
\r
376 sigma = GetExpectedSigma(pt,type,charge);
\r
378 // electrons (standard Gaussian calculation of probabilities)
\r
379 if(type == AliPID::kElectron){
\r
380 if (TMath::Abs(nsigma) > fRange) {
\r
381 pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
\r
384 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
\r
389 // E/p < 0.5 --> return probability below E/p = 0.5
\r
390 if ( nsigma == -999){
\r
391 pEMCAL[species] = GetLowProb(pt,type,charge);
\r
393 // E/p > 1.5 --> return probability above E/p = 1.5
\r
394 else if ( nsigma == 999){
\r
395 pEMCAL[species] = GetHighProb(pt,type,charge);
\r
397 // in parametrized region --> calculate probability for corresponding Gauss curve
\r
399 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
\r
401 // normalize to total probability == 1
\r
402 pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
\r
407 // return the electron probability
\r
408 proba = pEMCAL[AliPID::kElectron];
\r