]>
Commit | Line | Data |
---|---|---|
b2138b40 | 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 | // // | |
00a38d07 | 48 | // NO Parametrization (outside pT range): --> return kFALSE // |
b2138b40 | 49 | ////////////////////////////////////////////////////////////////////////// |
50 | ||
51 | #include <TF1.h> | |
52 | #include <TMath.h> | |
53 | ||
54 | #include "AliEMCALPIDResponse.h" //class header | |
55 | ||
56 | #include "AliLog.h" | |
57 | ||
58 | ClassImp(AliEMCALPIDResponse) | |
59 | ||
60 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
61 | AliEMCALPIDResponse::AliEMCALPIDResponse(): | |
62 | TObject(), | |
63 | fNorm(NULL), | |
64 | fkPIDParams(NULL) | |
65 | { | |
66 | // | |
67 | // The default constructor | |
68 | // | |
69 | ||
70 | ||
71 | fNorm = new TF1("fNorm","gaus",-20,20); | |
72 | } | |
73 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
74 | AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other): | |
75 | TObject(other), | |
76 | fNorm(other.fNorm), | |
77 | fkPIDParams(other.fkPIDParams) | |
78 | { | |
79 | // | |
80 | // The copy constructor | |
81 | // | |
82 | ||
83 | } | |
84 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
85 | AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other) | |
86 | { | |
87 | // | |
88 | // The assignment operator | |
89 | // | |
90 | ||
91 | if(this == &other) return *this; | |
92 | ||
93 | // Make copy | |
94 | TObject::operator=(other); | |
95 | fNorm = other.fNorm; | |
96 | fkPIDParams = other.fkPIDParams; | |
97 | ||
98 | ||
99 | return *this; | |
100 | } | |
101 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
102 | AliEMCALPIDResponse::~AliEMCALPIDResponse() { | |
103 | ||
104 | delete fNorm; | |
105 | ||
106 | } | |
107 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
108 | Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n, Int_t charge) const { | |
109 | // | |
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 | |
112 | // | |
113 | // | |
114 | ||
115 | Double_t norm = 1.; | |
116 | ||
117 | // Check the charge | |
118 | if( charge != -1 && charge != 1){ | |
119 | return norm; | |
120 | } | |
121 | ||
122 | // Get the parameters for this particle type and pt | |
1f631618 | 123 | const TVectorD *params = GetParams(n, pt, charge); |
b2138b40 | 124 | |
125 | // IF not in momentum range, NULL is returned --> return default value | |
126 | if(!params) return norm; | |
127 | ||
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 | |
134 | ||
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); | |
138 | ||
139 | return norm; | |
140 | } | |
141 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
142 | Double_t AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt, Float_t eop, AliPID::EParticleType n, Int_t charge) const { | |
143 | ||
144 | Double_t nsigma = -999.; | |
145 | ||
146 | // Check the charge | |
147 | if( charge != -1 && charge != 1){ | |
148 | return nsigma; | |
149 | } | |
150 | ||
151 | // Get the parameters for this particle type and pt | |
1f631618 | 152 | const TVectorD *params = GetParams(n, pt, charge); |
b2138b40 | 153 | |
154 | // IF not in momentum range, NULL is returned --> return default value | |
155 | if(!params) return nsigma; | |
156 | ||
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 | |
161 | ||
162 | // if electron | |
163 | if(n == AliPID::kElectron){ | |
164 | if(sigma != 0) nsigma = (eop - mean) / sigma; | |
165 | } | |
166 | ||
167 | // if NON electron | |
168 | else{ | |
169 | if ( eop < eopMin ) | |
170 | nsigma = -99; // not parametrized (smaller than eopMin) | |
171 | else if ( eop > eopMax ) | |
172 | nsigma = 99.; // not parametrized (bigger than eopMax) | |
173 | else{ | |
174 | if(sigma != 0) nsigma = (eop - mean) / sigma; | |
175 | } | |
176 | } | |
177 | ||
178 | return nsigma; | |
179 | ||
180 | } | |
181 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
00a38d07 | 182 | Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const { |
b2138b40 | 183 | // |
184 | // | |
185 | Double_t fRange = 5.0; // hardcoded (???) | |
186 | Double_t nsigma = 0.0; | |
b2138b40 | 187 | |
188 | ||
189 | // Check the charge | |
190 | if( charge != -1 && charge != 1){ | |
00a38d07 | 191 | return kFALSE; |
b2138b40 | 192 | } |
193 | ||
194 | ||
195 | // default value (will be returned, if pt below threshold) | |
00a38d07 | 196 | for (Int_t species = 0; species < nSpecies; species++) { |
197 | pEMCAL[species] = 1./nSpecies; | |
b2138b40 | 198 | } |
199 | ||
200 | // set E/p range | |
201 | if(eop < 0.05) eop = 0.05; | |
202 | if(eop > 2.00) eop = 2.00; | |
203 | ||
00a38d07 | 204 | for (Int_t species = 0; species < nSpecies; species++) { |
b2138b40 | 205 | |
206 | AliPID::EParticleType type = AliPID::EParticleType(species); | |
207 | ||
208 | // Get the parameters for this particle type and pt | |
1f631618 | 209 | const TVectorD *params = GetParams(species, pt, charge); |
b2138b40 | 210 | |
00a38d07 | 211 | // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE |
212 | if(!params) return kFALSE; | |
b2138b40 | 213 | |
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 | |
217 | ||
218 | // get nsigma value for each particle type at this E/p value | |
219 | nsigma = GetNumberOfSigmas(pt,eop,type,charge); | |
220 | ||
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); | |
225 | } | |
226 | else{ | |
227 | pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma); | |
228 | } | |
229 | } | |
230 | //NON electrons | |
231 | else{ | |
232 | // E/p < eopMin --> return probability below E/p = eopMin | |
233 | if ( nsigma == -99){ | |
234 | pEMCAL[species] = probLow; | |
235 | } | |
236 | // E/p > eopMax --> return probability above E/p = eopMax | |
237 | else if ( nsigma == 99){ | |
238 | pEMCAL[species] = probHigh; | |
239 | } | |
240 | // in parametrized region --> calculate probability for corresponding Gauss curve | |
241 | else{ | |
242 | pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma); | |
243 | ||
244 | // normalize to total probability == 1 | |
245 | pEMCAL[species]*=GetExpectedNorm(pt,type,charge); | |
246 | } | |
247 | } | |
b2138b40 | 248 | } |
249 | ||
00a38d07 | 250 | return kTRUE; |
b2138b40 | 251 | |
252 | } | |
253 | //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | |
1f631618 | 254 | const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const { |
b2138b40 | 255 | // |
256 | // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt | |
257 | // | |
258 | // 0 = momMin | |
259 | // 1 = momMax | |
260 | // 2 = mean of Gaus | |
261 | // 3 = sigma of Gaus | |
262 | // 4 = eopLow | |
263 | // 5 = eopHig | |
264 | // 6 = probLow (not used for electrons) | |
265 | // 7 = probHigh (not used for electrons) | |
266 | // | |
267 | ||
268 | if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL; | |
1f631618 | 269 | if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons |
b2138b40 | 270 | |
271 | TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle)); | |
272 | if(!particlePar) return NULL; | |
1f631618 | 273 | |
b2138b40 | 274 | TIter parIter(particlePar); |
275 | const TVectorD *parameters = NULL; | |
276 | Double_t momMin = 0.; | |
277 | Double_t momMax = 0.; | |
278 | ||
279 | while((parameters = static_cast<const TVectorD *>(parIter()))){ | |
280 | ||
281 | momMin = (*parameters)[0]; | |
282 | momMax = (*parameters)[1]; | |
283 | ||
284 | if( fPt > momMin && fPt < momMax ) return parameters; | |
285 | ||
286 | } | |
287 | AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt)); | |
288 | ||
289 | return parameters; | |
290 | } |