]>
Commit | Line | Data |
---|---|---|
da32329d AM |
1 | |
2 | /* | |
3 | <one line to give the program's name and a brief idea of what it does.> | |
4 | Copyright (C) <year> <name of author> | |
5 | ||
6 | p This program is free software: you can redistribute it and/or modify | |
7 | it under the terms of the GNU General Public License as published by | |
8 | the Free Software Foundation, either version 3 of the License, or | |
9 | (at your option) any later version. | |
10 | ||
11 | This program is distributed in the hope that it will be useful, | |
12 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
14 | GNU General Public License for more details. | |
15 | ||
16 | You should have received a copy of the GNU General Public License | |
17 | along with this program. If not, see <http://www.gnu.org/licenses/>. | |
18 | ||
19 | */ | |
20 | ||
21 | #include "spectrum.h" | |
22 | #include <cmath> | |
23 | #include "beambeamsystem.h" | |
24 | #include <randomgenerator.h> | |
25 | #include <iostream> | |
26 | ||
27 | spectrum::spectrum(beamBeamSystem *bbs) : | |
28 | _bMin(5.0) | |
29 | ,_bMax(128000.0) | |
30 | ,_nBbins(6400) | |
31 | ,_probOfBreakup(_nBbins) | |
32 | ,_beamBeamSystem(bbs) | |
33 | ,_nK(10000) | |
34 | ,_fnSingle(_nK) | |
35 | ,_fnDouble(_nK) | |
36 | ,_fnSingleCumulative(_nK+1) | |
37 | ,_fnDoubleCumulative(_nK+1) | |
38 | ,_fnDoubleInt(_nK) | |
39 | ,_fnDoubleIntCumulative(_nK+1) | |
40 | ,_eGamma(_nK+1) | |
41 | ,_eGammaMin(6.0) | |
42 | ,_eGammaMax(600000.0) | |
43 | ,_zTarget(82) | |
44 | ,_aTarget(278) | |
45 | ,_hadBreakProbCalculated(false) | |
46 | { | |
47 | _eGamma.resize(_nK+1); | |
48 | _probOfBreakup.resize(_nBbins); | |
49 | } | |
50 | ||
51 | int spectrum::generateKsingle() | |
52 | { | |
53 | ||
54 | _fnSingle.resize(_nK); | |
55 | _fnSingleCumulative.resize(_nK+1); | |
56 | ||
57 | double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK); | |
58 | ||
59 | double egamma = _eGammaMin; | |
60 | for (int i = 0; i < _nK+1; i++) | |
61 | { | |
62 | _eGamma[i] = egamma; | |
63 | egamma = egamma * eg_inc; | |
64 | } | |
65 | egamma = _eGammaMin; | |
66 | ||
67 | double fnorm = 0; | |
68 | ||
69 | ||
70 | if (_hadBreakProbCalculated == false) | |
71 | { | |
72 | _hadBreakProbCalculated = generateBreakupProbabilities(); | |
73 | } | |
74 | double binc = exp((log(_bMax/_bMin))/(double)_nBbins); | |
75 | ||
76 | for (int i = 0; i < _nK; i++) | |
77 | { | |
78 | double b = _bMin; | |
79 | ||
80 | double bint = 0.0; | |
81 | ||
82 | double f1 = 0; | |
83 | double f2 = 0; | |
84 | ||
85 | for (int j = 0; j < _nBbins - 1; j++) | |
86 | { | |
87 | double bold = b; | |
88 | if (j == 0) | |
89 | { | |
90 | //f1 = fBeamBeamSystem->getBeam1().nofe(egamma, b)*GetSigma(egamma)*fProbOfBreakup[j]*b; | |
91 | f1 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j]*b; | |
92 | //std::cout << fProbOfBreakup[j] << std::endl; | |
93 | } | |
94 | else | |
95 | { | |
96 | f1 = f2; | |
97 | } | |
98 | b = b*binc; | |
99 | // f2 = fBeamBeamSystem->getBeam1().nofe(egamma, b)*GetSigma(egamma)*fProbOfBreakup[j+1]*b;; | |
100 | f2 = getTransformedNofe(egamma, b)*getSigma(egamma)*_probOfBreakup[j+1]*b;; | |
101 | bint = bint + 0.5*(f1+f2)*(b-bold); | |
102 | } | |
103 | bint = 2.0*starlightConstants::pi*bint; | |
104 | if (i == 0) | |
105 | { | |
106 | fnorm = 1.0/bint; | |
107 | } | |
108 | _fnSingle[i] = bint*(_eGamma[i+1]-_eGamma[i]); | |
109 | ||
110 | egamma = egamma*eg_inc; | |
111 | } | |
112 | ||
113 | _fnSingleCumulative[0] = 0.00; | |
114 | for (int i = 0; i < _nK; i++) | |
115 | { | |
116 | _fnSingleCumulative[i+1] = _fnSingleCumulative[i]+_fnSingle[i]; | |
117 | } | |
118 | ||
119 | double fnormfactor = 1.00/_fnSingleCumulative[_nK]; | |
120 | for (int i = 0; i < _nK; i++) | |
121 | { | |
122 | _fnSingleCumulative[i+1] = fnormfactor*_fnSingleCumulative[i+1]; | |
123 | } | |
124 | ||
125 | return 0; | |
126 | ||
127 | } | |
128 | ||
129 | int spectrum::generateKdouble() | |
130 | { | |
131 | //Quick fix for now TODO: Fix it! | |
132 | _nK = 100; | |
133 | ||
134 | _fnDouble.resize(_nK); | |
135 | _fnDoubleInt.resize(_nK); | |
136 | _fnDoubleIntCumulative.resize(_nK+1); | |
137 | _fnDoubleCumulative.resize(_nK+1); | |
138 | for (int i = 0; i < _nK; i++) | |
139 | { | |
140 | _fnDouble[i].resize(_nK); | |
141 | _fnDoubleCumulative[i].resize(_nK+1); | |
142 | } | |
143 | _fnDoubleCumulative[_nK].resize(_nK+1); | |
144 | ||
145 | double eg_inc = exp(log(_eGammaMax/_eGammaMin)/(double)_nK); | |
146 | double egamma1 = _eGammaMin; | |
147 | double egamma2 = _eGammaMin; | |
148 | ||
149 | for (int i = 0; i < _nK+1; i++) | |
150 | { | |
151 | _eGamma[i] = egamma1; | |
152 | egamma1 = egamma1 * eg_inc; | |
153 | } | |
154 | egamma1 = _eGammaMin; | |
155 | ||
156 | double fnorm = 0; | |
157 | ||
158 | if (_hadBreakProbCalculated == false) | |
159 | { | |
160 | _hadBreakProbCalculated = generateBreakupProbabilities(); | |
161 | } | |
162 | ||
163 | double binc = exp((log(_bMax/_bMin))/(double)_nBbins); | |
164 | ||
165 | int nbbins = _nBbins; | |
166 | ||
167 | //double b_min = _bMin; | |
168 | //double b_max = _bMax; | |
169 | ||
170 | for (int i = 0; i < _nK; i++) | |
171 | { | |
172 | ||
173 | egamma2 = _eGammaMin; | |
174 | //double sum_over_k = 0.0; | |
175 | ||
176 | for (int j = 0; j < _nK; j++) | |
177 | { | |
178 | double bint = 0.0; | |
179 | double b = _bMin; | |
180 | double f1 = 0; | |
181 | double f2 = 0; | |
182 | ||
183 | for (int k = 0; k < nbbins - 1; k++) | |
184 | { | |
185 | double bold = b; | |
186 | ||
187 | if (k == 0) | |
188 | { | |
189 | // f1 = fBeamBeamSystem->getBeam1().nofe(egamma1, b) * fBeamBeamSystem->getBeam2().nofe(egamma2, b) | |
190 | // * GetSigma(egamma1) * GetSigma(egamma2) *fProbOfBreakup[j]*b; | |
191 | f1 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b) | |
192 | * getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k]*b; } | |
193 | else | |
194 | { | |
195 | f1 = f2; | |
196 | } | |
197 | b = b*binc; | |
198 | // f2 = fBeamBeamSystem->getBeam1().nofe(egamma1, b) * fBeamBeamSystem->getBeam2().nofe(egamma2, b) | |
199 | // * GetSigma(egamma1) * GetSigma(egamma2) *fProbOfBreakup[j+1]*b; | |
200 | f2 = getTransformedNofe(egamma1, b) * getTransformedNofe(egamma2, b) | |
201 | * getSigma(egamma1) * getSigma(egamma2) *_probOfBreakup[k+1]*b; | |
202 | bint = bint + 0.5*(f1+f2)*(b-bold); | |
203 | } | |
204 | bint = 2.0*starlightConstants::pi*bint; | |
205 | _fnDouble[i][j] = bint * (_eGamma[i+1] - _eGamma[i]) * (_eGamma[j+1] - _eGamma[j]); | |
206 | egamma2 = egamma2 * eg_inc; | |
207 | } | |
208 | egamma1 = egamma1 * eg_inc; | |
209 | } | |
210 | ||
211 | for (int i = 0; i < _nK; i++) | |
212 | { | |
213 | _fnDoubleInt[i] = 0.0; | |
214 | for (int j = 0; j < _nK; j++) | |
215 | { | |
216 | _fnDoubleInt[i] = _fnDoubleInt[i] + _fnDouble[i][j]; | |
217 | } | |
218 | } | |
219 | ||
220 | _fnDoubleIntCumulative[0] = 0.0; | |
221 | for (int i = 1; i < _nK+1; i++) | |
222 | { | |
223 | _fnDoubleIntCumulative[i] = _fnDoubleIntCumulative[i-1] + _fnDoubleInt[i-1]; | |
224 | } | |
225 | ||
226 | fnorm = 1.0/_fnDoubleIntCumulative[_nK]; | |
227 | for (int i = 0; i < _nK+1; i++) | |
228 | { | |
229 | _fnDoubleIntCumulative[i] = fnorm * _fnDoubleIntCumulative[i]; | |
230 | } | |
231 | ||
232 | return 0; | |
233 | } | |
234 | ||
235 | double spectrum::drawKsingle() | |
236 | { | |
237 | double xtest = 0; | |
238 | int itest = 0; | |
239 | double egamma = 0.0; | |
240 | ||
241 | xtest = randyInstance.Rndom(); | |
242 | while (xtest > _fnSingleCumulative[itest]) | |
243 | { | |
244 | itest++; | |
245 | } | |
246 | itest = itest - 1; | |
247 | ||
248 | if (itest >= _nK || itest < 0) | |
249 | { | |
250 | std::cerr << "ERROR: itest: " << itest << std::endl; | |
251 | ||
252 | } | |
253 | ||
254 | double delta_f = xtest - _fnSingleCumulative[itest]; | |
255 | if (delta_f <= 0.0) | |
256 | { | |
257 | std::cout << "WARNING: delta_f: " << delta_f << std::endl; | |
258 | return -1; | |
259 | } | |
260 | double dE = _eGamma[itest+1] - _eGamma[itest]; | |
261 | double dF = _fnSingleCumulative[itest+1] - _fnSingleCumulative[itest]; | |
262 | ||
263 | double delta_e = delta_f*dE/dF; | |
264 | ||
265 | if (delta_e > (_eGamma[itest+1] - _eGamma[itest])) | |
266 | { | |
267 | std::cerr << "ERROR: delta_E: " << delta_e << std::endl; | |
268 | } | |
269 | ||
270 | egamma = _eGamma[itest] + delta_e; | |
271 | return egamma; | |
272 | } | |
273 | ||
274 | void spectrum::drawKdouble(float& egamma1, float& egamma2) | |
275 | { | |
276 | double xtest1 = 0.0; | |
277 | double xtest2 = 0.0; | |
278 | int itest1 = 0; | |
279 | int itest2 = 0; | |
280 | ||
281 | xtest1 = randyInstance.Rndom(); | |
282 | ||
283 | while (xtest1 > _fnDoubleIntCumulative[itest1]) | |
284 | { | |
285 | itest1++; | |
286 | } | |
287 | itest1 = itest1 - 1; | |
288 | ||
289 | if (itest1 >= _nK || itest1 < 0) | |
290 | { | |
291 | std::cerr << "ERROR: itest1: " << itest1 << std::endl; | |
292 | } | |
293 | double delta_f = xtest1 - _fnDoubleIntCumulative[itest1]; | |
294 | ||
295 | if (delta_f <= 0.0) | |
296 | { | |
297 | std::cout << "WARNING: delta_f: " << delta_f << std::endl; | |
298 | } | |
299 | ||
300 | double dE = _eGamma[itest1+1] - _eGamma[itest1]; | |
301 | double dF = _fnDoubleIntCumulative[itest1+1] - _fnDoubleIntCumulative[itest1]; | |
302 | ||
303 | double delta_e = delta_f*dE/dF; | |
304 | ||
305 | if (delta_e > (_eGamma[itest1+1] - _eGamma[itest1])) | |
306 | { | |
307 | std::cerr << "ERROR: delta_E: " << delta_e << std::endl; | |
308 | } | |
309 | ||
310 | egamma1 = _eGamma[itest1] + delta_e; | |
311 | ||
312 | // Second gamma | |
313 | ||
314 | // double reldw = delta_e/(fEGamma[itest1+1] - fEGamma[itest1]); | |
315 | // std::vector<double> fn_second(fNK); | |
316 | std::vector<double> fn_second_cumulative(_nK+1); | |
317 | ||
318 | // for(int i = 0; i < fNK; i++) | |
319 | // { | |
320 | // fn_second[i] = fFnDouble[itest1][i] + (fFnDouble[itest1+1][i] - fFnDouble[itest1][i])*reldw; | |
321 | // } | |
322 | ||
323 | fn_second_cumulative[0] = 0.0; | |
324 | for(int i = 1; i < _nK+1; i++) | |
325 | { | |
326 | // fn_second_cumulative[i] = fn_second_cumulative[i-1] + fn_second[i-1]; //TODO:check indexing | |
327 | fn_second_cumulative[i] = fn_second_cumulative[i-1] + _fnDouble[itest1][i-1]; | |
328 | } | |
329 | ||
330 | double norm_factor = 1.0/fn_second_cumulative[_nK]; | |
331 | for(int i = 0; i < _nK+1; i++) | |
332 | { | |
333 | fn_second_cumulative[i] = norm_factor*fn_second_cumulative[i]; | |
334 | } | |
335 | ||
336 | xtest2 = randyInstance.Rndom(); | |
337 | ||
338 | while (xtest2 > fn_second_cumulative[itest2]) | |
339 | { | |
340 | itest2++; | |
341 | } | |
342 | itest2 = itest2 - 1; | |
343 | ||
344 | if (itest2 >= _nK || itest2 < 0) | |
345 | { | |
346 | std::cerr << "ERROR: itest2: " << itest2 << std::endl; | |
347 | } | |
348 | delta_f = xtest2 - fn_second_cumulative[itest2]; | |
349 | ||
350 | if (delta_f <= 0.0) | |
351 | { | |
352 | std::cout << "WARNING: delta_f: " << delta_f << std::endl; | |
353 | } | |
354 | ||
355 | dE = _eGamma[itest2+1] - _eGamma[itest2]; | |
356 | dF = fn_second_cumulative[itest2+1] - fn_second_cumulative[itest2]; | |
357 | ||
358 | delta_e = delta_f*dE/dF; | |
359 | ||
360 | if (delta_e > (_eGamma[itest2+1] - _eGamma[itest2])) | |
361 | { | |
362 | std::cerr << "ERROR: delta_E: " << delta_e << std::endl; | |
363 | } | |
364 | ||
365 | egamma2 = _eGamma[itest2] + delta_e; | |
366 | ||
367 | return; | |
368 | } | |
369 | ||
370 | ||
371 | bool spectrum::generateBreakupProbabilities() | |
372 | { | |
373 | ||
374 | int nbbins = _nBbins; | |
375 | ||
376 | double b_min = _bMin; | |
377 | //double b_max = _bMax; | |
378 | ||
379 | // double binc = (log(b_max/b_min))/(double)nbbins; | |
380 | double binc = exp((log(_bMax/_bMin))/(double)_nBbins); | |
381 | ||
382 | double b = b_min; | |
383 | ||
384 | _probOfBreakup.resize(nbbins); | |
385 | ||
386 | for (int i = 0; i < nbbins; i++) | |
387 | { | |
388 | double bimp = b; | |
389 | double rhad = 0; | |
390 | rhad = _beamBeamSystem->probabilityOfBreakup(bimp); | |
391 | _probOfBreakup[i] = exp(-rhad); | |
392 | b = b*binc; | |
393 | } | |
394 | return true; | |
395 | } | |
396 | ||
397 | double spectrum::getFnSingle(double egamma) const | |
398 | { | |
399 | double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK)); | |
400 | int i1 = log(egamma/_eGammaMin)/log(eginc); | |
401 | int i2 = i1 + 1; | |
402 | double fnSingle = 0.0; | |
403 | ||
404 | if (i1 < 0 || i2 > _nK) | |
405 | { | |
406 | std::cout << "I1, I2 out of bounds. Egamma = " << egamma << std::endl; | |
407 | std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl; | |
408 | fnSingle = 0.0; | |
409 | } | |
410 | else | |
411 | { | |
412 | double dE = _eGamma[i2] - _eGamma[i1]; | |
413 | double eFrac = (egamma - _eGamma[i1])/dE; | |
414 | ||
415 | if (eFrac < 0.0 || eFrac > 1.0) | |
416 | { | |
417 | std::cout << "WARNING: Efrac = " << eFrac << std::endl; | |
418 | } | |
419 | fnSingle = (1.0 - eFrac)*_fnSingle[i1] + eFrac*_fnSingle[i2]; | |
420 | } | |
421 | return fnSingle; | |
422 | } | |
423 | ||
424 | double spectrum::getFnDouble(double egamma1, double egamma2) const | |
425 | { | |
426 | double eginc = exp(log(_eGammaMax/_eGammaMin)/static_cast<double>(_nK)); | |
427 | int i1 = log(egamma1/_eGammaMin)/log(eginc); | |
428 | int i2 = i1 + 1; | |
429 | double fnDouble = 0.0; | |
430 | ||
431 | if (i1 < 0 || i2 > _nK) | |
432 | { | |
433 | std::cout << "I1, I2 out of bounds. Egamma1 = " << egamma1 << std::endl; | |
434 | std::cout << "I1, I2 = " << i1 << ", " << i2 << std::endl; | |
435 | fnDouble = 0.0; | |
436 | return fnDouble; | |
437 | } | |
438 | ||
439 | int j1 = log(egamma2/_eGammaMin)/log(eginc); | |
440 | int j2 = j1 + 1; | |
441 | ||
442 | if (j1 < 0 || j2 > _nK) | |
443 | { | |
444 | std::cout << "J1, J2 out of bounds. Egamma2 = " << egamma2 << std::endl; | |
445 | std::cout << "J1, J2 = " << j1 << ", " << j2 << std::endl; | |
446 | fnDouble = 0.0; | |
447 | return fnDouble; | |
448 | } | |
449 | ||
450 | double dE1 = _eGamma[i2] - _eGamma[i1]; | |
451 | double eFrac1 = (egamma1 - _eGamma[i1])/dE1; | |
452 | ||
453 | if (eFrac1 < 0.0 || eFrac1 > 1.0) | |
454 | { | |
455 | std::cout << "WARNING: Efrac1 = " << eFrac1 << std::endl; | |
456 | } | |
457 | ||
458 | double ptemp1 = (1.0 - eFrac1)*_fnDouble[i1][j1] + eFrac1*_fnDouble[i2][j1]; | |
459 | double ptemp2 = (1.0 - eFrac1)*_fnDouble[i1][j2] + eFrac1*_fnDouble[i2][j2]; | |
460 | ||
461 | double dE2 = _eGamma[j2] - _eGamma[j1]; | |
462 | double eFrac2 = (egamma2 - _eGamma[j1])/dE2; | |
463 | ||
464 | if (eFrac2 < 0.0 || eFrac2 > 1.0) | |
465 | { | |
466 | std::cout << "WARNING: Efrac2 = " << eFrac2 << std::endl; | |
467 | } | |
468 | ||
469 | fnDouble = (1.0 - eFrac2)*ptemp1 + eFrac2*ptemp2; | |
470 | ||
471 | return fnDouble; | |
472 | ||
473 | } | |
474 | ||
475 | double spectrum::getTransformedNofe(double egamma, double b) | |
476 | { | |
477 | double factor = 1.0/(2.0*_beamBeamSystem->beamLorentzGamma()); | |
478 | double res = factor * _beamBeamSystem->beam1().photonFlux(b, egamma*factor); | |
479 | ||
480 | return res; | |
481 | } | |
482 | ||
483 | ||
484 | ||
485 |