Garfield++ v1r0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
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