Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
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 : //////////////////////////////////////////////////////////////////////////
16 : // //
17 : // AliEMCALPIDResponse //
18 : // //
19 : // EMCAL class to perfom PID //
20 : // This is a prototype and still under development //
21 : // //
22 : // ---------------------------------------------------------------------//
23 : // GetNumberOfSigmas(): //
24 : // //
25 : // Electrons: Number of Sigmas for E/p value //
26 : // Parametrization of LHC11a (after recalibration) //
27 : // //
28 : // NON electrons: //
29 : // Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5) //
30 : // --> return +/- 99 //
31 : // Otherwise //
32 : // --> return nsigma (parametrization of LHC10e) //
33 : // //
34 : // NO Parametrization (outside pT range): --> return -999 //
35 : // //
36 : // ---------------------------------------------------------------------//
37 : // ComputeEMCALProbability(): //
38 : // //
39 : // Electrons: Probability from Gaussian distribution //
40 : // //
41 : // NON electrons: //
42 : // Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5) //
43 : // --> probability to find particles below or above thr. //
44 : // Otherwise //
45 : // --> Probability from Gaussian distribution //
46 : // (proper normalization to each other?) //
47 : // //
48 : // NO Parametrization (outside pT range): --> return kFALSE //
49 : //////////////////////////////////////////////////////////////////////////
50 :
51 : #include <TF1.h>
52 : #include <TMath.h>
53 :
54 : #include "AliEMCALPIDResponse.h" //class header
55 :
56 : #include "AliLog.h"
57 :
58 176 : ClassImp(AliEMCALPIDResponse)
59 :
60 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
61 : AliEMCALPIDResponse::AliEMCALPIDResponse():
62 10 : TObject(),
63 10 : fNorm(NULL),
64 10 : fCurrCentrality(-1.),
65 10 : fkPIDParams(NULL)
66 50 : {
67 : //
68 : // The default constructor
69 : //
70 :
71 :
72 30 : fNorm = new TF1("fNorm","gaus",-20,20);
73 20 : }
74 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
75 : AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
76 0 : TObject(other),
77 0 : fNorm(other.fNorm),
78 0 : fCurrCentrality(other.fCurrCentrality),
79 0 : fkPIDParams(other.fkPIDParams)
80 0 : {
81 : //
82 : // The copy constructor
83 : //
84 :
85 0 : }
86 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87 : AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
88 : {
89 : //
90 : // The assignment operator
91 : //
92 :
93 0 : if(this == &other) return *this;
94 :
95 : // Make copy
96 0 : TObject::operator=(other);
97 0 : fNorm = other.fNorm;
98 0 : fCurrCentrality = other.fCurrCentrality;
99 0 : fkPIDParams = other.fkPIDParams;
100 :
101 :
102 0 : return *this;
103 0 : }
104 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
105 40 : AliEMCALPIDResponse::~AliEMCALPIDResponse() {
106 :
107 20 : delete fNorm;
108 :
109 20 : }
110 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111 : Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n, Int_t charge) const {
112 : //
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
115 : //
116 : //
117 :
118 : Double_t norm = 1.;
119 :
120 : // Check the charge
121 0 : if( charge != -1 && charge != 1){
122 0 : return norm;
123 : }
124 :
125 : // Get the parameters for this particle type and pt
126 0 : const TVectorD *params = GetParams(n, pt, charge);
127 :
128 : // IF not in momentum range, NULL is returned --> return default value
129 0 : if(!params) return norm;
130 :
131 0 : Double_t mean = (*params)[2]; // mean value of Gausiian parametrization
132 0 : Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
133 0 : Double_t eopMin = (*params)[4]; // min E/p value for parametrization
134 0 : Double_t eopMax = (*params)[5]; // max E/p value for parametrization
135 0 : Double_t probLow = (*params)[6]; // probability to be below eopMin
136 0 : Double_t probHigh = (*params)[7]; // probability to be above eopMax
137 :
138 : // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
139 0 : fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*sigma*sigma),mean,sigma);
140 0 : norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);
141 :
142 : return norm;
143 0 : }
144 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
145 : Double_t AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt, Float_t eop, AliPID::EParticleType n, Int_t charge) const {
146 :
147 : Double_t nsigma = -999.;
148 :
149 : // Check the charge
150 0 : if( charge != -1 && charge != 1){
151 0 : return nsigma;
152 : }
153 :
154 : // Get the parameters for this particle type and pt
155 0 : const TVectorD *params = GetParams(n, pt, charge);
156 :
157 : // IF not in momentum range, NULL is returned --> return default value
158 0 : if(!params) return nsigma;
159 :
160 0 : Double_t mean = (*params)[2]; // mean value of Gausiian parametrization
161 0 : Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
162 0 : Double_t eopMin = (*params)[4]; // min E/p value for parametrization
163 0 : Double_t eopMax = (*params)[5]; // max E/p value for parametrization
164 :
165 : // if electron
166 0 : if(n == AliPID::kElectron){
167 0 : if(sigma != 0) nsigma = (eop - mean) / sigma;
168 : }
169 :
170 : // if NON electron
171 : else{
172 0 : if ( eop < eopMin )
173 0 : nsigma = -99; // not parametrized (smaller than eopMin)
174 0 : else if ( eop > eopMax )
175 0 : nsigma = 99.; // not parametrized (bigger than eopMax)
176 : else{
177 0 : if(sigma != 0) nsigma = (eop - mean) / sigma;
178 : }
179 : }
180 :
181 : return nsigma;
182 :
183 0 : }
184 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
185 : Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
186 : //
187 : //
188 : Double_t fRange = 5.0; // hardcoded (???)
189 : Double_t nsigma = 0.0;
190 :
191 :
192 : // Check the charge
193 0 : if( charge != -1 && charge != 1){
194 0 : return kFALSE;
195 : }
196 :
197 :
198 : // default value (will be returned, if pt below threshold)
199 0 : for (Int_t species = 0; species < nSpecies; species++) {
200 0 : pEMCAL[species] = 1./nSpecies;
201 : }
202 :
203 : // set E/p range
204 0 : if(eop < 0.05) eop = 0.05;
205 0 : if(eop > 2.00) eop = 2.00;
206 :
207 0 : for (Int_t species = 0; species < nSpecies; species++) {
208 :
209 : AliPID::EParticleType type = AliPID::EParticleType(species);
210 :
211 : // Get the parameters for this particle type and pt
212 0 : const TVectorD *params = GetParams(species, pt, charge);
213 :
214 : // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
215 0 : if(!params) return kFALSE;
216 :
217 0 : Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
218 0 : Double_t probLow = (*params)[6]; // probability to be below eopMin
219 0 : Double_t probHigh = (*params)[7]; // probability to be above eopMax
220 :
221 : // get nsigma value for each particle type at this E/p value
222 0 : nsigma = GetNumberOfSigmas(pt,eop,type,charge);
223 :
224 : // electrons (standard Gaussian calculation of probabilities)
225 0 : if(type == AliPID::kElectron){
226 0 : if (TMath::Abs(nsigma) > fRange) {
227 0 : pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
228 0 : }
229 : else{
230 0 : pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
231 : }
232 : }
233 : //NON electrons
234 : else{
235 : // E/p < eopMin --> return probability below E/p = eopMin
236 0 : if ( nsigma == -99){
237 0 : pEMCAL[species] = probLow;
238 0 : }
239 : // E/p > eopMax --> return probability above E/p = eopMax
240 0 : else if ( nsigma == 99){
241 0 : pEMCAL[species] = probHigh;
242 0 : }
243 : // in parametrized region --> calculate probability for corresponding Gauss curve
244 : else{
245 0 : pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
246 :
247 : // normalize to total probability == 1
248 0 : pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
249 : }
250 : }
251 0 : }
252 :
253 0 : return kTRUE;
254 :
255 0 : }
256 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
257 : const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const {
258 : //
259 : // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
260 : //
261 : // 0 = momMin
262 : // 1 = momMax
263 : // 2 = mean of Gaus
264 : // 3 = sigma of Gaus
265 : // 4 = eopLow
266 : // 5 = eopHig
267 : // 6 = probLow (not used for electrons)
268 : // 7 = probHigh (not used for electrons)
269 : //
270 : // for PbPb the parametrization is done centrality dependent (marked by TString "Centrality")
271 : // so first the correct centrality bin has to be found
272 :
273 : // **** Centrality bins (hard coded for the moment)
274 : const Int_t nCent = 7;
275 0 : Int_t centBins[nCent+1] = {0,10,20,30,40,50,70,90};
276 :
277 0 : if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
278 0 : if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
279 :
280 0 : TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
281 0 : if(!particlePar) return NULL;
282 :
283 : const TVectorD *parameters = NULL;
284 : Double_t momMin = 0.;
285 : Double_t momMax = 0.;
286 :
287 : // is the centrality dependent parametrization used
288 0 : TString arrayName = particlePar->GetName();
289 :
290 : // centrality dependent parametrization
291 0 : if(arrayName.Contains("Centrality")){
292 :
293 0 : for(Int_t iCent = 0; iCent < nCent; iCent++){
294 :
295 0 : if( fCurrCentrality > centBins[iCent] && fCurrCentrality < centBins[iCent+1] ){
296 :
297 0 : TObjArray * centPar = dynamic_cast<TObjArray *>(particlePar->At(iCent));
298 0 : if(!centPar) return NULL;
299 :
300 0 : TIter centIter(centPar);
301 : parameters = NULL;
302 : momMin = 0.;
303 : momMax = 0.;
304 :
305 0 : while((parameters = static_cast<const TVectorD *>(centIter()))){
306 :
307 0 : momMin = (*parameters)[0];
308 0 : momMax = (*parameters)[1];
309 :
310 0 : if( fPt > momMin && fPt < momMax ) return parameters;
311 :
312 : }
313 0 : }
314 : }
315 : }
316 :
317 : // NO centrality dependent parametrization
318 : else{
319 :
320 0 : TIter parIter(particlePar);
321 0 : while((parameters = static_cast<const TVectorD *>(parIter()))){
322 :
323 0 : momMin = (*parameters)[0];
324 0 : momMax = (*parameters)[1];
325 :
326 0 : if( fPt > momMin && fPt < momMax ) return parameters;
327 :
328 : }
329 0 : }
330 0 : AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));
331 :
332 0 : return parameters;
333 0 : }
|
