1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 //////////////////////////////////////////////////////////////////////////
17 // AliEMCALPIDResponse //
19 // EMCAL class to perfom PID //
20 // This is a prototype and still under development //
22 // ---------------------------------------------------------------------//
23 // GetNumberOfSigmas(): //
25 // Electrons: Number of Sigmas for E/p value //
26 // Parametrization of LHC11a (after recalibration) //
29 // Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5) //
30 // --> return +/- 99 //
32 // --> return nsigma (parametrization of LHC10e) //
34 // NO Parametrization (outside pT range): --> return -999 //
36 // ---------------------------------------------------------------------//
37 // ComputeEMCALProbability(): //
39 // Electrons: Probability from Gaussian distribution //
42 // Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5) //
43 // --> probability to find particles below or above thr. //
45 // --> Probability from Gaussian distribution //
46 // (proper normalization to each other?) //
48 // NO Parametrization (outside pT range): --> return kFALSE //
49 //////////////////////////////////////////////////////////////////////////
54 #include "AliEMCALPIDResponse.h" //class header
58 ClassImp(AliEMCALPIDResponse)
60 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
61 AliEMCALPIDResponse::AliEMCALPIDResponse():
67 // The default constructor
71 fNorm = new TF1("fNorm","gaus",-20,20);
73 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
74 AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
77 fkPIDParams(other.fkPIDParams)
80 // The copy constructor
84 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85 AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
88 // The assignment operator
91 if(this == &other) return *this;
94 TObject::operator=(other);
96 fkPIDParams = other.fkPIDParams;
101 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
102 AliEMCALPIDResponse::~AliEMCALPIDResponse() {
107 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
108 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n, Int_t charge) const {
110 // Calculates the expected sigma of the PID signal as the function of
111 // the information stored in the track, for the specified particle type
118 if( charge != -1 && charge != 1){
122 // Get the parameters for this particle type and pt
123 const TVectorD *params = GetParams(n, pt, charge);
125 // IF not in momentum range, NULL is returned --> return default value
126 if(!params) return norm;
128 Double_t mean = (*params)[2]; // mean value of Gausiian parametrization
129 Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
130 Double_t eopMin = (*params)[4]; // min E/p value for parametrization
131 Double_t eopMax = (*params)[5]; // max E/p value for parametrization
132 Double_t probLow = (*params)[6]; // probability to be below eopMin
133 Double_t probHigh = (*params)[7]; // probability to be above eopMax
135 // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
136 fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*sigma*sigma),mean,sigma);
137 norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);
141 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
142 Double_t AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt, Float_t eop, AliPID::EParticleType n, Int_t charge) const {
144 Double_t nsigma = -999.;
147 if( charge != -1 && charge != 1){
151 // Get the parameters for this particle type and pt
152 const TVectorD *params = GetParams(n, pt, charge);
154 // IF not in momentum range, NULL is returned --> return default value
155 if(!params) return nsigma;
157 Double_t mean = (*params)[2]; // mean value of Gausiian parametrization
158 Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
159 Double_t eopMin = (*params)[4]; // min E/p value for parametrization
160 Double_t eopMax = (*params)[5]; // max E/p value for parametrization
163 if(n == AliPID::kElectron){
164 if(sigma != 0) nsigma = (eop - mean) / sigma;
170 nsigma = -99; // not parametrized (smaller than eopMin)
171 else if ( eop > eopMax )
172 nsigma = 99.; // not parametrized (bigger than eopMax)
174 if(sigma != 0) nsigma = (eop - mean) / sigma;
181 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
182 Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
185 Double_t fRange = 5.0; // hardcoded (???)
186 Double_t nsigma = 0.0;
190 if( charge != -1 && charge != 1){
195 // default value (will be returned, if pt below threshold)
196 for (Int_t species = 0; species < nSpecies; species++) {
197 pEMCAL[species] = 1./nSpecies;
201 if(eop < 0.05) eop = 0.05;
202 if(eop > 2.00) eop = 2.00;
204 for (Int_t species = 0; species < nSpecies; species++) {
206 AliPID::EParticleType type = AliPID::EParticleType(species);
208 // Get the parameters for this particle type and pt
209 const TVectorD *params = GetParams(species, pt, charge);
211 // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
212 if(!params) return kFALSE;
214 Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
215 Double_t probLow = (*params)[6]; // probability to be below eopMin
216 Double_t probHigh = (*params)[7]; // probability to be above eopMax
218 // get nsigma value for each particle type at this E/p value
219 nsigma = GetNumberOfSigmas(pt,eop,type,charge);
221 // electrons (standard Gaussian calculation of probabilities)
222 if(type == AliPID::kElectron){
223 if (TMath::Abs(nsigma) > fRange) {
224 pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
227 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
232 // E/p < eopMin --> return probability below E/p = eopMin
234 pEMCAL[species] = probLow;
236 // E/p > eopMax --> return probability above E/p = eopMax
237 else if ( nsigma == 99){
238 pEMCAL[species] = probHigh;
240 // in parametrized region --> calculate probability for corresponding Gauss curve
242 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
244 // normalize to total probability == 1
245 pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
253 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
254 const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const {
256 // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
264 // 6 = probLow (not used for electrons)
265 // 7 = probHigh (not used for electrons)
268 if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
269 if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
271 TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
272 if(!particlePar) return NULL;
274 TIter parIter(particlePar);
275 const TVectorD *parameters = NULL;
276 Double_t momMin = 0.;
277 Double_t momMax = 0.;
279 while((parameters = static_cast<const TVectorD *>(parIter()))){
281 momMin = (*parameters)[0];
282 momMax = (*parameters)[1];
284 if( fPt > momMin && fPt < momMax ) return parameters;
287 AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));