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():
68 // The default constructor
72 fNorm = new TF1("fNorm","gaus",-20,20);
74 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
75 AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
78 fCurrCentrality(other.fCurrCentrality),
79 fkPIDParams(other.fkPIDParams)
82 // The copy constructor
86 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87 AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
90 // The assignment operator
93 if(this == &other) return *this;
96 TObject::operator=(other);
98 fCurrCentrality = other.fCurrCentrality;
99 fkPIDParams = other.fkPIDParams;
104 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
105 AliEMCALPIDResponse::~AliEMCALPIDResponse() {
110 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111 Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n, Int_t charge) const {
113 // Calculates the expected sigma of the PID signal as the function of
114 // the information stored in the track, for the specified particle type
121 if( charge != -1 && charge != 1){
125 // Get the parameters for this particle type and pt
126 const TVectorD *params = GetParams(n, pt, charge);
128 // IF not in momentum range, NULL is returned --> return default value
129 if(!params) return norm;
131 Double_t mean = (*params)[2]; // mean value of Gausiian parametrization
132 Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
133 Double_t eopMin = (*params)[4]; // min E/p value for parametrization
134 Double_t eopMax = (*params)[5]; // max E/p value for parametrization
135 Double_t probLow = (*params)[6]; // probability to be below eopMin
136 Double_t probHigh = (*params)[7]; // probability to be above eopMax
138 // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
139 fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*sigma*sigma),mean,sigma);
140 norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);
144 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
145 Double_t AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt, Float_t eop, AliPID::EParticleType n, Int_t charge) const {
147 Double_t nsigma = -999.;
150 if( charge != -1 && charge != 1){
154 // Get the parameters for this particle type and pt
155 const TVectorD *params = GetParams(n, pt, charge);
157 // IF not in momentum range, NULL is returned --> return default value
158 if(!params) return nsigma;
160 Double_t mean = (*params)[2]; // mean value of Gausiian parametrization
161 Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
162 Double_t eopMin = (*params)[4]; // min E/p value for parametrization
163 Double_t eopMax = (*params)[5]; // max E/p value for parametrization
166 if(n == AliPID::kElectron){
167 if(sigma != 0) nsigma = (eop - mean) / sigma;
173 nsigma = -99; // not parametrized (smaller than eopMin)
174 else if ( eop > eopMax )
175 nsigma = 99.; // not parametrized (bigger than eopMax)
177 if(sigma != 0) nsigma = (eop - mean) / sigma;
184 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
185 Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
188 Double_t fRange = 5.0; // hardcoded (???)
189 Double_t nsigma = 0.0;
193 if( charge != -1 && charge != 1){
198 // default value (will be returned, if pt below threshold)
199 for (Int_t species = 0; species < nSpecies; species++) {
200 pEMCAL[species] = 1./nSpecies;
204 if(eop < 0.05) eop = 0.05;
205 if(eop > 2.00) eop = 2.00;
207 for (Int_t species = 0; species < nSpecies; species++) {
209 AliPID::EParticleType type = AliPID::EParticleType(species);
211 // Get the parameters for this particle type and pt
212 const TVectorD *params = GetParams(species, pt, charge);
214 // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
215 if(!params) return kFALSE;
217 Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
218 Double_t probLow = (*params)[6]; // probability to be below eopMin
219 Double_t probHigh = (*params)[7]; // probability to be above eopMax
221 // get nsigma value for each particle type at this E/p value
222 nsigma = GetNumberOfSigmas(pt,eop,type,charge);
224 // electrons (standard Gaussian calculation of probabilities)
225 if(type == AliPID::kElectron){
226 if (TMath::Abs(nsigma) > fRange) {
227 pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
230 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
235 // E/p < eopMin --> return probability below E/p = eopMin
237 pEMCAL[species] = probLow;
239 // E/p > eopMax --> return probability above E/p = eopMax
240 else if ( nsigma == 99){
241 pEMCAL[species] = probHigh;
243 // in parametrized region --> calculate probability for corresponding Gauss curve
245 pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
247 // normalize to total probability == 1
248 pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
256 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
257 const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const {
259 // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
267 // 6 = probLow (not used for electrons)
268 // 7 = probHigh (not used for electrons)
270 // for PbPb the parametrization is done centrality dependent (marked by TString "Centrality")
271 // so first the correct centrality bin has to be found
273 // **** Centrality bins (hard coded for the moment)
274 const Int_t nCent = 7;
275 Int_t centBins[nCent+1] = {0,10,20,30,40,50,70,90};
277 if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
278 if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
280 TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
281 if(!particlePar) return NULL;
283 const TVectorD *parameters = NULL;
284 Double_t momMin = 0.;
285 Double_t momMax = 0.;
287 // is the centrality dependent parametrization used
288 TString arrayName = particlePar->GetName();
290 // centrality dependent parametrization
291 if(arrayName.Contains("Centrality")){
293 for(Int_t iCent = 0; iCent < nCent; iCent++){
295 if( fCurrCentrality > centBins[iCent] && fCurrCentrality < centBins[iCent+1] ){
297 TObjArray * centPar = dynamic_cast<TObjArray *>(particlePar->At(iCent));
298 if(!centPar) return NULL;
300 TIter centIter(centPar);
305 while((parameters = static_cast<const TVectorD *>(centIter()))){
307 momMin = (*parameters)[0];
308 momMax = (*parameters)[1];
310 if( fPt > momMin && fPt < momMax ) return parameters;
317 // NO centrality dependent parametrization
320 TIter parIter(particlePar);
321 while((parameters = static_cast<const TVectorD *>(parIter()))){
323 momMin = (*parameters)[0];
324 momMax = (*parameters)[1];
326 if( fPt > momMin && fPt < momMax ) return parameters;
330 AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));