File: | scratch/gluex/scan-build-work/hdgeant4/src/CobremsGeneration.cc |
Location: | line 638, column 16 |
Description: | Value stored to 'kmin' is never read |
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_VERBOSITY0 0 |
52 | #define BOOST_PYTHON_WRAPPING1 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 | |
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 | // 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(false); |
84 | |
85 | #if COBREMS_GENERATOR_VERBOSITY0 > 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); // 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 | |
106 | void 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 | |
156 | double 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 | |
186 | void 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 | |
218 | void 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 | |
274 | void 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 | |
382 | double 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 | |
398 | double 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 | |
411 | void 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 | |
426 | CobremsGeneration::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 | |
448 | CobremsGeneration &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 | |
471 | CobremsGeneration::~CobremsGeneration() { } |
472 | |
473 | double 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 | |
483 | double 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 | |
491 | double 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 | |
512 | double 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 | |
520 | double 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 | |
535 | double 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 | |
562 | double 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_VERBOSITY0 > 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_VERBOSITY0 > 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 | // 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 | |
691 | double 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 | |
739 | double 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 | |
767 | double 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 | |
780 | double 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 | |
792 | double 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 | |
805 | double 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 | |
817 | double 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 | |
841 | double 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 | |
860 | double 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 | |
921 | void 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 | |
961 | double 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 | |
973 | double 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 | |
1002 | double 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 | |
1012 | double 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 | |
1036 | double 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_WRAPPING1 |
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)void init_module_libcobrems(); extern "C" __attribute__ ((__visibility__ ("default"))) void initlibcobrems() { boost::python::detail:: init_module( "libcobrems",&init_module_libcobrems); } void init_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 |