Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
EnTransfCS.cpp
Go to the documentation of this file.
1#include <iomanip>
2#include <fstream>
3#include "wcpplib/clhep_units/WSystemOfUnits.h"
4#include "wcpplib/math/lorgamma.h"
5#include "wcpplib/math/tline.h"
6#include "heed++/code/EnTransfCS.h"
7#include "heed++/code/HeedMatterDef.h"
8/*
92003, I. Smirnov
10*/
11
12namespace Heed {
13
14EnTransfCS::EnTransfCS(double fparticle_mass, double fgamma_1,
15 int fs_primary_electron, HeedMatterDef* fhmd,
16 long fparticle_charge)
17 : particle_mass(fparticle_mass),
18 particle_charge(fparticle_charge),
19 beta(lorbeta(fgamma_1)),
20 beta2(0.0),
21 beta12(0.0),
22 gamma(fgamma_1 + 1.0),
23 gamma_1(fgamma_1),
24 s_simple_form(1),
25 s_primary_electron(fs_primary_electron),
26 hmd(fhmd),
27 quanC(0.0)
28#ifndef EXCLUDE_MEAN
29 ,
30 meanC(0.0),
31 meanC1(0.0),
32 meaneleC(0.0),
33 meaneleC1(0.0)
34#endif
35 {
36 mfunnamep("EnTransfCS::EnTransfCS(double fparticle_mass, double fbeta, "
37 "HeedMatterDef* fhmd)");
38 beta2 = beta * beta;
39 beta12 = 1.0 - beta2;
40 particle_tkin = particle_mass * gamma_1;
41 particle_ener = particle_mass * gamma;
42 if (s_primary_electron == 1) {
43 maximal_energy_trans = 0.5 * particle_tkin;
44 } else {
45 double rm2 = particle_mass * particle_mass;
46 double rme = ELMAS;
47 if (beta12 > 1.0e-10) {
48 maximal_energy_trans =
49 2.0 * rm2 * ELMAS * beta2 /
50 ((rm2 + rme * rme + 2.0 * rme * gamma * particle_mass) * (beta12));
51 if (maximal_energy_trans > particle_tkin) {
52 maximal_energy_trans = particle_tkin;
53 }
54 } else {
55 maximal_energy_trans = particle_tkin;
56 }
57 }
58 long qe = hmd->energy_mesh->get_q();
59 long ne;
60 log1C.put_qel(qe, 0.0);
61 log2C.put_qel(qe, 0.0);
62 chereC.put_qel(qe, 0.0);
63 chereCangle.put_qel(qe, 0.0);
64 Rreser.put_qel(qe, 0.0);
65#ifdef DEBUG_EnTransfCS
66 treser.put_qel(qe, 0.0);
67#endif
68 addaC.put_qel(qe, 0.0);
69#ifndef EXCLUDE_A_VALUES
70 addaC_a.put_qel(qe, 0.0);
71#endif
72
73 long qa = hmd->matter->qatom();
74 cher.put_qel(qa);
75 frezer.put_qel(qa);
76 adda.put_qel(qa);
77 fadda.put_qel(qa);
78 quan.put_qel(qa);
79#ifndef EXCLUDE_A_VALUES
80 cher_a.put_qel(qa);
81 adda_a.put_qel(qa);
82 fadda_a.put_qel(qa);
83 quan_a.put_qel(qa);
84#endif
85
86#ifndef EXCLUDE_VAL_FADDA
87 val_fadda.put_qel(qa);
88#ifndef EXCLUDE_A_VALUES
89 val_fadda_a.put_qel(qa);
90#endif
91#endif
92
93#ifndef EXCLUDE_MEAN
94 mean.put_qel(qa);
95#ifndef EXCLUDE_A_VALUES
96 mean_a.put_qel(qa);
97#endif
98#endif
99
100 for (long na = 0; na < qa; na++) {
101 long qs = hmd->apacs[na]->get_qshell();
102 cher[na].put_qel(qs);
103 frezer[na].put_qel(qs);
104 adda[na].put_qel(qs);
105 fadda[na].put_qel(qs);
106 quan[na].put_qel(qs);
107#ifndef EXCLUDE_A_VALUES
108 cher_a[na].put_qel(qs);
109 adda_a[na].put_qel(qs);
110 fadda_a[na].put_qel(qs);
111 quan_a[na].put_qel(qs);
112#endif
113#ifndef EXCLUDE_VAL_FADDA
114 val_fadda[na].put_qel(qs);
115#ifndef EXCLUDE_A_VALUES
116 val_fadda_a[na].put_qel(qs);
117#endif
118#endif
119
120#ifndef EXCLUDE_MEAN
121 mean[na].put_qel(qs);
122#ifndef EXCLUDE_A_VALUES
123 mean_a[na].put_qel(qs);
124#endif
125#endif
126 for (long ns = 0; ns < qs; ns++) {
127 cher[na][ns].put_qel(qe, 0.0);
128 frezer[na][ns].put_qel(qe, 0.0);
129 adda[na][ns].put_qel(qe, 0.0);
130 fadda[na][ns].put_qel(qe, 0.0);
131#ifndef EXCLUDE_A_VALUES
132 cher_a[na][ns].put_qel(qe, 0.0);
133 adda_a[na][ns].put_qel(qe, 0.0);
134 fadda_a[na][ns].put_qel(qe, 0.0);
135#endif
136 }
137 }
138 double r;
139 for (ne = 0; ne < qe; ne++) {
140 r = -hmd->epsi1[ne] + (1.0 + hmd->epsi1[ne]) * beta12;
141 r = r * r + beta2 * beta2 * hmd->epsi2[ne] * hmd->epsi2[ne];
142 r = 1.0 / sqrt(r);
143 r = log(r);
144 log1C[ne] = r;
145 }
146 for (ne = 0; ne < qe; ne++) {
147 r = 2.0 * 0.511 * beta2 / hmd->energy_mesh->get_ec(ne);
148 if (r > 0.0) {
149 r = log(r);
150 } else {
151 r = 0.0;
152 }
153 log2C[ne] = r;
154 }
155 double r0, rr12, rr22, r1, r2, r3;
156 double coefpa =
157 double(particle_charge * particle_charge) / (FSCON * beta2 * M_PI);
158 for (ne = 0; ne < qe; ne++) {
159 r0 = 1.0 + hmd->epsi1[ne];
160 r = -hmd->epsi1[ne] + r0 * beta12;
161 rr12 = r0 * r0;
162 rr22 = hmd->epsi2[ne] * hmd->epsi2[ne];
163 r1 = (-r0 * r + beta2 * rr22) / (rr12 + rr22);
164 r2 = hmd->epsi2[ne] * beta2 / r;
165 r3 = atan(r2);
166 if (r < 0) r3 = M_PI + r3;
167 chereCangle[ne] = r3;
168 chereC[ne] = (coefpa / hmd->eldens) * r1 * r3;
169 }
170
171 for (ne = 0; ne < qe; ne++) {
172 double ec = hmd->energy_mesh->get_ec(ne);
173 if (s_simple_form == 1) {
174 if (s_primary_electron == 0) {
175 Rreser[ne] =
176 1.0 / (ec * ec) * (1.0 - beta2 * ec / maximal_energy_trans);
177 } else {
178 Rreser[ne] = 1.0 / (ec * ec);
179 }
180 } else {
181 if (s_primary_electron == 0) {
182 Rreser[ne] =
183 1.0 / (ec * ec) * (1.0 - beta2 * ec / maximal_energy_trans +
184 ec * ec / (2.0 * particle_ener * particle_ener));
185 } else {
186 double delta = ec / particle_mass;
187 double pg2 = gamma * gamma;
188 Rreser[ne] =
189 beta2 / (particle_mass * particle_mass) * 1.0 / (pg2 - 1.0) *
190 (gamma_1 * gamma_1 * pg2 / (pow(delta * (gamma_1 - delta), 2.0)) -
191 (2.0 * pg2 + 2.0 * gamma - 1.0) / (delta * (gamma_1 - delta)) +
192 1.0);
193 }
194 }
195 }
196
197 double Z_mean = hmd->matter->Z_mean();
198 for (long na = 0; na < qa; na++) {
199 long qs = hmd->apacs[na]->get_qshell();
200 double at_weight_quan = hmd->matter->weight_quan(na);
201 for (long ns = 0; ns < qs; ns++) {
202 DynLinArr<double>& acher = cher[na][ns];
203#ifndef EXCLUDE_A_VALUES
204 DynLinArr<double>& acher_a = cher_a[na][ns];
205#endif
206 DynLinArr<double>& afrezer = frezer[na][ns];
207
208 for (ne = 0; ne < qe; ne++) {
209 double e1 = hmd->energy_mesh->get_e(ne);
210 double e2 = hmd->energy_mesh->get_e(ne + 1);
211 double ics = 0.;
212 if (s_use_mixture_thresholds == 1) {
213 ics = hmd->apacs[na]->get_integral_TICS(
214 ns, e1, e2, hmd->min_ioniz_pot) / (e2 - e1) * C1_MEV2_MBN;
215 } else {
216 ics = hmd->apacs[na]->get_integral_ICS(ns, e1, e2) / (e2 - e1) *
217 C1_MEV2_MBN;
218 }
219
220#ifndef EXCLUDE_A_VALUES
221 double acs = hmd->apacs[na]->get_integral_ACS(ns, e1, e2) / (e2 - e1) *
222 C1_MEV2_MBN;
223#endif
224 check_econd11a(ics, < 0,
225 "na=" << na << " ns=" << ns << " ne=" << ne << '\n',
226 mcerr);
227 if (hmd->ACS[ne] > 0.0) {
228 acher[ne] =
229 chereC[ne] * at_weight_quan * ics / (hmd->ACS[ne] * C1_MEV2_MBN);
230#ifndef EXCLUDE_A_VALUES
231 acher_a[ne] =
232 chereC[ne] * at_weight_quan * acs / (hmd->ACS[ne] * C1_MEV2_MBN);
233#endif
234 } else {
235 acher[ne] = 0.0;
236#ifndef EXCLUDE_A_VALUES
237 acher_a[ne] = 0.0;
238#endif
239 }
240 }
241 // Calculate the integral.
242 double s = 0.;
243 for (ne = 0; ne < qe; ne++) {
244 double e1 = hmd->energy_mesh->get_e(ne);
245 double ec = hmd->energy_mesh->get_ec(ne);
246 double e2 = hmd->energy_mesh->get_e(ne + 1);
247 r = hmd->apacs[na]->get_integral_ACS(ns, e1, e2) * C1_MEV2_MBN *
248 at_weight_quan;
249 // Here it must be ACS to satisfy sum rule for rezerford
250 check_econd11a(r, < 0.0, "na=" << na << " ns=" << ns << " na=" << na,
251 mcerr);
252 if (ec > hmd->min_ioniz_pot && ec < maximal_energy_trans) {
253 afrezer[ne] = (s + 0.5 * r) * coefpa * Rreser[ne] / Z_mean;
254 check_econd11a(afrezer[ne], < 0,
255 "na=" << na << " ns=" << ns << " na=" << na, mcerr);
256 } else {
257 afrezer[ne] = 0.0;
258 }
259 s += r;
260#ifdef DEBUG_EnTransfCS
261 treser[ne] += afrezer[ne];
262#endif
263 }
264 }
265 }
266 for (ne = 0; ne < qe; ++ne) {
267 double s = 0.0;
268#ifndef EXCLUDE_A_VALUES
269 double s_a = 0.0;
270#endif
271 double e1 = hmd->energy_mesh->get_e(ne);
272 double ec = hmd->energy_mesh->get_ec(ne);
273 double e2 = hmd->energy_mesh->get_e(ne + 1);
274 double sqepsi = pow((1 + hmd->epsi1[ne]), 2.0) + pow(hmd->epsi2[ne], 2.0);
275 for (long na = 0; na < qa; na++) {
276 double at_weight_quan = hmd->matter->weight_quan(na);
277 long qs = hmd->apacs[na]->get_qshell();
278 for (long ns = 0; ns < qs; ns++) {
279 double ics;
280 if (s_use_mixture_thresholds == 1) {
281 ics = hmd->apacs[na]->get_integral_TICS(
282 ns, e1, e2, hmd->min_ioniz_pot) / (e2 - e1) * C1_MEV2_MBN;
283 } else {
284 ics = hmd->apacs[na]->get_integral_ICS(ns, e1, e2) / (e2 - e1) *
285 C1_MEV2_MBN;
286 }
287#ifndef EXCLUDE_A_VALUES
288 double acs = hmd->apacs[na]->get_integral_ACS(ns, e1, e2) / (e2 - e1) *
289 C1_MEV2_MBN;
290#endif
291 double r1 =
292 at_weight_quan * log1C[ne] * coefpa * ics / (ec * Z_mean * sqepsi);
293 double r2 =
294 at_weight_quan * log2C[ne] * coefpa * ics / (ec * Z_mean * sqepsi);
295 double& r_adda = adda[na][ns][ne];
296 double& r_frezer = frezer[na][ns][ne];
297 r_adda = r1 + r2 + r_frezer;
298 if (r_adda < 0.0) r_adda = 0.0;
299
300#ifndef EXCLUDE_A_VALUES
301 double r1_a =
302 at_weight_quan * log1C[ne] * coefpa * acs / (ec * Z_mean * sqepsi);
303 double r2_a =
304 at_weight_quan * log2C[ne] * coefpa * acs / (ec * Z_mean * sqepsi);
305 double& r_adda_a = adda_a[na][ns][ne];
306 r_adda_a = r1_a + r2_a + frezer;
307 if (r_adda_a < 0.0) r_adda_a = 0.0;
308#endif
309 if (ec > hmd->min_ioniz_pot) {
310 r_adda += cher[na][ns][ne];
311 if (r_adda < 0.0) {
312 funnw.whdr(mcout);
313 mcout << "negative adda\n";
314 mcout << "na=" << na << " ns=" << ns << " ne=" << ne
315 << " adda[na][ns][ne] = " << adda[na][ns][ne] << '\n';
316 adda[na][ns][ne] = 0.0;
317 }
318 }
319#ifndef EXCLUDE_A_VALUES
320 adda_a[na][ns][ne] += cher[na][ns][ne];
321 check_econd11a(adda_a[na][ns][ne], < 0,
322 "na=" << na << " ns=" << ns << " na=" << na, mcerr);
323#endif
324 s += r_adda;
325#ifndef EXCLUDE_A_VALUES
326 s_a += r_adda_a;
327#endif
328 }
329 }
330 addaC[ne] = s;
331#ifndef EXCLUDE_A_VALUES
332 addaC_a[ne] = s_a;
333#endif
334 }
335
336 const double* aetemp = hmd->energy_mesh->get_ae();
337 PointCoorMesh<double, const double*> pcm_e(qe + 1, &(aetemp));
338 double emin = hmd->energy_mesh->get_emin();
339 double emax = hmd->energy_mesh->get_emax();
340
341 quanC = t_integ_step_ar<double, DynLinArr<double>,
342 PointCoorMesh<double, const double*> >(
343 pcm_e, addaC, emin, emax, 0) * hmd->xeldens;
344
345#ifndef EXCLUDE_A_VALUES
346 quanC_a = t_integ_step_ar<double, DynLinArr<double>,
347 PointCoorMesh<double, const double*> >(
348 pcm_e, addaC_a, emin, emax, 0) * hmd->xeldens;
349#endif
350
351#ifndef EXCLUDE_MEAN
352 meanC = t_integ_step_ar<double, DynLinArr<double>,
353 PointCoorMesh<double, const double*> >(
354 pcm_e, addaC, emin, emax, 1) * hmd->xeldens;
355
356#ifndef EXCLUDE_A_VALUES
357 meanC_a = t_integ_step_ar<double, DynLinArr<double>,
358 PointCoorMesh<double, const double*> >(
359 pcm_e, addaC_a, emin, emax, 1) * hmd->xeldens;
360#endif
361
362 meanCleft = meanC; // meanCleft does not have sense in this approach
363
364 if (s_simple_form == 1) {
365 if (s_primary_electron == 0) {
366 meanC1 = meanC;
367 if (maximal_energy_trans > hmd->energy_mesh->get_e(qe)) {
368 double e1 = hmd->energy_mesh->get_e(qe);
369 double e2 = maximal_energy_trans;
370 meanC1 += double(particle_charge * particle_charge) * 2.0 * M_PI /
371 (FSCON * FSCON * ELMAS * beta2) * hmd->xeldens *
372 (log(e2 / e1) - beta2 / maximal_energy_trans * (e2 - e1));
373 }
374 } else {
375 meanC1 = meanC;
376 if (maximal_energy_trans > hmd->energy_mesh->get_e(qe)) {
377 double e1 = hmd->energy_mesh->get_e(qe);
378 double e2 = maximal_energy_trans;
379 meanC1 +=
380 double(particle_charge * particle_charge) * 2.0 * M_PI /
381 (pow(FSCON, 2.0) * ELMAS * beta2) * hmd->xeldens * log(e2 / e1);
382 }
383 }
384 } else {
385 if (s_primary_electron == 0) {
386 meanC1 = meanC;
387 if (maximal_energy_trans > hmd->energy_mesh->get_e(qe)) {
388 double e1 = hmd->energy_mesh->get_e(qe);
389 double e2 = maximal_energy_trans;
390 meanC1 += double(particle_charge * particle_charge) * 2.0 * M_PI /
391 (FSCON * FSCON * ELMAS * beta2) * hmd->xeldens *
392 (log(e2 / e1) - beta2 / maximal_energy_trans * (e2 - e1) +
393 (e2 * e2 - e1 * e1) / (4.0 * particle_ener * particle_ener));
394 }
395#ifndef EXCLUDE_A_VALUES
396 meanC1_a = meanC_a;
397 if (maximal_energy_trans > hmd->energy_mesh->get_e(qe)) {
398 double e1 = hmd->energy_mesh->get_e(qe);
399 double e2 = maximal_energy_trans;
400 meanC1_a +=
401 double(particle_charge * particle_charge) * 2.0 * M_PI /
402 (pow(FSCON, 2.0) * ELMAS * beta2) * hmd->xeldens *
403 (log(e2 / e1) - beta2 / maximal_energy_trans * (e2 - e1) +
404 (e2 * e2 - e1 * e1) / (4.0 * particle_ener * particle_ener));
405 }
406#endif
407 }
408 }
409
410 meaneleC = meanC / hmd->W;
411 meaneleC1 = meanC1 / hmd->W;
412#endif
413
414 for (long na = 0; na < qa; na++) {
415 long qs = hmd->apacs[na]->get_qshell();
416 for (long ns = 0; ns < qs; ns++) {
417 quan[na][ns] = t_integ_step_ar<double, DynLinArr<double>,
418 PointCoorMesh<double, const double*> >(
419 pcm_e, adda[na][ns], emin, emax, 0) * hmd->xeldens;
420#ifndef EXCLUDE_A_VALUES
421 quan_a[na][ns] = t_integ_step_ar<double, DynLinArr<double>,
422 PointCoorMesh<double, const double*> >(
423 pcm_e, adda_a[na][ns], emin, emax, 0) * hmd->xeldens;
424#endif
425#ifndef EXCLUDE_MEAN
426 mean[na][ns] = t_integ_step_ar<double, DynLinArr<double>,
427 PointCoorMesh<double, const double*> >(
428 pcm_e, adda[na][ns], emin, emax, 1) * hmd->xeldens;
429#ifndef EXCLUDE_A_VALUES
430 mean_a[na][ns] = t_integ_step_ar<double, DynLinArr<double>,
431 PointCoorMesh<double, const double*> >(
432 pcm_e, adda_a[na][ns], emin, emax, 1) * hmd->xeldens;
433#endif
434#endif
435 }
436 }
437
438 for (long na = 0; na < qa; na++) {
439 long qs = hmd->apacs[na]->get_qshell();
440 for (long ns = 0; ns < qs; ns++) {
441 if (quan[na][ns] > 0.0)
442#ifndef EXCLUDE_VAL_FADDA
443 val_fadda[na][ns] =
444#endif
445 t_hispre_step_ar<double, DynLinArr<double>,
446 PointCoorMesh<double, const double*> >(
447 pcm_e, adda[na][ns], fadda[na][ns]);
448
449#ifndef EXCLUDE_A_VALUES
450 if (quan_a[na][ns] > 0.0)
451#ifndef EXCLUDE_VAL_FADDA
452 val_fadda_a[na][ns] =
453#endif
454 t_hispre_step_ar<double, DynLinArr<double>,
455 PointCoorMesh<double, const double*> >(
456 pcm_e, adda_a[na][ns], fadda_a[na][ns]);
457#endif
458 }
459 }
460
461 length_y0 = DynLinArr<double>(qe, 0.0);
462 for (ne = 0; ne < qe; ne++) {
463 double k0 = hmd->energy_mesh->get_ec(ne) / PLANKCLIGHT;
464 double det_value = 1.0 / (gamma * gamma) - hmd->epsi1[ne] * beta2;
465 if (det_value <= 0.0) {
466 length_y0[ne] = 0.0;
467 } else {
468 length_y0[ne] = beta / k0 * 1.0 / sqrt(det_value);
469 }
470 }
471
472#ifdef CLEAR_ARRAYS
473 log1C.clear();
474 log2C.clear();
475 chereC.clear();
476 chereCangle.clear();
477 Rreser.clear();
478#ifdef DEBUG_EnTransfCS
479 treser.clear();
480#endif
481 std::ofstream dcsfile;
482 dcsfile.open("dcs.txt", std::ios::out);
483 dcsfile << "# energy [MeV] vs. diff. cs per electron [Mbarn / MeV]\n";
484 for (int i = 0; i < qe; ++i) {
485 dcsfile << hmd->energy_mesh->get_ec(i) << " " << addaC[i] / C1_MEV2_MBN
486 << "\n";
487 }
488 dcsfile.close();
489
490 addaC.clear();
491#ifndef EXCLUDE_A_VALUES
492 addaC_a.clear();
493#endif
494 cher.clear();
495#ifndef EXCLUDE_A_VALUES
496 cher_a.clear();
497#endif
498 frezer.clear();
499 adda.clear();
500#ifndef EXCLUDE_A_VALUES
501 adda_a.clear();
502#endif
503#ifndef EXCLUDE_A_VALUES
504 fadda_a.clear();
505#endif
506#ifndef EXCLUDE_VAL_FADDA
507 val_fadda.clear();
508#ifndef EXCLUDE_A_VALUES
509 val_fadda_a.clear();
510#endif
511#endif
512#ifndef EXCLUDE_MEAN
513 mean.clear();
514#ifndef EXCLUDE_A_VALUES
515 mean_a.clear();
516#endif
517#endif
518
519#endif
520
521}
522
523void EnTransfCS::print(std::ostream& file, int l) const {
524 if (l <= 0) return;
525 Ifile << "EnTransfCS(l=" << l << "):\n";
526 indn.n += 2;
527 Ifile << "particle_mass=" << particle_mass
528 << " particle_tkin=" << particle_tkin
529 << " particle_ener=" << particle_ener
530 << " particle_charge=" << particle_charge << std::endl;
531 Ifile << "beta=" << beta << " beta2=" << beta2 << " beta12=" << beta12
532 << " gamma=" << gamma << std::endl;
533 Ifile << "maximal_energy_trans=" << maximal_energy_trans << std::endl;
534 Ifile << "s_primary_electron=" << s_primary_electron << std::endl;
535 Ifile << "hmd:\n";
536 hmd->print(file, 1);
537#ifndef EXCLUDE_MEAN
538#ifndef EXCLUDE_A_VALUES
539 Ifile << "quanC=" << quanC << " quanC_a=" << quanC_a << '\n';
540 Ifile << "meanC=" << meanC << " meanC_a=" << meanC_a << '\n';
541 Ifile << " meanCleft=" << meanCleft << " meaneleC=" << meaneleC << '\n';
542 Ifile << "meanC1=" << meanC1 << " meanC1_a=" << meanC1_a << '\n';
543#else
544 Ifile << "quanC=" << quanC << '\n';
545 Ifile << "meanC=" << meanC << '\n';
546 Ifile << " meanCleft=" << meanCleft << " meaneleC=" << meaneleC << '\n';
547 Ifile << "meanC1=" << meanC1 << '\n';
548#endif
549 Ifile << " meaneleC1=" << meaneleC1 << '\n';
550#else
551#ifndef EXCLUDE_A_VALUES
552 Ifile << "quanC=" << quanC << " quanC_a=" << quanC_a << '\n';
553#else
554 Ifile << "quanC=" << quanC << '\n';
555#endif
556#endif
557 if (l > 2) {
558 long qe = hmd->energy_mesh->get_q();
559 long ne;
560 if (l > 4) {
561#ifdef DEBUG_EnTransfCS
562 Ifile << " enerc, log1C, log2C, chereC, addaC, "
563 "chereCangle Rreser treser length_y0\n";
564#else
565 Ifile << " enerc, log1C, log2C, chereC, addaC, "
566 "chereCangle Rreser length_y0\n";
567#endif
568 for (ne = 0; ne < qe; ne++) {
569 Ifile << std::setw(12) << hmd->energy_mesh->get_ec(ne) << std::setw(12)
570 << log1C[ne] << std::setw(12) << log2C[ne] << std::setw(12)
571 << chereC[ne] << std::setw(12) << addaC[ne] << std::setw(12)
572 << chereCangle[ne] << std::setw(12) << Rreser[ne]
573#ifdef DEBUG_EnTransfCS
574 << std::setw(12) << treser[ne]
575#endif
576 << std::setw(12) << length_y0[ne] << '\n';
577 }
578 }
579 if (l > 3) {
580 long qa = hmd->matter->qatom();
581 long na;
582 Iprintn(file, hmd->matter->qatom());
583 for (na = 0; na < qa; na++) {
584 Iprintn(file, na);
585 long qs = hmd->apacs[na]->get_qshell();
586 long ns;
587 Iprintn(file, hmd->apacs[na]->get_qshell());
588 for (ns = 0; ns < qs; ns++) {
589 Iprintn(file, ns);
590#ifndef EXCLUDE_MEAN
591 Ifile << "quan =" << std::setw(13) << quan[na][ns]
592 << " mean =" << std::setw(13) << mean[na][ns] << '\n';
593#ifndef EXCLUDE_A_VALUES
594 Ifile << "quan_a =" << std::setw(13) << quan_a[na][ns]
595 << " mean_a=" << std::setw(13) << mean_a[na][ns] << '\n';
596#endif
597#else
598 Ifile << "quan =" << std::setw(13) << quan[na][ns] << '\n';
599#ifndef EXCLUDE_A_VALUES
600 Ifile << "quan_a =" << std::setw(13) << quan_a[na][ns] << '\n';
601#endif
602#endif
603#ifndef EXCLUDE_VAL_FADDA
604 Ifile << "val_fadda=" << std::setw(13) << val_fadda[na][ns]
605#ifndef EXCLUDE_A_VALUES
606 << " val_fadda_a=" << std::setw(13) << val_fadda_a[na][ns]
607#endif
608 << '\n';
609#endif
610 if (l > 5) {
611 Ifile << " enerc, cher, cher_a, frezer, adda, "
612 " adda_a, fadda, fadda_a\n";
613 for (ne = 0; ne < qe; ne++) {
614 Ifile << std::setw(12) << hmd->energy_mesh->get_ec(ne)
615 // << std::setw(12) << flog1[na][ns][ne]
616 // << std::setw(12) << flog2[na][ns][ne]
617 << std::setw(12) << cher[na][ns][ne]
618#ifndef EXCLUDE_A_VALUES
619 << std::setw(12) << cher_a[na][ns][ne]
620#endif
621 // << std::setw(12) << rezer[na][ns][ne]
622 << std::setw(12) << frezer[na][ns][ne] << std::setw(12)
623 << adda[na][ns][ne]
624#ifndef EXCLUDE_A_VALUES
625 << std::setw(12) << adda_a[na][ns][ne]
626#endif
627 << std::setw(12) << fadda[na][ns][ne]
628#ifndef EXCLUDE_A_VALUES
629 << std::setw(12) << fadda_a[na][ns][ne]
630#endif
631 << '\n';
632 }
633 }
634 }
635 }
636 }
637 }
638 indn.n -= 2;
639
640}
641
642std::ostream& operator<<(std::ostream& file, const EnTransfCSType& f) {
643 mfunname("std::ostream& operator << (std::ostream& file, const "
644 "EnTransfCSType& f)");
645 if (f.etcs.get() == NULL) {
646 Ifile << "EnTransfCSType: type is not initialized\n";
647 } else {
648 Ifile << "EnTransfCSType: =";
649 f.etcs->print(file, 1);
650 }
651 return file;
652}
653
654}
pow
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:336
sqrt
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313
EXCLUDE_A_VALUES
#define EXCLUDE_A_VALUES
Definition: EnTransfCS.h:20
check_econd11a
#define check_econd11a(a, signb, add, stream)
Definition: FunNameStack.h:395
mfunnamep
#define mfunnamep(string)
Definition: FunNameStack.h:77
mfunname
#define mfunname(string)
Definition: FunNameStack.h:67
DynLinArr::clear
void clear(void)
Definition: AbsArr.h:450
DynLinArr::put_qel
void put_qel(long fqel)
Definition: AbsArr.h:774
Heed::EnTransfCSType::etcs
PassivePtr< EnTransfCS > etcs
Definition: EnTransfCS.h:151
Heed::EnTransfCS::particle_charge
long particle_charge
Charge in units of electron charge (used square, sign does not matter).
Definition: EnTransfCS.h:43
Heed::EnTransfCS::quan
DynLinArr< DynLinArr< double > > quan
Definition: EnTransfCS.h:134
Heed::EnTransfCS::Rreser
DynLinArr< double > Rreser
Definition: EnTransfCS.h:70
Heed::EnTransfCS::particle_ener
double particle_ener
Total energy [MeV].
Definition: EnTransfCS.h:41
Heed::EnTransfCS::EnTransfCS
EnTransfCS(void)
Constructors.
Definition: EnTransfCS.h:29
Heed::EnTransfCS::mean
DynLinArr< DynLinArr< double > > mean
Definition: EnTransfCS.h:139
Heed::EnTransfCS::cher
DynLinArr< DynLinArr< DynLinArr< double > > > cher
Definition: EnTransfCS.h:111
Heed::EnTransfCS::particle_mass
double particle_mass
Particle mass [MeV].
Definition: EnTransfCS.h:37
Heed::EnTransfCS::chereCangle
DynLinArr< double > chereCangle
Definition: EnTransfCS.h:69
Heed::EnTransfCS::s_primary_electron
int s_primary_electron
Definition: EnTransfCS.h:60
Heed::EnTransfCS::maximal_energy_trans
double maximal_energy_trans
Max. energy transfer [MeV].
Definition: EnTransfCS.h:52
Heed::EnTransfCS::adda
DynLinArr< DynLinArr< DynLinArr< double > > > adda
Sum.
Definition: EnTransfCS.h:115
Heed::EnTransfCS::hmd
PassivePtr< HeedMatterDef > hmd
Definition: EnTransfCS.h:62
Heed::EnTransfCS::log2C
DynLinArr< double > log2C
Definition: EnTransfCS.h:67
Heed::EnTransfCS::length_y0
DynLinArr< double > length_y0
Definition: EnTransfCS.h:145
Heed::EnTransfCS::particle_tkin
double particle_tkin
Kinetic energy [MeV].
Definition: EnTransfCS.h:39
Heed::EnTransfCS::quanC
double quanC
Integrated (ionization) cross-section.
Definition: EnTransfCS.h:79
Heed::EnTransfCS::fadda
DynLinArr< DynLinArr< DynLinArr< double > > > fadda
Integral, normalised to unity.
Definition: EnTransfCS.h:117
Heed::EnTransfCS::addaC
DynLinArr< double > addaC
Sum of (ionization) differential cross-section terms.
Definition: EnTransfCS.h:77
Heed::EnTransfCS::print
virtual void print(std::ostream &file, int l) const
Definition: EnTransfCS.cpp:523
Heed::EnTransfCS::meanCleft
double meanCleft
Definition: EnTransfCS.h:95
Heed::EnTransfCS::chereC
DynLinArr< double > chereC
Definition: EnTransfCS.h:68
Heed::EnTransfCS::frezer
DynLinArr< DynLinArr< DynLinArr< double > > > frezer
Rutherford term.
Definition: EnTransfCS.h:113
Heed::EnTransfCS::log1C
DynLinArr< double > log1C
Definition: EnTransfCS.h:66
Heed
Definition: BGMesh.cpp:3
Heed::operator<<
std::ostream & operator<<(std::ostream &file, const BGMesh &bgm)
Definition: BGMesh.cpp:22
Heed::lorbeta
double lorbeta(const double gamma_1)
Definition: lorgamma.cpp:22
Heed::ELMAS
const double ELMAS
Definition: heed++/code/PhysicalConstants.h:13
Heed::s_use_mixture_thresholds
const int s_use_mixture_thresholds
Definition: HeedMatterDef.h:26
Heed::C1_MEV2_MBN
const double C1_MEV2_MBN
Definition: heed++/code/PhysicalConstants.h:22
Heed::FSCON
const double FSCON
Definition: heed++/code/PhysicalConstants.h:14
Heed::PLANKCLIGHT
const double PLANKCLIGHT
Definition: heed++/code/PhysicalConstants.h:30
indn
indentation indn
Definition: prstream.cpp:13
mcout
#define mcout
Definition: prstream.h:133
Ifile
#define Ifile
Definition: prstream.h:207
mcerr
#define mcerr
Definition: prstream.h:135
Iprintn
#define Iprintn(file, name)
Definition: prstream.h:216
t_integ_step_ar
T t_integ_step_ar(const M &mesh, const D &y, T x1, T x2, int xpower)
Definition: tline.h:1810