1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 | |
38 | |
39 | |
40 | |
41 | |
42 | |
43 | |
44 | |
45 | |
46 | |
47 | |
48 | |
49 | |
50 | |
51 | #define COBREMS_GENERATOR_VERBOSITY1 1 |
52 | |
53 | |
54 | #include <iostream> |
55 | #include <stdio.h> |
56 | #include <CobremsGeneration.hh> |
57 | #include <boost/math/special_functions/expint.hpp> |
58 | #include <boost/math/special_functions/erf.hpp> |
59 | |
60 | const double CobremsGeneration::dpi = 3.1415926535897; |
61 | const double CobremsGeneration::me = 0.510998910e-3; |
62 | const double CobremsGeneration::alpha = 7.2973525698e-3; |
63 | const double CobremsGeneration::hbarc = 0.1973269718e-15; |
64 | |
65 | CobremsGeneration::CobremsGeneration(double Emax_GeV, double Epeak_GeV) |
66 | { |
67 | |
68 | |
69 | |
70 | fBeamEnergy = Emax_GeV; |
71 | fBeamErms = 6.0e-4; |
72 | fBeamEmittance = 2.5e-9; |
73 | fCollimatorSpotrms = 0.0005; |
74 | fCollimatorDistance = 76.0; |
75 | fCollimatorDiameter = 0.005; |
76 | fTargetThickness = 50e-6; |
77 | fTargetThetay = 0.050; |
78 | fTargetThetaz = 0; |
79 | setTargetCrystal("diamond"); |
80 | setCoherentEdge(Epeak_GeV); |
81 | fPhotonEnergyMin = 0.120; |
82 | setPolarizedFlag(false); |
83 | setCollimatedFlag(true); |
84 | |
85 | #if COBREMS_GENERATOR_VERBOSITY1 > 0 |
86 | std::cout << std::endl |
87 | << "Initialization for coherent bremsstralung calculation" |
88 | << std::endl |
89 | << " electron beam energy: " << Emax_GeV << " GeV" |
90 | << std::endl |
91 | << " primary coherent edge: " << Epeak_GeV << " GeV" |
92 | << std::endl; |
93 | #endif |
94 | } |
95 | |
96 | void CobremsGeneration::updateTargetOrientation() |
97 | { |
98 | resetTargetOrientation(); |
99 | RotateTarget(0, dpi/2, 0); |
100 | RotateTarget(0, 0, dpi/4); |
101 | RotateTarget(0, 0, -fTargetThetaz); |
102 | RotateTarget(0, -fTargetThetay, 0); |
103 | RotateTarget(-fTargetThetax, 0, 0); |
104 | } |
105 | |
106 | void CobremsGeneration::setTargetCrystal(std::string crystal) |
107 | { |
108 | |
109 | |
110 | if (crystal == "diamond") { |
111 | fTargetCrystal.name = "diamond"; |
112 | fTargetCrystal.Z = 6; |
113 | fTargetCrystal.A = 12.01; |
114 | fTargetCrystal.density = 3.534; |
115 | fTargetCrystal.lattice_constant = 3.5668e-10; |
116 | fTargetCrystal.Debye_Waller_const = 0.40e9; |
117 | } |
118 | else if (crystal == "silicon") { |
119 | fTargetCrystal.name = "silicon"; |
120 | fTargetCrystal.Z = 14; |
121 | fTargetCrystal.A = 28.09; |
122 | fTargetCrystal.density = 2.320; |
123 | fTargetCrystal.lattice_constant = 5.431e-10; |
124 | fTargetCrystal.Debye_Waller_const = 1.5e9; |
125 | } |
126 | else { |
127 | std::cerr << "Error in CobremsGeneration::setTargetCrystal - " |
128 | << "unknown crystal " << crystal << " requested, " |
129 | << "cannot continue." << std::endl; |
130 | exit(1); |
131 | } |
132 | |
133 | |
134 | fTargetCrystal.nsites = 8; |
135 | fTargetCrystal.ucell_site.clear(); |
136 | fTargetCrystal.ucell_site.push_back(lattice_vector(0.0, 0.0, 0.0)); |
137 | fTargetCrystal.ucell_site.push_back(lattice_vector(0.0, 0.5, 0.5)); |
138 | fTargetCrystal.ucell_site.push_back(lattice_vector(0.5, 0.0, 0.5)); |
139 | fTargetCrystal.ucell_site.push_back(lattice_vector(0.5, 0.5, 0.0)); |
140 | fTargetCrystal.ucell_site.push_back(lattice_vector(0.25, 0.25, 0.25)); |
141 | fTargetCrystal.ucell_site.push_back(lattice_vector(0.25, 0.75, 0.75)); |
142 | fTargetCrystal.ucell_site.push_back(lattice_vector(0.75, 0.25, 0.75)); |
143 | fTargetCrystal.ucell_site.push_back(lattice_vector(0.75, 0.75, 0.25)); |
144 | fTargetCrystal.primaryHKL = lattice_vector(2,2,0); |
145 | |
146 | |
147 | fTargetCrystal.betaFF = 111 * pow(fTargetCrystal.Z, -1/3.) / me; |
148 | |
149 | |
150 | fTargetCrystal.mosaic_spread = 20e-6; |
151 | |
152 | |
153 | fTargetCrystal.radiation_length = getTargetRadiationLength_Schiff(); |
154 | } |
155 | |
156 | double CobremsGeneration::getTargetDebyeWallerConstant(double DebyeT_K, |
157 | double T_K) |
158 | { |
159 | |
160 | |
161 | |
162 | |
163 | |
164 | |
165 | |
166 | |
167 | |
168 | |
169 | |
170 | |
171 | |
172 | |
173 | double kBoltzmann = 8.617e-14; |
174 | double amassGeV = fTargetCrystal.A * 0.932; |
175 | double A0 = 3 / (2 * amassGeV * kBoltzmann * DebyeT_K); |
176 | double Tnormal = (T_K + 0.1) / DebyeT_K; |
177 | int niter = 50; |
178 | double f = 0; |
179 | for (int iter=0; iter < niter; ++iter) { |
180 | double x = (iter + 0.5) / niter; |
181 | f += x * (1 + 2 / (exp(x / Tnormal) - 1)) / niter; |
182 | } |
183 | return A0 * f; |
184 | } |
185 | |
186 | void CobremsGeneration::printBeamlineInfo() |
187 | { |
188 | |
189 | |
190 | std::cout << " electron beam energy: " << fBeamEnergy << " GeV" |
191 | << std::endl |
192 | << " electron beam emittance: " |
193 | << fBeamEmittance * 1e9 << " mm.urad" |
194 | << std::endl |
195 | << " radiator crystal: " << fTargetCrystal.name |
196 | << ", thickness " << fTargetThickness * 1e6 << " um" |
197 | << std::endl |
198 | << " radiation length: " |
199 | << fTargetCrystal.radiation_length * 100 << " cm," |
200 | << " mosaic spread: " |
201 | << fTargetCrystal.mosaic_spread * 1e6 << " urad" |
202 | << std::endl |
203 | << " photon beam collimator half-angle: " |
204 | << fCollimatorDiameter / (2 * fCollimatorDistance) |
205 | * fBeamEnergy / me << " (m/E)" |
206 | << std::endl |
207 | << " collimator diameter: " |
208 | << fCollimatorDiameter * 100 << " cm" |
209 | << std::endl |
210 | << " crystal orientation: theta_x " |
211 | << fTargetThetax * 1e3 << " mrad" |
212 | << std::endl |
213 | << " theta_y " |
214 | << fTargetThetay * 1e3 << " mrad" |
215 | << std::endl << std::endl; |
216 | } |
217 | |
218 | void CobremsGeneration::printTargetCrystalInfo() |
219 | { |
220 | |
221 | |
222 | double kBoltzmann = 8.617e-14; |
223 | double amass = fTargetCrystal.A * 0.932; |
224 | double DebyeTheta0 = 3 / (4 * amass * fTargetCrystal.Debye_Waller_const); |
225 | double DebyeTheta300 = DebyeTheta0; |
226 | for (int i=0; i < 5; i++) { |
227 | DebyeTheta300 = DebyeTheta0 * (1 + |
228 | pow(2 * dpi * 300 * kBoltzmann / DebyeTheta300, 2) / 6); |
229 | } |
230 | |
231 | std::cout << "CobremsGeneration crystal type is " << fTargetCrystal.name |
232 | << std::endl |
233 | << " atomic number Z=" << fTargetCrystal.Z |
234 | << ", atomic weight A=" << fTargetCrystal.A << " amu" |
235 | << std::endl |
236 | << " mass density: " << fTargetCrystal.density << " g/cm^3" |
237 | << std::endl |
238 | << " radiation length: " << fTargetCrystal.radiation_length * 100 |
239 | << " cm" << std::endl |
240 | << " Debye-Waller constant: " << fTargetCrystal.Debye_Waller_const |
241 | << " /GeV^2 (" << DebyeTheta300 / kBoltzmann << " K)" |
242 | << std::endl |
243 | << " mosaic spread: " << fTargetCrystal.mosaic_spread * 1e6 |
244 | << " urad" << std::endl |
245 | << " atomic form-factor cutoff momentum: " |
246 | << sqrt(1 / fTargetCrystal.betaFF) * 1e6 << " keV" |
247 | << std::endl |
248 | << " primary lattice reflection h,k,l=" |
249 | << fTargetCrystal.primaryHKL.x << "," |
250 | << fTargetCrystal.primaryHKL.y << "," |
251 | << fTargetCrystal.primaryHKL.z |
252 | << std::endl |
253 | << " lattice constant: " << fTargetCrystal.lattice_constant * 1e9 |
254 | << " nm" << std::endl |
255 | << " occupied sites of the crystal lattice unit cell are:" |
256 | << std::endl; |
257 | for (unsigned int i=0; i < fTargetCrystal.ucell_site.size(); ++i) { |
258 | char s[100]; |
259 | snprintf(s, 100, "%4.2f %4.2f %4.2f", fTargetCrystal.ucell_site[i].x, |
260 | fTargetCrystal.ucell_site[i].y, |
261 | fTargetCrystal.ucell_site[i].z); |
262 | std::cout << " " << i + 1 << ": " << s << std::endl; |
263 | } |
264 | std::cout << " Crystal orientation matrix is:" << std::endl; |
265 | for (int i=0; i < 3; ++i) { |
266 | char s[100]; |
267 | snprintf(s, 100, "%15.12f %15.12f %15.12f", fTargetRmatrix[i][0], |
268 | fTargetRmatrix[i][1], |
269 | fTargetRmatrix[i][2]); |
270 | std::cout << " " << s << std::endl; |
271 | } |
272 | } |
273 | |
274 | void CobremsGeneration::applyBeamCrystalConvolution(int nbins, double *xvalues, |
275 | double *yvalues) |
276 | { |
277 | |
278 | |
279 | |
280 | |
281 | |
282 | |
283 | |
284 | |
285 | |
286 | |
287 | |
288 | |
289 | |
290 | |
291 | |
292 | |
293 | |
294 | |
295 | |
296 | |
297 | |
298 | |
299 | |
300 | |
301 | |
302 | |
303 | |
304 | |
305 | |
306 | |
307 | double x0 = xvalues[0]; |
308 | double x1 = xvalues[nbins - 1]; |
309 | double var0 = pow(fTargetCrystal.mosaic_spread, 2) + |
310 | pow(fBeamEmittance / fCollimatorSpotrms, 2); |
311 | double varMS = Sigma2MS(fTargetThickness); |
312 | |
313 | |
314 | |
315 | |
316 | |
317 | |
318 | double a = fTargetCrystal.lattice_constant; |
319 | double qabs = sqrt(8.0) * hbarc * 2*dpi / a; |
320 | double xfact = 2 * fBeamEnergy * qabs / (me*me); |
321 | double *norm = new double[nbins]; |
322 | double *result = new double[nbins]; |
323 | for (int j=0; j < nbins; ++j) { |
324 | norm[j] = 0; |
325 | result[j] = 0; |
326 | for (int i=0; i < nbins; ++i) { |
327 | double dx = (x1 - x0) * (j - i) / nbins; |
328 | double x = x0 + (x1 - x0) * (j + 0.5) / nbins; |
329 | double dalph = dx / xfact / pow(1 - x + 1e-99, 2); |
330 | double term; |
331 | if (varMS / var0 > 1e-4) { |
332 | term = dalph / varMS * |
333 | (boost::math::erf(dalph / sqrt(2 * (var0 + varMS))) - |
334 | boost::math::erf(dalph / sqrt(2 * var0))) + |
335 | sqrt(2 / dpi) / varMS * |
336 | (exp(-dalph*dalph / (2 * (var0 + varMS))) * |
337 | sqrt(var0 + varMS) - |
338 | exp(-dalph*dalph / (2 * var0)) * sqrt(var0)); |
339 | } |
340 | else { |
341 | term = exp(-dalph*dalph / (2 * var0)) / sqrt(2 * dpi * var0); |
342 | } |
343 | norm[j] += term; |
344 | } |
345 | } |
346 | |
347 | for (int i=0; i < nbins; ++i) { |
348 | for (int j=0; j < nbins; ++j) { |
349 | double dx = (x1 - x0) * (j - i) / nbins; |
350 | double x = x0 + (x1 - x0) * (j + 0.5) / nbins; |
351 | double dalph = dx / xfact / pow(1 - x + 1e-99, 2); |
352 | double term; |
353 | if (varMS / var0 > 1e-4) { |
354 | term = dalph / varMS * |
355 | (boost::math::erf(dalph / sqrt(2 * (var0 + varMS))) - |
356 | boost::math::erf(dalph / sqrt(2 * var0))) + |
357 | sqrt(2 / dpi) / varMS * |
358 | (exp(-dalph*dalph / (2 * (var0 + varMS))) * |
359 | sqrt(var0 + varMS) - |
360 | exp(-dalph*dalph/ (2 * var0)) * sqrt(var0)); |
361 | } |
362 | else { |
363 | term = exp(-dalph*dalph / (2 * var0)) / sqrt(2 * dpi * var0); |
364 | } |
365 | result[i] += term * yvalues[j] / norm[j]; |
366 | } |
367 | } |
368 | |
369 | for (int i=0; i < nbins; ++i) { |
370 | if (fabs(result[i]) > 1e-35) { |
371 | yvalues[i] = result[i]; |
372 | } |
373 | else { |
374 | yvalues[i] = 0; |
375 | } |
376 | } |
377 | |
378 | delete [] norm; |
379 | delete [] result; |
380 | } |
381 | |
382 | double CobremsGeneration::getTargetRadiationLength_PDG() |
383 | { |
384 | |
385 | |
386 | double Z = fTargetCrystal.Z; |
387 | double N = fTargetCrystal.nsites; |
388 | double a = fTargetCrystal.lattice_constant; |
389 | double c = alpha * Z; |
390 | double s = 4 * N * pow(alpha, 3) * pow(hbarc/(a*me), 2) / a * |
391 | (Z*Z * (log(184.15 * pow(Z, -1/3.)) - |
392 | c*c * (1 / (1 + c*c) + 0.20206 - 0.0369 * c*c + |
393 | 0.0083 * pow(c, 4) - 0.002 * pow(c, 6))) + |
394 | Z * log(1194 * pow(Z, -2/3.))); |
395 | return 1/s; |
396 | } |
397 | |
398 | double CobremsGeneration::getTargetRadiationLength_Schiff() |
399 | { |
400 | |
401 | |
402 | double Z = fTargetCrystal.Z; |
403 | double N = fTargetCrystal.nsites; |
404 | double a = fTargetCrystal.lattice_constant; |
405 | double zeta = log(1440 * pow(Z, -2/3.)) / log(183 * pow(Z, -1/3.)); |
406 | double s = 4 * N * pow(alpha, 3) * pow(hbarc/(a*me), 2) / a * |
407 | Z * (Z + zeta) * log(183 * pow(Z, -1/3.)); |
408 | return 1/s; |
409 | } |
410 | |
411 | void CobremsGeneration::setCoherentEdge(double Epeak_GeV) |
412 | { |
413 | |
414 | |
415 | |
416 | |
417 | double edge = Epeak_GeV; |
418 | double qtotal = hbarc * (2 * dpi / fTargetCrystal.lattice_constant); |
419 | lattice_vector hkl = fTargetCrystal.primaryHKL; |
420 | qtotal *= sqrt(hkl.x * hkl.x + hkl.y * hkl.y + hkl.z * hkl.z); |
421 | double qlong = edge * me*me / (2 * fBeamEnergy * (fBeamEnergy - edge)); |
422 | fTargetThetax = -qlong / qtotal; |
423 | updateTargetOrientation(); |
424 | } |
425 | |
426 | CobremsGeneration::CobremsGeneration(const CobremsGeneration &src) |
427 | { |
428 | |
429 | |
430 | fTargetCrystal = src.fTargetCrystal; |
431 | fTargetThickness = src.fTargetThickness; |
432 | fTargetThetax = src.fTargetThetax; |
433 | fTargetThetay = src.fTargetThetay; |
434 | fTargetThetaz = src.fTargetThetaz; |
435 | for (int i=0; i < 3; ++i) |
436 | for (int j=0; j < 3; ++j) |
437 | fTargetRmatrix[i][j] = src.fTargetRmatrix[i][j]; |
438 | fBeamEnergy = src.fBeamEnergy; |
439 | fBeamErms = src.fBeamErms; |
440 | fBeamEmittance = src.fBeamEmittance; |
441 | fCollimatorSpotrms = src.fCollimatorSpotrms; |
442 | fCollimatorDistance = src.fCollimatorDistance; |
443 | fCollimatorDiameter = src.fCollimatorDiameter; |
444 | fQ2theta2 = src.fQ2theta2; |
445 | fQ2weight = src.fQ2weight; |
446 | } |
447 | |
448 | CobremsGeneration &CobremsGeneration::operator=(const CobremsGeneration &src) |
449 | { |
450 | |
451 | |
452 | fTargetCrystal = src.fTargetCrystal; |
453 | fTargetThickness = src.fTargetThickness; |
454 | fTargetThetax = src.fTargetThetax; |
455 | fTargetThetay = src.fTargetThetay; |
456 | fTargetThetaz = src.fTargetThetaz; |
457 | for (int i=0; i < 3; ++i) |
458 | for (int j=0; j < 3; ++j) |
459 | fTargetRmatrix[i][j] = src.fTargetRmatrix[i][j]; |
460 | fBeamEnergy = src.fBeamEnergy; |
461 | fBeamErms = src.fBeamErms; |
462 | fBeamEmittance = src.fBeamEmittance; |
463 | fCollimatorSpotrms = src.fCollimatorSpotrms; |
464 | fCollimatorDistance = src.fCollimatorDistance; |
465 | fCollimatorDiameter = src.fCollimatorDiameter; |
466 | fQ2theta2 = src.fQ2theta2; |
467 | fQ2weight = src.fQ2weight; |
468 | return *this; |
469 | } |
470 | |
471 | CobremsGeneration::~CobremsGeneration() { } |
472 | |
473 | double CobremsGeneration::CoherentEnhancement(double x) |
474 | { |
475 | |
476 | |
477 | |
478 | double yc = Rate_dNcdx(x); |
479 | double yi = Rate_dNidx(x); |
480 | return (yi + yc) / (yi + 1e-99); |
481 | } |
482 | |
483 | double CobremsGeneration::Rate_dNtdx(double x) |
484 | { |
485 | |
486 | |
487 | |
488 | return Rate_dNcdx(x) + Rate_dNidx(x); |
489 | } |
490 | |
491 | double CobremsGeneration::Rate_dNtdx(double x, |
492 | double distance_m, double diameter_m) |
493 | { |
494 | |
495 | |
496 | |
497 | |
498 | |
499 | |
500 | double dist = fCollimatorDistance; |
501 | double diam = fCollimatorDiameter; |
502 | fCollimatorDistance = (distance_m > 0)? distance_m : fCollimatorDistance; |
503 | fCollimatorDiameter = (diameter_m > 0)? diameter_m : (diameter_m < 0)? |
504 | -2 * distance_m * diameter_m * me / fBeamEnergy : |
505 | fCollimatorDiameter; |
506 | double rate = Rate_dNtdx(x); |
507 | fCollimatorDistance = dist; |
508 | fCollimatorDiameter = diam; |
509 | return rate; |
510 | } |
511 | |
512 | double CobremsGeneration::Rate_dNtdk(double k_GeV) |
513 | { |
514 | |
515 | |
516 | |
517 | return Rate_dNtdx(k_GeV / fBeamEnergy) / fBeamEnergy; |
518 | } |
519 | |
520 | double CobremsGeneration::Rate_dNcdx(double x) |
521 | { |
522 | |
523 | |
524 | |
525 | double rate = 0; |
526 | int npoints = 2; |
527 | for (int n=0; n < npoints; ++n) { |
528 | double phi = (n + 0.5) * (dpi/2) / npoints; |
529 | rate += Rate_dNcdxdp(x, phi); |
530 | } |
531 | rate *= 2*dpi / npoints; |
532 | return rate; |
533 | } |
534 | |
535 | double CobremsGeneration::Rate_dNcdx(double x, |
536 | double distance_m, double diameter_m) |
537 | { |
538 | |
539 | |
540 | |
541 | |
542 | |
543 | |
544 | double dist = fCollimatorDistance; |
545 | double diam = fCollimatorDiameter; |
546 | fCollimatorDistance = (distance_m > 0)? distance_m : fCollimatorDistance; |
547 | fCollimatorDiameter = (diameter_m > 0)? diameter_m : (diameter_m < 0)? |
548 | -2 * distance_m * diameter_m * me / fBeamEnergy : |
549 | fCollimatorDiameter; |
550 | double rate = 0; |
551 | int npoints = 2; |
552 | for (int n=0; n < npoints; ++n) { |
553 | double phi = (n + 0.5) * (dpi/2) / npoints; |
554 | rate += Rate_dNcdxdp(x, phi); |
555 | } |
556 | rate *= 2*dpi / npoints; |
557 | fCollimatorDistance = dist; |
558 | fCollimatorDiameter = diam; |
559 | return rate; |
560 | } |
561 | |
562 | double CobremsGeneration::Rate_dNcdxdp(double x, double phi) |
563 | { |
564 | |
565 | |
566 | |
567 | |
568 | |
569 | |
570 | |
571 | |
572 | double Z = fTargetCrystal.Z; |
573 | double a = fTargetCrystal.lattice_constant; |
574 | double sigma0 = 16 * dpi * fTargetThickness * Z*Z * pow(alpha, 3) * |
575 | fBeamEnergy * hbarc/(a*a) * pow(hbarc / (a * me), 4); |
576 | |
577 | fQ2theta2.clear(); |
578 | fQ2weight.clear(); |
579 | double qzmin = 99; |
580 | int hmin, kmin, lmin; |
581 | double sum = 0; |
582 | |
583 | for (int h = -4; h <= 4; ++h) { |
584 | for (int k = -10; k <= 10; ++k) { |
585 | for (int l = -10; l <= 10; ++l) { |
586 | if (h/2 * 2 == h) { |
587 | if (k/2 * 2 != k || l/2 * 2 != l || |
588 | (h + k + l)/4 * 4 != h + k + l) |
589 | { |
590 | continue; |
591 | } |
592 | } |
593 | else if (k/2 * 2 == k || l/2 * 2 == l) { |
594 | continue; |
595 | } |
596 | double ReS = 0; |
597 | double ImS = 0; |
598 | for (int i=0; i < fTargetCrystal.nsites; ++i) { |
599 | double qdota = 2 * dpi * (h * fTargetCrystal.ucell_site[i].x + |
600 | k * fTargetCrystal.ucell_site[i].y + |
601 | l * fTargetCrystal.ucell_site[i].z); |
602 | ReS += cos(qdota); |
603 | ImS += sin(qdota); |
604 | } |
605 | double S2 = ReS*ReS + ImS*ImS; |
606 | if (S2 < 1e-4) |
607 | continue; |
608 | double qnorm = hbarc * 2 * dpi / a; |
609 | double q[3]; |
610 | q[0] = qnorm * (fTargetRmatrix[0][0] * h + |
611 | fTargetRmatrix[0][1] * k + |
612 | fTargetRmatrix[0][2] * l); |
613 | q[1] = qnorm * (fTargetRmatrix[1][0] * h + |
614 | fTargetRmatrix[1][1] * k + |
615 | fTargetRmatrix[1][2] * l); |
616 | q[2] = qnorm * (fTargetRmatrix[2][0] * h + |
617 | fTargetRmatrix[2][1] * k + |
618 | fTargetRmatrix[2][2] * l); |
619 | double q2 = q[0]*q[0] + q[1]*q[1] + q[2]*q[2]; |
620 | double qT2 = q[0]*q[0] + q[1]*q[1]; |
621 | double xmax = 2 * fBeamEnergy * q[2]; |
622 | xmax /= xmax + me*me; |
623 | if (x > xmax || xmax > 1) { |
624 | continue; |
625 | } |
626 | |
627 | #if COBREMS_GENERATOR_VERBOSITY1 > 2 |
628 | else { |
629 | std::cout << h << "," << k << "," << l << "," |
630 | << S2 << "," << q2 << "," << xmax |
631 | << std::endl; |
632 | } |
633 | #endif |
634 | |
635 | if (q[2] < qzmin) { |
636 | qzmin = q[2]; |
637 | hmin = h; |
638 | kmin = k; |
| Value stored to 'kmin' is never read |
639 | lmin = l; |
640 | } |
641 | double theta2 = (1 - x) * xmax / (x * (1 - xmax) + 1e-99) - 1; |
642 | double betaFF2 = pow(fTargetCrystal.betaFF, 2); |
643 | double FF = 1 / (1 + q2 * betaFF2); |
644 | sum += sigma0 * qT2 * S2 * pow(FF * betaFF2, 2) * |
645 | exp(-q2 * fTargetCrystal.Debye_Waller_const) * |
646 | ((1 - x) / pow(x * (1 + theta2) + 1e-99, 2)) * |
647 | ((1 + pow(1 - x, 2)) - 8 * (theta2 / pow(1 + theta2, 2) * |
648 | (1 - x) * pow(cos(phi), 2))) * |
649 | ((fCollimatedFlag)? Acceptance(theta2) : 1) * |
650 | ((fPolarizedFlag)? Polarization(x, theta2, phi) : 1); |
651 | fQ2theta2.push_back(theta2); |
652 | fQ2weight.push_back(sum); |
653 | } |
654 | } |
655 | } |
656 | |
657 | #if COBREMS_GENERATOR_VERBOSITY1 > 1 |
658 | if (qzmin < 99) { |
659 | #else |
660 | if (false) { |
661 | #endif |
662 | std::cout << hmin << "," << kmin << "," << lmin |
663 | << " is the best plane at x=" << x |
664 | << std::endl; |
665 | } |
666 | |
667 | return sum; |
668 | } |
669 | |
670 | double CobremsGeneration::Rate_dNidx(double x) |
671 | { |
672 | |
673 | |
674 | |
675 | if (x > 1) |
676 | return 0; |
677 | |
678 | |
679 | |
680 | int niter = 50; |
681 | double dNidx = 0; |
682 | double du = 1. / niter; |
683 | for (int iter = 0; iter < niter; ++iter) { |
684 | double u = (iter + 0.5) / niter; |
685 | double theta2 = (1 - u) / u; |
686 | dNidx += Rate_dNidxdt2(x, theta2) * du/(u*u); |
687 | } |
688 | return dNidx; |
689 | } |
690 | |
691 | double CobremsGeneration::Rate_dNBidx(double x) |
692 | { |
693 | |
694 | |
695 | |
696 | |
697 | |
698 | |
699 | |
700 | |
701 | |
702 | |
703 | |
704 | |
705 | |
706 | |
707 | |
708 | |
709 | |
710 | |
711 | |
712 | |
713 | |
714 | |
715 | |
716 | |
717 | |
718 | |
719 | |
720 | |
721 | |
722 | |
723 | double Z = fTargetCrystal.Z; |
724 | double betaFF = fTargetCrystal.betaFF; |
725 | double a = fTargetCrystal.lattice_constant; |
726 | double AoverB2 = fTargetCrystal.Debye_Waller_const / (betaFF * betaFF); |
727 | double Tfact = -(1 + AoverB2) * exp(AoverB2) * |
728 | boost::math::expint(1, AoverB2); |
729 | double psiC1 = 2 * (2 * log(betaFF * me) + Tfact + 2); |
730 | double psiC2 = psiC1 - 2/3.; |
731 | double zeta = log(1440 * pow(Z, -2/3.)) / log(183 * pow(Z, -1/3.)); |
732 | double dNBidx = fTargetCrystal.nsites * fTargetThickness * |
733 | Z * (Z + zeta) * pow(alpha, 3) * |
734 | pow(hbarc / (a*me), 2) / (a * x) * |
735 | (psiC1 * (1 + pow(1 - x, 2)) - psiC2 * (1 - x) * 2/3.); |
736 | return dNBidx; |
737 | } |
738 | |
739 | double CobremsGeneration::Rate_dNidxdt2(double x, double theta2) |
740 | { |
741 | |
742 | |
743 | |
744 | |
745 | |
746 | |
747 | |
748 | double delta = 1.02; |
749 | double Z = fTargetCrystal.Z; |
750 | double betaFF = fTargetCrystal.betaFF; |
751 | double a = fTargetCrystal.lattice_constant; |
752 | double zeta = log(1440 * pow(Z, -2/3.)) / log(183 * pow(Z, -1/3.)); |
753 | double MSchiff = 1 / (pow((me * x) / (2*fBeamEnergy * (1 - x) + 1e-99), 2) + |
754 | 1 / pow(betaFF * me * (1 + theta2), 2)); |
755 | double dNidxdt2 = 2 * fTargetCrystal.nsites * fTargetThickness * Z * |
756 | (Z + zeta) * pow(alpha, 3) * pow(hbarc/(a*me), 2) / (a*x) * |
757 | ( ((1 + pow(1 - x, 2)) - 4 * theta2 * (1 - x) / |
758 | pow(1 + theta2, 2)) / |
759 | pow(1 + theta2, 2) * |
760 | (log(MSchiff) - 2 * delta * Z / (Z + zeta)) + |
761 | 16 * theta2 * (1 - x) / pow(1 + theta2, 4) - |
762 | pow(2 - x, 2) / pow(1 + theta2, 2) ) * |
763 | ((fCollimatedFlag)? Acceptance(theta2) : 1); |
764 | return dNidxdt2; |
765 | } |
766 | |
767 | double CobremsGeneration::Rate_para(double x, double theta2, double phi) |
768 | { |
769 | |
770 | |
771 | |
772 | |
773 | |
774 | |
775 | return 0.5 * pow((2 - x) * (1 + theta2), 2) - |
776 | 8 * theta2 * (1 - x) * pow(cos(phi), 2) - |
777 | 8 * pow(theta2, 2) * (1 - x) * pow(cos(phi) * sin(phi), 2); |
778 | } |
779 | |
780 | double CobremsGeneration::Rate_ortho(double x, double theta2, double phi) |
781 | { |
782 | |
783 | |
784 | |
785 | |
786 | |
787 | |
788 | return 0.5 * pow(x * (1 + theta2), 2) + |
789 | 8 * pow(theta2, 2) * (1 - x) * pow(cos(phi) * sin(phi), 2); |
790 | } |
791 | |
792 | double CobremsGeneration::Polarization(double x, double theta2) |
793 | { |
794 | |
795 | |
796 | |
797 | |
798 | |
799 | |
800 | |
801 | return 2 * (1 - x) / (pow(1 + theta2, 2) * (pow(1 - x + 1e-99, 2) + 1) - |
802 | 4 * theta2 * (1 - x)); |
803 | } |
804 | |
805 | double CobremsGeneration::Polarization(double x, double theta2, double phi) |
806 | { |
807 | |
808 | |
809 | |
810 | |
811 | |
812 | double Rpara = Rate_para(x, theta2, phi); |
813 | double Rperp = Rate_ortho(x, theta2, phi); |
814 | return (Rpara - Rperp) / (Rpara + Rperp); |
815 | } |
816 | |
817 | double CobremsGeneration::AbremsPolarization(double x, double theta2, double phi) |
818 | { |
819 | |
820 | |
821 | |
822 | |
823 | |
824 | |
825 | |
826 | double Acoeff[3][4] = {{0.93000, 0.64250, 0.66598, 1.62506}, |
827 | {0.73000, 1.05648, 0.84643, 1.97061}, |
828 | {0.87610, 0.57510, 0.74918, 1.52849}}; |
829 | double a[3]; |
830 | for (int n=0; n < 3; ++n) { |
831 | double A = pow(Acoeff[n][0], 2) + |
832 | pow(Acoeff[n][1], 2) * pow(x, 2) + |
833 | pow(Acoeff[n][2], 2) * pow(x, 4) + |
834 | pow(Acoeff[n][3], 2) * pow(x, 16); |
835 | a[n] = A*A; |
836 | } |
837 | double ppol = theta2 / (a[0] + a[1] * theta2 + a[2] * theta2*theta2); |
838 | return ppol * cos(2 * phi); |
839 | } |
840 | |
841 | double CobremsGeneration::Acceptance(double theta2, double phi, |
842 | double xshift_m, double yshift_m) |
843 | { |
844 | |
845 | |
846 | |
847 | |
848 | |
849 | |
850 | |
851 | |
852 | |
853 | double theta = sqrt(theta2) * (me/fBeamEnergy); |
854 | double xc = fCollimatorDistance * tan(theta) * cos(phi) + xshift_m; |
855 | double yc = fCollimatorDistance * tan(theta) * sin(phi) + yshift_m; |
856 | double thetaprime = atan2(sqrt(xc*xc + yc*yc), fCollimatorDistance); |
857 | return Acceptance(pow(thetaprime * fBeamEnergy/me, 2)); |
858 | } |
859 | |
860 | double CobremsGeneration::Acceptance(double theta2) |
861 | { |
862 | |
863 | |
864 | |
865 | |
866 | |
867 | |
868 | |
869 | |
870 | |
871 | double acceptance = 0; |
872 | double niter = 50; |
873 | double theta = sqrt(theta2); |
874 | double thetaC = fCollimatorDiameter / (2 * fCollimatorDistance) * |
875 | fBeamEnergy / me; |
876 | double var0 = pow((fCollimatorSpotrms / fCollimatorDistance) * |
877 | fBeamEnergy / me, 2); |
878 | double varMS = Sigma2MS(fTargetThickness) * pow(fBeamEnergy / me, 2); |
879 | if (theta < thetaC) { |
880 | double u1 = thetaC - theta; |
881 | if (u1*u1 / (var0 + varMS) > 20) { |
882 | return 1; |
883 | } |
884 | for (int iter = 0; iter < niter; ++iter) { |
885 | double u = u1 * (iter + 0.5) / niter; |
886 | double u2 = u * u; |
887 | double du2 = 2 * u * u1 / niter; |
888 | double pu; |
889 | if (varMS / var0 > 1e-4) { |
890 | pu = (boost::math::expint(1, u2 / (2 * (var0 + varMS))) - |
891 | boost::math::expint(1, u2 / (2 * var0))) / (2 * varMS); |
892 | } |
893 | else { |
894 | pu = exp(-u2 / (2 * var0)) / (2 * var0); |
895 | } |
896 | acceptance += pu * du2; |
897 | } |
898 | } |
899 | double u0 = fabs(theta - thetaC); |
900 | double u1 = fabs(theta + thetaC); |
901 | for (int iter = 0; iter < niter; ++iter) { |
902 | double u = u0 + (u1 - u0) * (iter + 0.5) / niter; |
903 | double u2 = u * u; |
904 | double du2 = 2 * u * (u1 - u0) / niter; |
905 | double pu; |
906 | if (varMS / var0 > 1e-4) { |
907 | pu = (boost::math::expint(1, u2 / (2 * (var0 + varMS))) - |
908 | boost::math::expint(1, u2 / (2 * var0))) / (2 * varMS); |
909 | } |
910 | else { |
911 | pu = exp(-u2 / (2 * var0)) / (2 * var0); |
912 | } |
913 | acceptance += pu * du2/dpi * |
914 | atan2(sqrt((theta2 - pow(thetaC - u, 2)) * |
915 | (pow(thetaC + u, 2) - theta2)), |
916 | theta2 - pow(thetaC, 2) + u2); |
917 | } |
918 | return acceptance; |
919 | } |
920 | |
921 | void CobremsGeneration::RotateTarget(double thetax, |
922 | double thetay, |
923 | double thetaz) |
924 | { |
925 | |
926 | |
927 | |
928 | |
929 | if (thetaz != 0) { |
930 | double sint = sin(thetaz); |
931 | double cost = cos(thetaz); |
932 | for (int i=0; i < 3; ++i) { |
933 | double x = fTargetRmatrix[0][i]; |
934 | double y = fTargetRmatrix[1][i]; |
935 | fTargetRmatrix[0][i] = cost * x + sint * y; |
936 | fTargetRmatrix[1][i] = cost * y - sint * x; |
937 | } |
938 | } |
939 | if (thetay != 0) { |
940 | double sint = -sin(thetay); |
941 | double cost = cos(thetay); |
942 | for (int i=0; i < 3; ++i) { |
943 | double x = fTargetRmatrix[0][i]; |
944 | double z = fTargetRmatrix[2][i]; |
945 | fTargetRmatrix[0][i] = cost * x + sint * z; |
946 | fTargetRmatrix[2][i] = cost * z - sint * x; |
947 | } |
948 | } |
949 | if (thetax != 0) { |
950 | double sint = sin(thetax); |
951 | double cost = cos(thetax); |
952 | for (int i=0; i < 3; ++i) { |
953 | double y = fTargetRmatrix[1][i]; |
954 | double z = fTargetRmatrix[2][i]; |
955 | fTargetRmatrix[1][i] = cost * y + sint * z; |
956 | fTargetRmatrix[2][i] = cost * z - sint * y; |
957 | } |
958 | } |
959 | } |
960 | |
961 | double CobremsGeneration::Sigma2MS(double thickness_m) |
962 | { |
963 | |
964 | |
965 | |
966 | |
967 | |
968 | |
969 | |
970 | return fabs(Sigma2MS_Geant(thickness_m)); |
971 | } |
972 | |
973 | double CobremsGeneration::Sigma2MS_Kaune(double thickness_m) |
974 | { |
975 | |
976 | |
977 | |
978 | |
979 | |
980 | |
981 | |
982 | |
983 | |
984 | |
985 | |
986 | |
987 | |
988 | |
989 | |
990 | |
991 | |
992 | |
993 | double carboncor = 4.2 / 4.6; |
994 | double Z = fTargetCrystal.Z; |
995 | double a = fTargetCrystal.lattice_constant; |
996 | return 8 * dpi * fTargetCrystal.nsites * pow(alpha * Z, 2) * |
997 | thickness_m * pow(hbarc / (fBeamEnergy * a), 2) / a * |
998 | log(183 * pow(Z, -1/3.)) * |
999 | carboncor; |
1000 | } |
1001 | |
1002 | double CobremsGeneration::Sigma2MS_PDG(double thickness_m) |
1003 | { |
1004 | |
1005 | |
1006 | |
1007 | |
1008 | double t = thickness_m / fTargetCrystal.radiation_length; |
1009 | return pow(13.6e-3 / fBeamEnergy, 2) * t * pow(1 + 0.038 * log(t), 2); |
1010 | } |
1011 | |
1012 | double CobremsGeneration::Sigma2MS_Geant(double thickness_m) |
1013 | { |
1014 | |
1015 | |
1016 | |
1017 | |
1018 | |
1019 | |
1020 | |
1021 | |
1022 | |
1023 | double rBohr = 0.52917721e-10; |
1024 | double F = 0.98; |
1025 | double Z = fTargetCrystal.Z; |
1026 | double chi2cc = pow(0.39612e-2, 2) * Z * (Z + 1) * |
1027 | fTargetCrystal.density / 12; |
1028 | double chi2c = chi2cc * thickness_m / pow(fBeamEnergy, 2); |
1029 | double chi2alpha = 1.13 * pow(hbarc / (fBeamEnergy * rBohr * 0.885), 2) * |
1030 | pow(Z, 2/3.) * (1 + 3.34 * pow(alpha * Z, 2)); |
1031 | double omega0 = chi2c / (1.167 * chi2alpha); |
1032 | double gnu = omega0 / (2 * (1 - F)); |
1033 | return chi2c / (1 + pow(F, 2)) * ((1 + gnu) / gnu * log(1 + gnu) -1); |
1034 | } |
1035 | |
1036 | double CobremsGeneration::Sigma2MS_Hanson(double thickness_m) |
1037 | { |
1038 | |
1039 | |
1040 | |
1041 | |
1042 | |
1043 | |
1044 | |
1045 | |
1046 | |
1047 | double Z = fTargetCrystal.Z; |
1048 | double ttingcm2 = thickness_m * 100 * fTargetCrystal.density; |
1049 | double EinMeV = fBeamEnergy * 1000; |
1050 | double theta2max = 0.157 * Z * (Z + 1) / fTargetCrystal.A * |
1051 | ttingcm2 / pow(EinMeV, 2); |
1052 | double theta2screen = theta2max * fTargetCrystal.A * |
1053 | (1 + 3.35 * pow(Z * alpha, 2)) / |
1054 | (7800 * (Z + 1) * pow(Z, 1/3.) * ttingcm2); |
1055 | double BminuslogB = log(theta2max / theta2screen) - 0.154; |
1056 | double Blast = 1; |
1057 | double B; |
1058 | for (int i=0; i < 999; ++i) { |
1059 | B = BminuslogB + log(Blast); |
1060 | if (B < 1.2) { |
1061 | B = 1.21; |
1062 | break; |
1063 | } |
1064 | else if (fabs(B - Blast) > 1e-6) { |
1065 | Blast = B; |
1066 | } |
1067 | else { |
1068 | break; |
1069 | } |
1070 | } |
1071 | return theta2max * (B - 1.2) / 2; |
1072 | } |
1073 | |
1074 | #ifdef BOOST_PYTHON_WRAPPING |
1075 | |
1076 | void CobremsGeneration::pyApplyBeamCrystalConvolution(int nbins, pyobject xarr, |
1077 | pyobject yarr) |
1078 | { |
1079 | using boost::python::extract; |
1080 | typedef boost::python::tuple pytuple; |
1081 | pytuple xtuple = extract<pytuple>(xarr.attr("buffer_info")()); |
1082 | pytuple ytuple = extract<pytuple>(yarr.attr("buffer_info")()); |
1083 | double *xbuf = reinterpret_cast<double*>((int)extract<int>(xtuple[0])); |
1084 | double *ybuf = reinterpret_cast<double*>((int)extract<int>(ytuple[0])); |
1085 | applyBeamCrystalConvolution(nbins, xbuf, ybuf); |
1086 | } |
1087 | |
1088 | double (CobremsGeneration::*Rate_dNtdx_1)(double) = &CobremsGeneration::Rate_dNtdx; |
1089 | double (CobremsGeneration::*Rate_dNtdx_3)(double, double, double) = &CobremsGeneration::Rate_dNtdx; |
1090 | double (CobremsGeneration::*Rate_dNcdx_1)(double) = &CobremsGeneration::Rate_dNcdx; |
1091 | double (CobremsGeneration::*Rate_dNcdx_3)(double, double, double) = &CobremsGeneration::Rate_dNcdx; |
1092 | double (CobremsGeneration::*Acceptance_1)(double) = &CobremsGeneration::Acceptance; |
1093 | double (CobremsGeneration::*Acceptance_4)(double, double, double, double) = &CobremsGeneration::Acceptance; |
1094 | double (CobremsGeneration::*Polarization_2)(double, double) = &CobremsGeneration::Polarization; |
1095 | double (CobremsGeneration::*Polarization_3)(double, double, double) = &CobremsGeneration::Polarization; |
1096 | |
1097 | BOOST_PYTHON_MODULE(libcobrems) |
1098 | { |
1099 | using boost::python::class_; |
1100 | using boost::python::enum_; |
1101 | using boost::python::def; |
1102 | |
1103 | class_<CobremsGeneration, CobremsGeneration*> |
1104 | ("CobremsGeneration", |
1105 | "coherent bremsstrahlung spectrum and polarization calculator, " |
1106 | "with methods for generating random Monte Carlo samples", |
1107 | boost::python::init<double, double>()) |
1108 | .def("setBeamEnergy", &CobremsGeneration::setBeamEnergy) |
1109 | .def("setBeamErms", &CobremsGeneration::setBeamErms) |
1110 | .def("setBeamEmittance", &CobremsGeneration::setBeamEmittance) |
1111 | .def("setCollimatorSpotrms", &CobremsGeneration::setCollimatorSpotrms) |
1112 | .def("setCollimatorDistance", &CobremsGeneration::setCollimatorDistance) |
1113 | .def("setCollimatorDiameter", &CobremsGeneration::setCollimatorDiameter) |
1114 | .def("setTargetThickness", &CobremsGeneration::setTargetThickness) |
1115 | .def("setTargetCrystal", &CobremsGeneration::setTargetCrystal) |
1116 | .def("setCoherentEdge", &CobremsGeneration::setCoherentEdge) |
1117 | .def("setTargetThetax", &CobremsGeneration::setTargetThetax) |
1118 | .def("setTargetThetay", &CobremsGeneration::setTargetThetay) |
1119 | .def("setTargetThetaz", &CobremsGeneration::setTargetThetaz) |
1120 | .def("setTargetOrientation", &CobremsGeneration::setTargetOrientation) |
1121 | .def("RotateTarget", &CobremsGeneration::RotateTarget) |
1122 | .def("getBeamEnergy", &CobremsGeneration::getBeamEnergy) |
1123 | .def("getBeamErms", &CobremsGeneration::getBeamErms) |
1124 | .def("getBeamEmittance", &CobremsGeneration::getBeamEmittance) |
1125 | .def("getCollimatorSpotrms", &CobremsGeneration::getCollimatorSpotrms) |
1126 | .def("getCollimatorDistance", &CobremsGeneration::getCollimatorDistance) |
1127 | .def("getCollimatorDiameter", &CobremsGeneration::getCollimatorDiameter) |
1128 | .def("getTargetThickness", &CobremsGeneration::getTargetThickness) |
1129 | .def("getTargetCrystal", &CobremsGeneration::getTargetCrystal) |
1130 | .def("getTargetCrystalNsites", &CobremsGeneration::getTargetCrystalNsites) |
1131 | .def("getTargetCrystalAtomicNumber", &CobremsGeneration::getTargetCrystalAtomicNumber) |
1132 | .def("getTargetCrystalAtomicWeight", &CobremsGeneration::getTargetCrystalAtomicWeight) |
1133 | .def("getTargetCrystalDensity", &CobremsGeneration::getTargetCrystalDensity) |
1134 | .def("getTargetCrystalLatticeConstant", &CobremsGeneration::getTargetCrystalLatticeConstant) |
1135 | .def("getTargetCrystalRadiationLength", &CobremsGeneration::getTargetCrystalRadiationLength) |
1136 | .def("getTargetCrystalDebyeWallerConst", &CobremsGeneration::getTargetCrystalDebyeWallerConst) |
1137 | .def("getTargetCrystalMosaicSpread", &CobremsGeneration::getTargetCrystalMosaicSpread) |
1138 | .def("getTargetCrystalBetaFF", &CobremsGeneration::getTargetCrystalBetaFF) |
1139 | .def("getTargetThetax", &CobremsGeneration::getTargetThetax) |
1140 | .def("getTargetThetay", &CobremsGeneration::getTargetThetay) |
1141 | .def("getTargetThetaz", &CobremsGeneration::getTargetThetaz) |
1142 | .def("getTargetRadiationLength_PDG", &CobremsGeneration::getTargetRadiationLength_PDG) |
1143 | .def("getTargetRadiationLength_Schiff", &CobremsGeneration::getTargetRadiationLength_Schiff) |
1144 | .def("getTargetDebyeWallerConstant", &CobremsGeneration::getTargetDebyeWallerConstant) |
1145 | .def("getCollimatedFlag", &CobremsGeneration::getCollimatedFlag) |
1146 | .def("setCollimatedFlag", &CobremsGeneration::setCollimatedFlag) |
1147 | .def("getPolarizedFlag", &CobremsGeneration::getPolarizedFlag) |
1148 | .def("setPolarizedFlag", &CobremsGeneration::setPolarizedFlag) |
1149 | .def("applyBeamCrystalConvolution", &CobremsGeneration::pyApplyBeamCrystalConvolution) |
1150 | .def("printBeamlineInfo", &CobremsGeneration::printBeamlineInfo) |
1151 | .def("printTargetCrystalInfo", &CobremsGeneration::printTargetCrystalInfo) |
1152 | .def("CoherentEnhancement", &CobremsGeneration::CoherentEnhancement) |
1153 | .def("Rate_dNtdx", Rate_dNtdx_1) |
1154 | .def("Rate_dNtdx", Rate_dNtdx_3) |
1155 | .def("Rate_dNtdk", &CobremsGeneration::Rate_dNtdk) |
1156 | .def("Rate_dNcdx", Rate_dNcdx_1) |
1157 | .def("Rate_dNcdx", Rate_dNcdx_3) |
1158 | .def("Rate_dNcdxdp", &CobremsGeneration::Rate_dNcdxdp) |
1159 | .def("Rate_dNidx", &CobremsGeneration::Rate_dNidx) |
1160 | .def("Rate_dNBidx", &CobremsGeneration::Rate_dNBidx) |
1161 | .def("Rate_dNidxdt2", &CobremsGeneration::Rate_dNidxdt2) |
1162 | .def("Rate_para", &CobremsGeneration::Rate_para) |
1163 | .def("Rate_ortho", &CobremsGeneration::Rate_ortho) |
1164 | .def("Polarization", Polarization_2) |
1165 | .def("Polarization", Polarization_3) |
1166 | .def("Acceptance", Acceptance_1) |
1167 | .def("Acceptance", Acceptance_4) |
1168 | .def("Sigma2MS", &CobremsGeneration::Sigma2MS) |
1169 | .def("Sigma2MS_Kaune", &CobremsGeneration::Sigma2MS_Kaune) |
1170 | .def("Sigma2MS_PDG", &CobremsGeneration::Sigma2MS_PDG) |
1171 | .def("Sigma2MS_Geant", &CobremsGeneration::Sigma2MS_Geant) |
1172 | .def("Sigma2MS_Hanson", &CobremsGeneration::Sigma2MS_Hanson) |
1173 | .def_readonly("dpi", &CobremsGeneration::dpi) |
1174 | .def_readonly("me", &CobremsGeneration::me) |
1175 | .def_readonly("alpha", &CobremsGeneration::alpha) |
1176 | .def_readonly("hbarc", &CobremsGeneration::hbarc) |
1177 | ; |
1178 | } |
1179 | |
1180 | #endif |