Bug Summary

File:libraries/UTILITIES/CobremsGeneration.cc
Location:line 639, column 16
Description:Value stored to 'lmin' is never read

Annotated Source Code

1//
2// CobremsGeneration - class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: july 27, 2015
6//
7// notes:
8//
9// This class computes the spectrum of bremsstrahlung radiation from a
10// crystal radiator. The formalism is that described in the following paper.
11//
12// W. Kaune, G. Miller, W. Oliver, R.W. Williams, and K.K. Young,
13// "Inclusive cross sections for pion and proton production by photons
14// using collimated coherent bremsstrahlung", Phys Rev D, vol 11,
15// no 3 (1975) pp. 478-494.
16//
17// The model for the photon beam contains the following parameters.
18// 1. electron beam
19// * beam energy: mean and rms spread
20// * spot on radiator: gaussian model, cylindrical symmetry
21// * emittance: gaussian model, cylindrical symmetry
22// 2. crystal target
23// * implemented for diamond, silicon
24// * uniform thickness across beam spot
25// * mosaic spread: gaussian model
26// * dipole atomic form factor
27// * Debye-Waller factor: defines coherent domain in q sum
28// 3. downstream collimator
29// * fixed distance from radiator
30// * sharp cutoff at collimator radius
31// * perfect alignment with beam axis assumed
32//
33// The crystal orientation is computed based on the requested coherent
34// edge position requested by the user. For a high-energy electron beam
35// this fixes one of the angles of the crystal with respect to the
36// electron beam. The other angle must be chosen based on other
37// considerations. A default value for this secondary angle parameter
38// is assigned below, based on the observation that
39// a) it is significantly larger than the primary edge-defining
40// angle, so that additional peaks from reciprocal lattice sites
41// (hkl) with the same h and different k values do not contribute
42// significantly to the spectrum below the endpoint, and
43// b) it is small enough, to render it unlikely that a random lattice
44// vector from a distant region in q-space will cross through the
45// coherent enhancement region as the primary (220) peak is moved
46// through its full range in x from 30% to 90% of the endpoint
47// for beamline parameters similar to those describing Hall D and GlueX.
48// Should the user wish to try other values for this angle, a public
49// method is provided for this purpose.
50
51#define COBREMS_GENERATOR_VERBOSITY1 1
52//#define BOOST_PYTHON_WRAPPING 1
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
60const double CobremsGeneration::dpi = 3.1415926535897;
61const double CobremsGeneration::me = 0.510998910e-3;
62const double CobremsGeneration::alpha = 7.2973525698e-3;
63const double CobremsGeneration::hbarc = 0.1973269718e-15;
64
65CobremsGeneration::CobremsGeneration(double Emax_GeV, double Epeak_GeV)
66{
67 // Unique constructor for this class, initialize for the given
68 // endpoint energy and peak position but these can be changed.
69
70 fBeamEnergy = Emax_GeV;
71 fBeamErms = 6.0e-4; // GeV
72 fBeamEmittance = 2.5e-9; // m r
73 fCollimatorSpotrms = 0.0005; // m
74 fCollimatorDistance = 76.0; // m
75 fCollimatorDiameter = 0.005; // m
76 fTargetThickness = 50e-6; // m
77 fTargetThetay = 0.050; // radians
78 fTargetThetaz = 0; // radians
79 setTargetCrystal("diamond");
80 setCoherentEdge(Epeak_GeV);
81 fPhotonEnergyMin = 0.120; // GeV
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
96void CobremsGeneration::updateTargetOrientation()
97{
98 resetTargetOrientation();
99 RotateTarget(0, dpi/2, 0); // point (1,0,0) along beam
100 RotateTarget(0, 0, dpi/4); // point (0,1,1) vertically
101 RotateTarget(0, 0, -fTargetThetaz);
102 RotateTarget(0, -fTargetThetay, 0);
103 RotateTarget(-fTargetThetax, 0, 0);
104}
105
106void CobremsGeneration::setTargetCrystal(std::string crystal)
107{
108 // declare the radiator target crystal type by name
109
110 if (crystal == "diamond") {
111 fTargetCrystal.name = "diamond";
112 fTargetCrystal.Z = 6;
113 fTargetCrystal.A = 12.01;
114 fTargetCrystal.density = 3.534; // g/cm^3
115 fTargetCrystal.lattice_constant = 3.5668e-10; // m
116 fTargetCrystal.Debye_Waller_const = 0.40e9; // 1/GeV^2
117 }
118 else if (crystal == "silicon") {
119 fTargetCrystal.name = "silicon";
120 fTargetCrystal.Z = 14;
121 fTargetCrystal.A = 28.09;
122 fTargetCrystal.density = 2.320; // g/cm^3
123 fTargetCrystal.lattice_constant = 5.431e-10; // m
124 fTargetCrystal.Debye_Waller_const = 1.5e9; // 1/GeV^2
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 // define the stanard unit cell of the diamond Bravais lattice
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 // approximate formula for atomic form factor beta
147 fTargetCrystal.betaFF = 111 * pow(fTargetCrystal.Z, -1/3.) / me;
148
149 // set mosaic spread to GlueX specification
150 fTargetCrystal.mosaic_spread = 20e-6;
151
152 // compute the radiation length
153 fTargetCrystal.radiation_length = getTargetRadiationLength_Schiff();
154}
155
156double CobremsGeneration::getTargetDebyeWallerConstant(double DebyeT_K,
157 double T_K)
158{
159 // Computes the Debye-Waller constant A for a simple model
160 // assuming an isotropic crystal -- see Kaune et.al.
161 //
162 // A(T) = A0 f(T)
163 // where
164 // A0 = 3 / (4 * atomicMass_GeV * DebyeTemperature_GeV)
165 // and
166 // f(T) = (2 / DebyeTemperature_GeV^2) *
167 // Integral_dw[0,DebyeTemperature_GeV]
168 // {w (1 + 2 / (exp(w/T) - 1))}
169 //
170 // T is the crystal temperature in GeV and A0 is the limiting
171 // value of A as T->0.
172
173 double kBoltzmann = 8.617e-14; // GeV/K
174 double amassGeV = fTargetCrystal.A * 0.932; // GeV
175 double A0 = 3 / (2 * amassGeV * kBoltzmann * DebyeT_K); // /GeV^2
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
186void CobremsGeneration::printBeamlineInfo()
187{
188 // Print a summary of the target crystal model parameters
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
218void CobremsGeneration::printTargetCrystalInfo()
219{
220 // Print a summary of the target crystal model parameters
221
222 double kBoltzmann = 8.617e-14; // GeV/K
223 double amass = fTargetCrystal.A * 0.932; // GeV
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
274void CobremsGeneration::applyBeamCrystalConvolution(int nbins, double *xvalues,
275 double *yvalues)
276{
277 // Electron beam emittance produces two effects in the coherent
278 // bremsstrahlung spectrum:
279 // 1) smears out the collimation acceptance function in production
280 // angle, so it varies smoothly to zero instead of being a step;
281 // 2) combines with mosaic spread of the target crystal to smear out
282 // the relation between photon energy fraction x and production
283 // angle theta for a given lattice reflection.
284 // The first one affects the left-hand (low energy) side of the
285 // coherent peaks in the coherent bremsstrahlung spectrum, while the
286 // second affects the right-hand (high energy) side, limiting the
287 // sharpness of the edges in either case.
288 //
289 // Effect (1) is taken into account in the way the acceptance function
290 // is computed on the final state, but effect (2) is more difficult to
291 // treat analytically because it involves smearing directions in the
292 // initial state. This is only relevant to the coherent part of the
293 // spectrum, where it acts by broadening the step on the high side of
294 // the coherent edge from a sharp drop to a gradually sloping curve.
295 // One way to take this into account in an effective manner is to treat
296 // the coherent spectrum at every photon energy bin as being dominated
297 // by a single reciprocal lattice vector, and smearing out the relation
298 // x=x(theta) that follows from the two-body nature of the scattering
299 // from that plane. Considering a 1D spectrum of beam intensity vs x,
300 // this leads to a convolution of the distribution with an x-dependent
301 // smearing function. The applyBeamConvolution method computes that
302 // smearing function for each value of x in the input xvalues array
303 // and applies it to the input spectrum represented by the yvalues
304 // array. The yvalues array is overwritten with the convoluted spectrum.
305 // For simplicity, the xvalues are assumed to be equally spaced.
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 // Here we have to guess which reciprocal lattice vector is dominantly
314 // for the coherent photons in each bin in x. For simplicity, I assume
315 // it is a (2,2,0) vector. Higher order vectors exhibit more smearing
316 // but this is a good approximation if the primary peaks in the spectrum
317 // come from (2,2,0) vectors.
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
382double CobremsGeneration::getTargetRadiationLength_PDG()
383{
384 // PDG formula for radiation length, converted to meters
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
398double CobremsGeneration::getTargetRadiationLength_Schiff()
399{
400 // Schiff formula for radiation length, converted to meters
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
411void CobremsGeneration::setCoherentEdge(double Epeak_GeV)
412{
413 // Adjust theta_x of the target to align the coherent edge at
414 // energy Epeak_GeV in the photon spectrum, then orient the
415 // crystal according to theta_x, theta_y, theta_z tip angles.
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
426CobremsGeneration::CobremsGeneration(const CobremsGeneration &src)
427{
428 // copy constructor
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
448CobremsGeneration &CobremsGeneration::operator=(const CobremsGeneration &src)
449{
450 // assignment operator
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
471CobremsGeneration::~CobremsGeneration() { }
472
473double CobremsGeneration::CoherentEnhancement(double x)
474{
475 // Returns ratio of total bremsstrahlung yield over incoherent yield
476 // for photon energy k = x*fBeamEnergy
477
478 double yc = Rate_dNcdx(x);
479 double yi = Rate_dNidx(x);
480 return (yi + yc) / (yi + 1e-99);
481}
482
483double CobremsGeneration::Rate_dNtdx(double x)
484{
485 // Returns total bremsstrahlung probability density differential in
486 // x (scaled photon energy) at photon energy k = x*fBeamEnergy.
487
488 return Rate_dNcdx(x) + Rate_dNidx(x);
489}
490
491double CobremsGeneration::Rate_dNtdx(double x,
492 double distance_m, double diameter_m)
493{
494 // Returns total bremsstrahlung probability density differential in x
495 // (scaled photon energy) at photon energy k = x*fBeamEnergy with
496 // user-specified variations in the collimator distance and diameter,
497 // for plotting. Special case: if diameter_m < 0 then interpret its
498 // absolute value as the collimator radius in characteristic units m/E.
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
512double CobremsGeneration::Rate_dNtdk(double k_GeV)
513{
514 // Returns total bremsstrahlung probability density differential
515 // in photon energy k (GeV).
516
517 return Rate_dNtdx(k_GeV / fBeamEnergy) / fBeamEnergy;
518}
519
520double CobremsGeneration::Rate_dNcdx(double x)
521{
522 // Returns the coherent bremsstrahlung probability density differential
523 // in x (scaled photon energy) at photon energy k = x*fBeamEnergy.
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
535double CobremsGeneration::Rate_dNcdx(double x,
536 double distance_m, double diameter_m)
537{
538 // Returns the coherent bremsstrahlung probability density differential
539 // in x (scaled photon energy) at photon energy k = x*fBeamEnergy with
540 // user-specified variations in the collimator distance and diameter.
541 // Special case: if diameter_m < 0 then interpret its absolute
542 // value as the collimator radius in characteristic units m/E.
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
562double CobremsGeneration::Rate_dNcdxdp(double x, double phi)
563{
564 // Returns the coherent bremsstrahlung probabililty density differential
565 // in x (scaled photon energy) and phi (azimuthal emission angle) for
566 // fixed photon energy k = x*fBeamEnergy and phi. If fPolarizedFlag is
567 // false (0, default) then the total yield is returned, otherwise it is
568 // only the polarized fraction. If fCollimatedFlag is false (0) then
569 // the total yield is returned, otherwise only the part that passes the
570 // collimator is counted (default).
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 // can restrict to h=0 for cpu speedup, if crystal alignment is "reasonable"
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;
639 lmin = l;
Value stored to 'lmin' is never read
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
670double CobremsGeneration::Rate_dNidx(double x)
671{
672 // Returns the incoherent bremsstrahlung probabililty density differential
673 // in x (scaled photon energy) at fixed photon energy k = x*fBeamEnergy.
674
675 if (x > 1)
676 return 0;
677
678 // Numerical integration in d(theta**2) over [0,inf]
679 // is mapped onto u=1/(1+theta^2) as (1/u^2) d(u) over [0,1]
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
691double CobremsGeneration::Rate_dNBidx(double x)
692{
693 // In the following paper, a closed form is given for the integral that
694 // is being performed analytically by dNidx. I include this second form
695 // here in case some time it might be useful as a cross check.
696 //
697 // "Coherent bremsstrahlung in crystals as a tool for producing high
698 // energy photon beams to be used in photoproduction experiments at
699 // CERN SPS", Nucl. Instr. Meth. 204 (1983) pp.299-310.
700 //
701 // Note: in this paper they have swapped subscripts for coherent and
702 // incoherent intensities. This is not very helpful to the reader!
703 //
704 // The result is some 15% lower radiation rate than the result of dNidx.
705 // I take the latter to be more detailed (because it gives a more
706 // realistic behaviour at the endpoint and agrees better with the PDG
707 // radiation length for carbon). Most of this deficiency is remedied
708 // by simply replacing Z**2 in the cross section with Z*(Z+zeta) as
709 // recommended by Kaune et.al., and followed by the PDG in their fit
710 // to radiation lengths.
711 //
712 // WARNING
713 // dNidx and dNBidx give the incoherent radiation rate for crystalline
714 // radiators. If you take the incoherent radiation formulae here and
715 // integrate them you will NOT obtain the radiation length for amorphous
716 // radiators; it will be overestimated by some 15%. The reason is that
717 // the part of the integral in q-space that is covered by the discrete
718 // sum has been subtracted to avoid double-counting with the coherent
719 // part. If you were to spin the crystal fast enough, the coherent
720 // spectrum should average out to yield the remaining 15% with a
721 // spectral shape resembling the Bethe-Heitler result.
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
739double CobremsGeneration::Rate_dNidxdt2(double x, double theta2)
740{
741 // Returns the incoherent bremsstrahlung probabililty density differential
742 // in x (scaled photon energy) and theta^2 at fixed photon energy
743 // k = x*fBeamEnergy and production angle theta. Argument theta2 is equal
744 // to theta^2 expressed in units of (me/fBeamEnergy)^2. If internal flag
745 // fCollimatedFlag is false (0) then the total yield is returned,
746 // otherwise only the part that passes the collimator is counted (default).
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
767double CobremsGeneration::Rate_para(double x, double theta2, double phi)
768{
769 // Returns the relative rate of in-plane polarized flux from coherent
770 // bremsstrahlung at production angles theta and phi and photon energy
771 // k = x*fBeamEnergy. The units are arbitrary, but the same as Rate_ortho
772 // (see below). The argument theta2 is the production polar angle theta^2
773 // expressed in units of (me/fBeamEnergy)^2.
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
780double CobremsGeneration::Rate_ortho(double x, double theta2, double phi)
781{
782 // Returns the relative rate of out-of-plane polarized flux from coherent
783 // bremsstrahlung at production angles theta and phi and photon energy k
784 // = x*fBeamEnergy. The units are arbitrary, but the same as Rate_para
785 // (see above). The argument theta2 is the production polar angle theta^2
786 // expressed in units of (me/fBeamEnergy)^2.
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
792double CobremsGeneration::Polarization(double x, double theta2)
793{
794 // Returns the degree of linear polarization in a coherent bremsstrahlung
795 // beam at photon energy k = x*fBeamEnergy and production angle theta.
796 // The formula evaluated below is the azimuthal average of the ratio
797 // (Rate_para - Rate_ortho) / (Rate_para + Rate_ortho)
798 // The argument theta2 is the production polar angle theta^2 expressed
799 // in units of (me/fBeamEnergy)^2.
800
801 return 2 * (1 - x) / (pow(1 + theta2, 2) * (pow(1 - x + 1e-99, 2) + 1) -
802 4 * theta2 * (1 - x));
803}
804
805double CobremsGeneration::Polarization(double x, double theta2, double phi)
806{
807 // Returns the degree of linear polarization in a coherent bremsstrahlung
808 // beam at photon energy k = x*fBeamEnergy and production angles theta, phi.
809 // The argument theta2 is the production polar angle theta^2 expressed
810 // in units of (me/fBeamEnergy)^2.
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
817double CobremsGeneration::AbremsPolarization(double x, double theta2, double phi)
818{
819 // Returns the degree of linear polarization in an ordinary atomic
820 // bremsstrahlung beam at photon energy k = x*fBeamEnergy and production
821 // angles theta,phi. The formula is a parameterization of the linear
822 // polarization evaluated using the Dirac++ QED Monte Carlo generator.
823 // The argument theta2 is the production polar angle theta^2 expressed
824 // in units of (me/fBeamEnergy)^2.
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
841double CobremsGeneration::Acceptance(double theta2, double phi,
842 double xshift_m, double yshift_m)
843{
844 // Returns the acceptance of the collimator for photons emitted at
845 // polar angle theta and azimuthal angle phi at the radiator. Both
846 // beam emittance and multiple-scattering in the target contribute
847 // to smearing of the angular acceptance at the the collimator edge.
848 // The argument theta2 is the production polar angle theta^2
849 // expressed in units of (me/fBeamEnergy)^2. Misalignment of the
850 // collimator with the beam axis is taken into account by the
851 // arguments xshift,yshift.
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
860double CobremsGeneration::Acceptance(double theta2)
861{
862 // Returns the acceptance of the collimator for photons emitted at
863 // polar angle theta at the radiator, under the assumption that the
864 // collimator axis is perfectly aligned with the incident electron
865 // beam axis back at the radiator. Both beam emittance and
866 // multiple-scattering in the target contribute to smearing of the
867 // angular acceptance at the the collimator edge. The argument theta2
868 // is the production polar angle theta^2 expressed in units of
869 // (me/fBeamEnergy)^2.
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
921void CobremsGeneration::RotateTarget(double thetax,
922 double thetay,
923 double thetaz)
924{
925 // Apply a sequence of rotations to the target crystal as
926 // Rmatrix(out) = Rx(thx) Ry(thy) Rz(thz) Rmatrix(in)
927 // with rotations understood in the passive sense.
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
961double CobremsGeneration::Sigma2MS(double thickness_m)
962{
963 // Returns the mean-square multiple-scattering angle of the
964 // electron beam inside the radiator crystal target, in radians.
965 // This method wraps one of the concrete implementations, see below.
966 // Some formulas, although valid for a reasonable range of target
967 // thickness, can go negative for extremely small target thicknesses.
968 // Here I protect against these unusual cases by taking the absolute value.
969
970 return fabs(Sigma2MS_Geant(thickness_m));
971}
972
973double CobremsGeneration::Sigma2MS_Kaune(double thickness_m)
974{
975 // Multiple scattering formula of Kaune et.al.
976 // with a correction factor from a multiple-scattering calculation
977 // taking into account the atomic and nuclear form factors for carbon.
978 //
979 // Note by RTJ, Oct. 13, 2008:
980 // I think this formula overestimates multiple scattering in thin targets
981 // like these diamond radiators, because it scales simply like sqrt(t).
982 // Although the leading behavior is sqrt(t/radlen), it should increase
983 // faster than that because of the 1/theta^2 tail of the Rutherford
984 // distribution that makes the central gaussian region swell with increasing
985 // number of scattering events. For comparison, I include below the PDG
986 // formula (sigma2MS_PDG), the Moliere formula used in the Geant3 simulation
987 // of gaussian multiple scattering (sigma2MS_Geant), and a Moliere fit for
988 // thin targets taken from reference Phys.Rev. vol.3 no.2, (1958), p.647
989 // (sigma2MS_Hanson). The latter two separate the gaussian part from the
990 // tails in different ways, but both agree that the central part is much
991 // more narrow than the formulation by Kaune et.al. below.
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
1002double CobremsGeneration::Sigma2MS_PDG(double thickness_m)
1003{
1004 // Evaluates the PDG formula for multiple scattering of the beam electron
1005 // inside the target crystal, with beta=1, charge=1. This formula is said
1006 // to be within 11% for t > 1e-3 rad.len.
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
1012double CobremsGeneration::Sigma2MS_Geant(double thickness_m)
1013{
1014 // Returns the Geant3 formula for the rms multiple-scattering angle
1015 // This formula is based on the theory of Moliere scattering. It contains
1016 // a cutoff parameter F that is used for the fractional integral of the
1017 // scattering probability distribution that is included in computing the
1018 // rms. This is needed because the complete distribution of scattering
1019 // angles connects smoothly from a central gaussian (small-angle
1020 // multiple-scattering regime) to a 1/theta^2 tail (large-angle Rutherford
1021 // scattering regime) through the so-called plural scattering region.
1022
1023 double rBohr = 0.52917721e-10; // m
1024 double F = 0.98; // probability cutoff in definition of sigma2MS
1025 double Z = fTargetCrystal.Z;
1026 double chi2cc = pow(0.39612e-2, 2) * Z * (Z + 1) *
1027 fTargetCrystal.density / 12; // GeV^2/m
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); // mean number of scatters
1032 double gnu = omega0 / (2 * (1 - F));
1033 return chi2c / (1 + pow(F, 2)) * ((1 + gnu) / gnu * log(1 + gnu) -1);
1034}
1035
1036double CobremsGeneration::Sigma2MS_Hanson(double thickness_m)
1037{
1038 // Formulation of the rms projected angle attributed to Hanson et.al.
1039 // in reference Phys.Rev. vol.3 no.2, (1958), p.647. This is just Moliere
1040 // theory used to give the 1/e angular width of the scattering distribution.
1041 // In the paper, though, they compare it with experiment for a variety of
1042 // metal foils down to 1e-4 rad.len. in thickness, and show excellent
1043 // agreement with the gaussian approximation out to 4 sigma or so. I
1044 // like this paper because of the excellent agreement between the theory
1045 // and experimental data.
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
1076void 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
1088double (CobremsGeneration::*Rate_dNtdx_1)(double) = &CobremsGeneration::Rate_dNtdx;
1089double (CobremsGeneration::*Rate_dNtdx_3)(double, double, double) = &CobremsGeneration::Rate_dNtdx;
1090double (CobremsGeneration::*Rate_dNcdx_1)(double) = &CobremsGeneration::Rate_dNcdx;
1091double (CobremsGeneration::*Rate_dNcdx_3)(double, double, double) = &CobremsGeneration::Rate_dNcdx;
1092double (CobremsGeneration::*Acceptance_1)(double) = &CobremsGeneration::Acceptance;
1093double (CobremsGeneration::*Acceptance_4)(double, double, double, double) = &CobremsGeneration::Acceptance;
1094double (CobremsGeneration::*Polarization_2)(double, double) = &CobremsGeneration::Polarization;
1095double (CobremsGeneration::*Polarization_3)(double, double, double) = &CobremsGeneration::Polarization;
1096
1097BOOST_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