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
PhotoAbsCS.cpp
Go to the documentation of this file.
1#include <cmath>
2#include <fstream>
3#include <iomanip>
10#include "heed++/code/EnergyMesh.h" // From this used only make_log_mesh_ec
11 /*
122004, I. Smirnov
13*/
14
15//#define DEBUG_PRINT_get_escape_particles
16//#define DEBUG_ignore_non_standard_channels
17
18//#define ALWAYS_LINEAR_INTERPOLATION // how the paper was computed
19
20namespace Heed {
21
22int sign_nonlinear_interpolation(double e1, double cs1, double e2, double cs2,
23 double threshold) {
24#ifdef ALWAYS_LINEAR_INTERPOLATION
25 return 0;
26#else
27 // normal case:
28 if (cs2 >= cs1 || cs2 <= 0 || e1 < 300.0e-6 || e1 < 1.5 * threshold) {
29 //if(cs2 >= cs1 || cs2 <= 0) // for debug
30 return 0;
31 } else {
32 const double pw = log(cs1 / cs2) / log(e1 / e2);
33 //Iprintn(mcout, pw);
34 if (pw >= -1.0) {
35 // good case for linear interpolation
36 return 0;
37 } else if (pw < -5.0) {
38 // unclear odd case, power would be dangerous
39 return 0;
40 } else {
41 // non-linear interpolation
42 return 1;
43 }
44 }
45 return 0;
46#endif
47}
48
49double my_integr_fun(double xp1, double yp1, double xp2, double yp2,
50 double xmin, double /*xmax*/, double x1, double x2) {
51 double res = 0.;
52 if (sign_nonlinear_interpolation(xp1, yp1, xp2, yp2, xmin) == 1) {
53 res = t_integ_power_2point<double>(xp1, yp1, xp2, yp2, x1, x2);
54 } else {
55 res = t_integ_straight_2point<double>(xp1, yp1, xp2, yp2, x1, x2, 0, 1);
56 }
57 return res;
58}
59
60double my_val_fun(double xp1, double yp1, double xp2, double yp2, double xmin,
61 double /*xmax*/, double x) {
62 double res = 0.;
63 //Iprint4n(mcout, xp1, yp1, xp2, yp2);
64 if (sign_nonlinear_interpolation(xp1, yp1, xp2, yp2, xmin) == 1) {
65 // Non-linear interpolation
66 res = t_value_power_2point<double>(xp1, yp1, xp2, yp2, x);
67 } else {
68 // Linear interpolation
69 res = t_value_straight_2point<double>(xp1, yp1, xp2, yp2, x, 1);
70 }
71 return res;
72}
73
75 double x1, double x2, double threshold) {
76 // Fit table by a straight line and integrate the area below it.
77 mfunname("double glin_integ_ar(DynLinArr< double > x ...)");
78
80 double s = t_integ_generic_point_ar<
82 pcmd, y, &my_integr_fun, x1, x2, 1, threshold, 0, DBL_MAX);
83 /*
84 // compare with old version
85 double old_res = old_glin_integ_ar(x, y, q, x1, x2, threshold);
86 if (!apeq(s, old_res)) {
87 Iprint2n(mcout, s, old_res);
88 Iprint4n(mcout, q, x1, x2, threshold);
89 mcout<<"x,y:\n";
90 for (long n = 0; n < q; ++n) {
91 Iprint3n(mcout, n, x[n], y[n]);
92 }
93 spexit(mcerr);
94 }
95 */
96 return s;
97}
98
99/*
100double glin_val_ar(DynLinArr<double> x, DynLinArr<double> y, long q, double ex,
101double threshold) {
102 // fit table by a straight line and integrate the area below it.
103 mfunname("double glin_val_ar(DynLinArr< double > x ...)");
104
105 PointCoorMesh< double, DynLinArr< double > > pcmd( q, &x);
106 double s = t_value_generic_point_ar<double, DynLinArr<double>,
107 PointCoorMesh<double, DynLinArr<double> > >
108 (pcmd, y, &my_val_fun, ex
109 1, threshold, 0, DBL_MAX);
110 return s;
111}
112*/
113/*
114double old_glin_integ_ar(DynLinArr< double > x, DynLinArr< double > y,
115 long q, double x1, double x2,
116 double threshold )
117// fit table by a straight line and integrate the area below it.
118{
119 mfunname("double old_glin_integ_ar(DynLinArr< double > x ...)");
120 //mcout<<"old_glin: q="<<q<<" x1="<<x1<<" x2="<<x2
121 // <<" threshold="<<threshold<<'\n';
122 long i;
123 double a,b;
124 double s=0;
125 if(q<=0)return 0;
126 check_econd11(threshold , <= 0.0, mcerr); // to avoid the use for
127 // which it is not designed.
128 // The table is extrapolated to threshold.
129 // Therefore it should have valid nonzero value.
130 if( x2 < x1 || x2 < threshold || x1 > x[q-1] ) return 0;
131 if(x1 < threshold) x1 = threshold;
132 if(x2 > x[q-1]) x2 = x[q-1];
133 // to find n1 and n2: indexes of next points after x1 and x2
134 long n1,n2;
135 for(i=0; i<q; i++ )
136 if(x[i] > x1 )
137 { n1=i; break; }
138 n2=q-1;
139 for(i=n1; i<q; i++ )
140 if( x[i] >= x2 )
141 { n2=i; break; }
142 //Iprintn(mcout, n2);
143 long nr; // index of current interval
144 if( x1 < x[0] )
145 //{ xt1=x[0]; nr=0; }
146 {
147 nr=0;
148 if(n2 == 0) n2 = 1; // new correction
149 }
150 else
151 { nr=n1-1; }
152 double xr1,xr2; // limits of integration at the current step
153 xr2=x1;
154 //Iprintn(mcout, n2);
155 //Iprintn(mcout, q);
156 //Iprint3n(mcout, nr, n2, xr2);
157 for( ; nr<n2; nr++)
158 {
159
160 //for debug:
161 check_econd12(nr , > , q-2 , mcerr);
162
163 //if(nr>q-2)
164 //{
165 //cerr<<"step_integ_ar:wrong nr=",nr,"\n";
166 //exit(1);
167 //}
168
169 xr1=xr2;
170 if( nr < q-1 )
171 xr2 = (x2 < x[nr+1]) ? x2 : x[nr+1];
172 else
173 xr2 = x2;
174 if(sign_nonlinear_interpolation(x[nr], y[nr],
175 x[nr+1], y[nr+1], threshold) == 1)
176 {
177 double pw = log(y[nr]/y[nr+1]) / log(x[nr]/x[nr+1]);
178 check_econd11(pw , >= -1.0 , mcout);
179 double k = y[nr] * pow( x[nr] , -pw );
180 double t = k/(1+pw)*( pow( xr2 , (pw+1) ) - pow( xr1 , (pw+1) ) );
181 check_econd11a(t , < 0.0, "non-linear interpolation\n"
182 <<" x[nr]="<<x[nr]
183 <<" y[nr]="<<y[nr]<<'\n'
184 <<" x[nr+1]="<<x[nr+1]
185 <<" y[nr+1]="<<y[nr+1]<<'\n'
186 <<" threshold="<<threshold<<'\n'
187 <<" pw="<<pw
188 <<" k="<<k
189 <<'\n', mcout);
190
191 //if(nr == 0)
192 //{
193 // mcout<<"non-linear interpolation\n"
194 // <<" x[nr]="<<x[nr]
195 // <<" y[nr]="<<y[nr]<<'\n'
196 // <<" x[nr+1]="<<x[nr+1]
197 // <<" y[nr+1]="<<y[nr+1]<<'\n'
198 // <<" threshold="<<threshold<<'\n'
199 // <<" pw="<<pw
200 // <<" k="<<k
201 // <<" t="<<t<<'\n';
202 //}
203 s += t;
204 }
205 else
206 {
207 a = (y[nr+1] - y[nr])/(x[nr+1] - x[nr]);
208 b = y[nr];
209 if(nr == 0)
210 {
211 double c1 = a * (xr1 - x[nr]) + b;
212 double c2 = a * (xr2 - x[nr]) + b;
213 if(c2 < 0.0) // new correction
214 {
215 return 0.0;
216 }
217 if(c1 < 0.0)
218 {
219 //mcout<<"linear interpolation\n"
220 // <<"c="<<c<<'\n'
221 // <<"x1="<<x1<<" x2="<<x2<<" q="<<q
222 // <<" nr="<<nr<<'\n'
223 // <<" x[nr]="<<x[nr]<<" y[nr]="<<y[nr]<<'\n'
224 // <<" x[nr+1]="<<x[nr+1]<<" y[nr+1]="<<y[nr+1]<<'\n'
225 // <<" xr1,2="<<xr1<<' '<<xr2<<'\n'
226 // <<" a,b = "<<a<<' '<<b<<'\n';
227 // find zero and
228 // change xr1
229 check_econd11(a , == 0.0 , mcerr);
230 xr1 = (a * x[nr] - b) / a;
231 // for debug:
232 double t = 0.5*a*(xr2*xr2 - xr1*xr1) + (b - a*x[nr])*(xr2 - xr1);
233 //mcout<<"linear interpolation\n"
234 // <<"xr1="<<xr1<<" a * (xr1- x[nr]) + b="<<a * (xr1- x[nr]) + b
235 // <<" t="<<t<<'\n';
236 //exit(0);
237 }
238 }
239 double t = 0.5*a*(xr2*xr2 - xr1*xr1) + (b - a*x[nr])*(xr2 - xr1);
240 check_econd11a(t , < 0.0,
241 "linear interpolation\n"
242 <<"x1="<<x1<<" x2="<<x2<<" q="<<q
243 <<" nr="<<nr<<'\n'
244 <<" x[nr]="<<x[nr]<<" y[nr]="<<y[nr]<<'\n'
245 <<" x[nr+1]="<<x[nr+1]<<" y[nr+1]="<<y[nr+1]<<'\n'
246 <<" xr1,2="<<xr1<<' '<<xr2<<'\n'
247 <<" a,b = "<<a<<' '<<b<<'\n' , mcout);
248 //if(nr == 0)
249 //{
250 // mcout<<"linear interpolation\n"
251 // <<" xr1="<<xr1
252 // <<" xr2="<<xr2<<'\n'
253 // <<" x[nr]="<<x[nr]
254 // <<" y[nr]="<<y[nr]<<'\n'
255 // <<" x[nr+1]="<<x[nr+1]
256 // <<" y[nr+1]="<<y[nr+1]<<'\n'
257 // <<" threshold="<<threshold<<'\n'
258 // <<" a="<<a
259 // <<" b="<<b
260 // <<" t="<<t<<'\n';
261 //}
262 s+= t;
263 }
264 }
265 return s;
266}
267*/
268
269PhotoAbsCS::PhotoAbsCS(void) : Z(0), threshold(0.0) { ; }
270
271PhotoAbsCS::PhotoAbsCS(const String& fname, int fZ, double fthreshold)
272 : name(fname), Z(fZ), threshold(fthreshold) {
273 ;
274}
275
276void PhotoAbsCS::print(std::ostream& file, int l) const {
277 if (l > 0) {
278 Ifile << "PhotoAbsCS: name=" << name << " Z = " << Z
279 << " threshold = " << threshold << std::endl;
280 }
281}
282
283OveragePhotoAbsCS::OveragePhotoAbsCS(PhotoAbsCS* apacs, double fwidth, // MeV
284 double fstep, // MeV
285 long fmax_q_step)
286 : real_pacs(apacs, do_clone),
287 width(fwidth),
288 max_q_step(fmax_q_step),
289 step(fstep) {
290 mfunname("OveragePhotoAbsCS::OveragePhotoAbsCS(...)");
291 check_econd11(apacs, == NULL, mcerr);
292 if (fwidth > 0.0) {
293 check_econd11(fstep, >= 0.6 * fwidth, mcerr); // 0.5 is bad but OK
294 }
295 /* I do not understand why the access is not allowed below.
296 So I call functions
297 name =
298 apacs->name;
299 //real_pacs->name;
300 Z = real_pacs->Z;
301 threshold = real_pacs->threshold;
302 */
303 name = real_pacs->get_name();
304 Z = real_pacs->get_Z();
305 threshold = real_pacs->get_threshold();
306}
307
308double OveragePhotoAbsCS::get_CS(double energy) const {
309 mfunname("double OveragePhotoAbsCS::get_CS(double energy) const");
310 //mcout<<"OveragePhotoAbsCS::get_CS is started\n";
311 //mcout<<"OveragePhotoAbsCS::get_CS:\n";
312 if (width == 0.0) {
313 // for no modification:
314 return real_pacs->get_CS(energy);
315 } else {
316 double w2 = width * 0.5;
317 double e1 = energy - w2;
318 if (e1 < 0.0) e1 = 0.0;
319 double res = real_pacs->get_integral_CS(e1, energy + w2) / width;
320 //Iprint4n(mcout, e1, energy + w2, width, res);
321 //return real_pacs->get_integral_CS(e1, energy + w2)/width;
322 return res;
323 }
324}
325
327 double energy2) const {
328 mfunname("double OveragePhotoAbsCS::get_integral_CS(double energy1, double "
329 "energy2) const");
330 //mcout<<"OveragePhotoAbsCS::get_integral_CS is started\n";
331 if (width == 0.0 || energy1 >= energy2) {
332 // for no modification:
333 return real_pacs->get_integral_CS(energy1, energy2);
334 } else {
335 long q = long((energy2 - energy1) / step);
336 if (q > max_q_step) {
337 return real_pacs->get_integral_CS(energy1, energy2);
338 } else {
339 //if(q == 0)q = 1;
340 q++;
341 double rstep = (energy2 - energy1) / q;
342 double x0 = energy1 + 0.5 * rstep;
343 double s = 0.0;
344 for (long n = 0; n < q; n++) {
345 s += get_CS(x0 + rstep * n);
346 }
347 s *= rstep;
348 return s;
349 }
350 }
351
352}
353
354void OveragePhotoAbsCS::scale(double fact) {
355 mfunname("void OveragePhotoAbsCS::scale(double fact)");
356 real_pacs->scale(fact);
357}
358
359void OveragePhotoAbsCS::print(std::ostream& file, int l) const {
360 mfunname("void PhotoAbsCS::print(std::ostream& file, int l) const");
361 Ifile << "OveragePhotoAbsCS: width = " << width << " step=" << step
362 << " max_q_step=" << max_q_step << '\n';
363 indn.n += 2;
364 real_pacs->print(file, l);
365 indn.n -= 2;
366}
367
369 : PhotoAbsCS("H", 1, 15.43e-6), prefactor(1.) {}
370
371double HydrogenPhotoAbsCS::get_CS(double energy) const {
372 if (energy < threshold) {
373 return 0.0;
374 } else {
375 if (energy != DBL_MAX) {
376 return 0.5 * // accounts one atom instead of two
377 prefactor * 0.0535 * (pow(100.0e-6 / energy, 3.228));
378 } else {
379 return 0.0;
380 }
381 }
382}
383
385 double energy2) const {
386 if (energy2 < threshold) {
387 return 0.0;
388 } else {
389 if (energy1 < threshold) {
390 energy1 = threshold; // local var.
391 }
392 if (energy2 == DBL_MAX) {
393 return 0.5 * // accounts one atom instead of two
394 prefactor * 0.0535 * pow(100.0e-6, 3.228) / 2.228 *
395 (1.0 / pow(energy1, 2.228));
396 } else {
397 return 0.5 * // accounts one atom instead of two
398 prefactor * 0.0535 * pow(100.0e-6, 3.228) / 2.228 *
399 (1.0 / pow(energy1, 2.228) - 1.0 / pow(energy2, 2.228));
400 }
401 }
402}
403
404void HydrogenPhotoAbsCS::scale(double fact) { prefactor = fact; }
405
406void HydrogenPhotoAbsCS::print(std::ostream& file, int l) const {
407 if (l > 0) {
408 Ifile << "HydrogenPhotoAbsCS: name=" << name << " Z = " << Z
409 << " threshold = " << threshold << std::endl;
410 }
411}
412
414 double fthreshold,
415 const String& ffile_name)
416 : PhotoAbsCS(fname, fZ, fthreshold), file_name(ffile_name) {
417 mfunnamep("SimpleTablePhotoAbsCS::SimpleTablePhotoAbsCS(...)");
418#ifdef USE_STLSTRING
419 std::ifstream file(file_name.c_str());
420#else
421 std::ifstream file(file_name);
422#endif
423 if (!file) {
424 funnw.ehdr(mcerr);
425 mcerr << "cannot open file " << file_name << std::endl;
426 spexit(mcerr);
427 }
428 long q = 0; // number of read values
429 ener = DynLinArr<double>(10, 0.0);
430 cs = DynLinArr<double>(10, 0.0);
431 do {
432 file >> ener[q];
433 //Iprintn(mcout, file.good());
434 //Iprintn(mcout, file.eof());
435 //Iprintn(mcout, file.fail());
436 //Iprintn(mcout, file.bad());
437 //file.clear(); // idiotic thing necessary for the new compiler
438 if (!file.good()) break;
439 check_econd11(ener[q], < 0.0, mcerr);
440 if (q > 0) {
441 check_econd12(ener[q], <, ener[q - 1], mcerr);
442 }
443 ener[q] *= 1.0e-6; // convertion to MeV
444 file >> cs[q];
445 if (!file.good()) break;
446 check_econd11(cs[q], < 0.0, mcerr);
447 q++;
448 if (q == ener.get_qel()) { // increasing the size of arrays
449 long q1 = q * 2;
450 ener.put_qel(q1);
451 cs.put_qel(q1);
452 }
453 } while (1);
454 if (q != ener.get_qel()) { // reducing the size if necessary
455 ener.put_qel(q);
456 cs.put_qel(q);
457 }
458
459}
461 double fthreshold,
462 const DynLinArr<double>& fener,
463 const DynLinArr<double>& fcs)
464 : PhotoAbsCS(fname, fZ, fthreshold),
465 file_name("none"),
466 ener(fener),
467 cs(fcs) {
468 mfunname("SimpleTablePhotoAbsCS::SimpleTablePhotoAbsCS(...)");
469 check_econd12(ener.get_qel(), !=, cs.get_qel(), mcerr);
470}
471
473 double fthreshold, int l,
474 double E0, double yw, double ya,
475 double P, double sigma)
476 : PhotoAbsCS(fname, fZ, fthreshold) {
477 mfunname("SimpleTablePhotoAbsCS::SimpleTablePhotoAbsCS");
478 long q = 1000;
479 ener = make_log_mesh_ec(2.0e-6, 2.0e-1, q);
480 cs.put_qel(q, 0.0);
481 long n;
482 for (n = 0; n < q; n++) {
483 double energy = ener[n];
484 if (energy < threshold)
485 cs[n] = 0.0;
486 else {
487 double Q = 5.5 + l - 0.5 * P;
488 double y = energy / E0;
489 double Fpasc = ((y - 1) * (y - 1) + yw * yw) * pow(y, (-Q)) *
490 pow((1.0 + sqrt(y / ya)), (-P));
491 Fpasc = Fpasc * sigma;
492 cs[n] = Fpasc;
493 //Iprint3n(mcout, energy, E0, threshold);
494 //Iprint4n(mcout, Q, y, sigma, Fpasc);
495 }
496 }
498}
499
501 const SimpleTablePhotoAbsCS& part,
502 double emax_repl) {
503 mfunname("SimpleTablePhotoAbsCS::SimpleTablePhotoAbsCS (const "
504 "SimpleTablePhotoAbsCS& total,...)");
505
506 *this = total; // to assure that all is preserved
507
508 long qe_i = total.ener.get_qel();
509
510 const DynLinArr<double>& ener_r = part.get_arr_ener();
511 const DynLinArr<double>& cs_r = part.get_arr_CS();
512 long qe_r = ener_r.get_qel();
513 BlkArr<double> new_ener; // arrays to grow
514 BlkArr<double> new_cs;
515 long qe = 0;
516 long ne;
517 for (ne = 0; ne < qe_r; ne++) // first write replacements
518 {
519 if (ener_r[ne] >= total.get_threshold() && ener_r[ne] <= emax_repl) {
520 qe++;
521 new_ener.put_qel(qe);
522 new_cs.put_qel(qe);
523 new_ener[qe - 1] = ener_r[ne];
524 new_cs[qe - 1] = cs_r[ne];
525 }
526 }
527 for (ne = 0; ne < qe_i; ne++) {
528 if (ener[ne] >= total.get_threshold() && ener[ne] > emax_repl) {
529 qe++;
530 new_ener.put_qel(qe);
531 new_cs.put_qel(qe);
532 new_ener[qe - 1] = total.ener[ne];
533 new_cs[qe - 1] = total.cs[ne];
534 }
535 }
536 ener.put_qel(qe);
537 cs.put_qel(qe);
538 for (ne = 0; ne < qe; ne++) {
539 ener[ne] = new_ener[ne];
540 cs[ne] = new_cs[ne];
541 //Iprint3n(mcout, ne, ener[ne], cs[ne]);
542
543 }
544}
545
547 long q = ener.get_qel();
548 long ne;
549 long nez;
550 for (ne = 0; ne < q; ne++) {
551 if (cs[ne] > 0.0) break;
552 }
553 if (ne > 0) {
554 long qn = q - ne;
555 DynLinArr<double> enern(qn);
556 DynLinArr<double> csn(qn);
557 for (nez = ne; nez < q; nez++) {
558 enern[nez - ne] = ener[nez];
559 csn[nez - ne] = cs[nez];
560 }
561 ener = enern;
562 cs = csn;
563 }
564}
566 long q = ener.get_qel();
567 long ne;
568 long nez;
569 for (ne = 0; ne < q; ne++) {
570 if (cs[ne] > level) break;
571 }
572 if (ne > 0) {
573 long qn = q - ne;
574 DynLinArr<double> enern(qn);
575 DynLinArr<double> csn(qn);
576 for (nez = ne; nez < q; nez++) {
577 enern[nez - ne] = ener[nez];
578 csn[nez - ne] = cs[nez];
579 }
580 ener = enern;
581 cs = csn;
582 }
583}
584
585double SimpleTablePhotoAbsCS::get_CS(double energy) const {
586 mfunname("double SimpleTablePhotoAbsCS::get_CS(double energy) const");
587 //mcout<<"SimpleTablePhotoAbsCS::get_CS is started\n";
588 //Iprintn(mcout, energy);
589 long q = ener.get_qel();
590 if (q == 0) return 0.0;
591 check_econd11(q, == 1, mcerr);
592 //Iprint3n(mcout, threshold, ener[0], energy);
593 if (energy < threshold)
594 return 0.0;
595 else {
596 if (energy <= ener[q - 1]) {
597 //mcout<<"Create mesh\n";
599 double s;
601 double, DynLinArr<double>,
603 pcmd, cs, &my_val_fun, energy, 1, threshold, 0, DBL_MAX);
604 return s;
605 } else {
606 if (energy == DBL_MAX)
607 return 0.0;
608 else
609 return cs[q - 1] * pow(energy, -2.75) / pow(ener[q - 1], -2.75);
610 }
611
612 /*
613 long i1, i2;
614 long q = ener.get_qel();
615 if(energy < ener[0])
616 {
617 i1 = 0;
618 i2 = 1;
619 return cs[i1] + (energy - ener[i1])*(cs[i2] - cs[i1])/
620 (ener[i2] - ener[i1]);
621 }
622 else
623 {
624 if(energy < ener[q-1])
625 {
626 // now find interval
627 i1 = 0;
628 i2 = q-1;
629 while( i2 - i1 != 1)
630 {
631 long i3 = (i1 + i2)/2;
632 if(energy < ener[i3])
633 i2 = i3;
634 else
635 i1 = i3;
636 }
637 //mcout<<"i12="<<i1<<' '<<i2
638 // <<" ener="<< ener[i1] <<' '<< ener[i2] <<std::endl;
639 // now make interpolation
640 // linear is not always good
641 //return cs[i1] + (energy -
642 // ener[i1])*(cs[i2] - cs[i1])/(ener[i2] - ener[i1]);
643 //if(cs[i2] >= cs[i1] || ener[i2] < 1.2*ener[i1] ||
644 // cs[i1] == 0 || cs[i2] == 0 )
645 if( sign_nonlinear_interpolation(ener[i1], cs[i1],
646 ener[i2], cs[i2], threshold) == 0 )
647 { // linear case
648 //mcout<<"linear\n";
649 //Iprintn(mcout, cs[i1] + (energy - ener[i1])*(cs[i2] -
650 cs[i1])/(ener[i2] - ener[i1]) );
651 return cs[i1] + (energy -
652 ener[i1])*(cs[i2] - cs[i1])/(ener[i2] - ener[i1]);
653 }
654 else
655 {
656 //mcout<<"power\n";
657 //power
658 double pw = log(cs[i1]/cs[i2]) / log(ener[i1]/ener[i2]);
659 //Iprintn(mcout, pw);
660 return cs[i1] * pow( energy , pw ) / pow( ener[i1] , pw );
661 }
662 }
663 else
664 {
665 return cs[q-1] * pow( energy , -2.75 ) / pow( ener[q-1] , -2.75 );
666 }
667 }
668 */
669 }
670}
671
673 double energy2) const {
674 mfunname("double SimpleTablePhotoAbsCS::get_integral_CS(...)");
675
676 long q = ener.get_qel();
677 if (q == 0) return 0.0;
678 check_econd11(q, == 1, mcerr);
679 //Iprint2n(mcout, energy2, threshold);
680 if (energy2 < threshold) return 0.0;
681 if (energy1 < threshold) energy1 = threshold;
682 double s = 0.0;
683 double energy21 = ener[q - 1];
684 if (energy1 < energy21) {
685 if (energy21 > energy2) energy21 = energy2;
686 check_econd12(energy1, >, energy21, mcerr);
687 s = glin_integ_ar(ener, cs, ener.get_qel(), energy1, energy21, threshold);
688 }
689 //print(mcout, 3);
690 //mcout << "energy1="<<energy1
691 // << " energy21="<<energy21
692 // << " ener[q-1]="<<ener[q-1]
693 // << " threshold="<<threshold
694 // << " s="<<s<<'\n';
695 check_econd11(s, < 0.0, mcout);
696 if (energy2 > ener[q - 1]) {
697 // add tail
698 if (energy2 == DBL_MAX) {
699 if (energy1 < ener[q - 1]) energy1 = ener[q - 1];
700 double c =
701 cs[q - 1] / (1.75 * pow(ener[q - 1], -2.75)) * pow(energy1, -1.75);
702 //check_econd11(c , < 0.0, mcout);
703 s += c;
704 } else {
705 if (energy1 < ener[q - 1]) energy1 = ener[q - 1];
706 double c = cs[q - 1] / (1.75 * pow(ener[q - 1], -2.75)) *
707 (pow(energy1, -1.75) - pow(energy2, -1.75));
708 //check_econd11(c , < 0.0, mcout);
709 s += c;
710 }
711 }
712 return s;
713}
714
716 mfunnamep("void SimpleTablePhotoAbsCS::scale(double fact)");
717 long q = ener.get_qel();
718 for (long n = 0; n < q; ++n) {
719 cs[n] *= fact;
720 }
721}
722
723void SimpleTablePhotoAbsCS::print(std::ostream& file, int l) const {
724 if (l > 0) {
725 Ifile << "SimpleTablePhotoAbsCS: name=" << name << " Z = " << Z
726 << std::endl;
727 Ifile << " threshold = " << threshold << " file_name=" << file_name
728 << std::endl;
729 if (l > 1) {
730 indn.n += 2;
731 for (long n = 0; n < ener.get_qel(); ++n) {
732 Ifile << "n=" << n << " ener=" << ener[n] << " cs=" << cs[n]
733 << std::endl;
734 }
735 indn.n -= 2;
736 }
737 }
738}
739
740PhenoPhotoAbsCS::PhenoPhotoAbsCS() : PhotoAbsCS("none", 0, 0.0), power(0.0) {}
741
742PhenoPhotoAbsCS::PhenoPhotoAbsCS(const String& fname, int fZ, double fthreshold,
743 double fpower)
744 : PhotoAbsCS(fname, fZ, fthreshold), power(fpower) {
745 mfunname("PhenoPhotoAbsCS::PhenoPhotoAbsCS");
746 check_econd11a(power, <= 2,
747 " cannot be so, otherwise the integral is infinite", mcerr);
748 factor =
749 pow(threshold, power - 1.) * Thomas_sum_rule_const_Mb * Z * (power - 1.);
750}
751
752double PhenoPhotoAbsCS::get_CS(double energy) const {
753 if (energy < threshold || energy == DBL_MAX) return 0.0;
754 return factor * (pow(energy, -power));
755}
756
757double PhenoPhotoAbsCS::get_integral_CS(double energy1, double energy2) const {
758 //Imcout<<"energy1="<<energy1<<" energy2="<<energy2<<'\n';
759 if (energy2 < threshold) return 0.0;
760 if (energy1 < threshold) energy1 = threshold;
761 double s;
762 if (energy2 == DBL_MAX) {
763 s = factor / (power - 1.) * (1. / pow(energy1, power - 1.));
764 } else {
765 s = factor / (power - 1.) *
766 (1. / pow(energy1, power - 1.) - 1. / pow(energy2, power - 1.));
767 }
768 return s;
769}
770
771void PhenoPhotoAbsCS::scale(double fact) {
772 mfunnamep("void PhenoPhotoAbsCS::scale(double fact)");
773 factor *= fact;
774}
775
776void PhenoPhotoAbsCS::print(std::ostream& file, int l) const {
777 if (l > 0) {
778 Ifile << "PhenoPhotoAbsCS: name=" << name << " Z = " << Z << std::endl;
779 Ifile << " threshold = " << threshold << " power=" << power
780 << " factor=" << factor << std::endl;
781 }
782}
783
784/*
785// drafts, unused
786mesh_for_FitBT = EnergyMesh(
787
788FitBTPhotoAbsCS::FitBTPhotoAbsCS(void): PhotoAbsCS("none",0, 0.0),
789 E0(0.0), yw(0.0), ya(0.0), P(0.0), sigma(0.0)
790{;}
791
792FitBTPhotoAbsCS::FitBTPhotoAbsCS
793(const String& fname, int fZ, double fthreshold,
794 int fl, double fE0, double fyw, double fya,
795 double fP, double fsigma):
796 PhotoAbsCS(fname, fZ, fthreshold),
797 l(fl), E0(fE0), yw(fyw), ya(fya), P(fP), sigma(fsigma)
798{
799 mfunname("FitBTPhotoAbsCS::FitBTPhotoAbsCS(...)");
800}
801
802double FitBTPhotoAbsCS::get_CS(double energy) const
803{
804 if(energy < threshold)
805 return 0.0;
806 else
807 {
808 if(energy == BDL_MAX)
809 return 0.0;
810 else
811 {
812 double Q = 5.5 + fl - 0.5 * fP;
813 double y = energy / fE0;
814 double Fpasc = ((y - 1) * (y - 1) + yw * yw)
815 * pow(y , (-Q) ) * pow((1.0 + sqrt(y / ya)) , (-P));
816 Fpasc = Fpasc * sigma0;
817 return Fpasc;
818 }
819}
820
821double PhenoPhotoAbsCS::get_integral_CS(double energy1, double energy2) const
822{
823 //Imcout<<"energy1="<<energy1<<" energy2="<<energy2<<'\n';
824 if(energy2 < threshold)
825 return 0.0;
826 else
827 {
828 if(energy1 < threshold) energy1 = threshold;
829 double s = factor / (power - 1.0) *
830 (1.0/pow (energy1 , power-1.0) - 1.0/pow (energy2 , power-1.0));
831 //Imcout<<"s="<<s<<'\n';
832 return s;
833 }
834}
835*/
836
837/*
838FitBTPhotoAbsCS::FitBTPhotoAbsCS
839(const String& fname, int fZ, double fthreshold, double fpower):
840 PhotoAbsCS(fname, fZ, 0.0)
841{
842 mfunname("FitBTPhotoAbsCS::FitBTPhotoAbsCS");
843
844 name = fname;
845#ifdef USE_STLSTRING
846 std::ifstream BT_file(name.c_str());
847#else
848 std::ifstream BT_file(name);
849#endif
850 if( !BT_file )
851 {
852 funnw.ehdr(mcerr);
853 mcerr<<"cannot open file "<<name<<std::endl;
854 spexit(mcerr);
855 }
856 int i = 0;
857 while( (i = findmark(BT_file, "$")) == 1 )
858 {
859 long iZ;
860 BT_file >> iZ;
861 if(iz == Z)
862 {
863
864
865 funnw.ehdr(mcerr);
866 mcerr<<"there is no element Z="<<fZ<<" in file
867"<<threshold_file_name<<std::endl;
868 spexit(mcerr);
869 mark1:
870
871*/
872
873//------------------------------------------------------------------------
874
876 double fchannel_prob_dens, const DynLinArr<double>& felectron_energy,
877 const DynLinArr<double>& fphoton_energy, int s_all_rest) {
878 mfunnamep("void AtomicSecondaryProducts::add_channel(...)");
879 check_econd21(fchannel_prob_dens, < 0.0 ||, > 1.0, mcerr);
880 long q_old = channel_prob_dens.get_qel();
881 long q_new = q_old + 1;
884 photon_energy.put_qel(q_new);
885 //Iprintn(mcout, s);
886 //check_econd11( s , > 1.0 , mcerr);
887 if (s_all_rest == 1) {
888 double s = 0.0;
889 for (long n = 0; n < q_old; ++n) {
890 s += channel_prob_dens[n];
891 }
892 check_econd21(s, < 0.0 ||, > 1.0, mcerr);
893 fchannel_prob_dens = 1.0 - s;
894 }
895 channel_prob_dens[q_old] = fchannel_prob_dens;
896 electron_energy[q_old] = felectron_energy;
897 photon_energy[q_old] = fphoton_energy;
898 double s = 0.0;
899 for (long n = 0; n < q_new; ++n) {
900 s += channel_prob_dens[n];
901 }
902 if (s > 1.0) {
903 funnw.ehdr(mcerr);
904 mcerr << "s > 1.0, s=" << s << '\n';
905 Iprintn(mcerr, q_new);
906 for (long n = 0; n < q_new; ++n) {
907 mcerr << "n=" << n << " channel_prob_dens[n]=" << channel_prob_dens[n]
908 << '\n';
909 }
910 spexit(mcerr);
911 }
912}
913
915 DynLinArr<double>& felectron_energy,
916 DynLinArr<double>& fphoton_energy) const {
917 mfunname("int AtomicSecondaryProducts::get_channel(...)");
918#ifdef DEBUG_PRINT_get_escape_particles
919 mcout << "AtomicSecondaryProducts::get_channel is started\n";
921#endif
922 int ir = 0;
923 if (channel_prob_dens.get_qel() > 0) {
924 double rn = SRANLUX();
925#ifdef DEBUG_PRINT_get_escape_particles
926 Iprintn(mcout, rn);
927#endif
928 if (channel_prob_dens.get_qel() == 1) {
929 if (rn < channel_prob_dens[0]) {
930 felectron_energy = electron_energy[0];
931 fphoton_energy = photon_energy[0];
932 ir = 1;
933 }
934 } else {
935 long q = channel_prob_dens.get_qel();
936 double s = 0.0;
937 for (long n = 0; n < q; ++n) {
938 s += channel_prob_dens[n];
939 if (rn <= s) {
940 felectron_energy = electron_energy[n];
941 fphoton_energy = photon_energy[n];
942 ir = 1;
943#ifdef DEBUG_PRINT_get_escape_particles
944 Iprint2n(mcout, n, s);
945#endif
946 break;
947 }
948 }
949 }
950 }
951#ifdef DEBUG_PRINT_get_escape_particles
952 mcout << "AtomicSecondaryProducts::get_channel is finishing\n";
953 Iprintn(mcout, ir);
954#endif
955 return ir;
956}
957
958void AtomicSecondaryProducts::print(std::ostream& file, int l) const {
959 if (l > 0) {
960 Ifile << "AtomicSecondaryProducts(l=" << l << "):\n";
961 long q = channel_prob_dens.get_qel();
962 Ifile << "number of channels=" << q << '\n';
963 indn.n += 2;
964 for (long n = 0; n < q; ++n) {
965 Ifile << "n_channel=" << n << " probability=" << channel_prob_dens[n]
966 << '\n';
967 indn.n += 2;
968 long qel = electron_energy[n].get_qel();
969 Ifile << "number of electrons=" << qel << '\n';
970 indn.n += 2;
971 for (long nel = 0; nel < qel; ++nel) {
972 Ifile << "nel=" << nel << " electron_energy=" << electron_energy[n][nel]
973 << '\n';
974 }
975 indn.n -= 2;
976 long qph = photon_energy[n].get_qel();
977 Ifile << "number of photons=" << qph << '\n';
978 indn.n += 2;
979 for (long nph = 0; nph < qph; ++nph) {
980 Ifile << "nph=" << nph << " photon_energy=" << photon_energy[n][nph]
981 << '\n';
982 }
983 indn.n -= 2;
984 indn.n -= 2;
985 }
986 indn.n -= 2;
987 }
988}
989
990AtomPhotoAbsCS::AtomPhotoAbsCS() : name("none"), Z(0), qshell(0) {}
991
992double AtomPhotoAbsCS::get_TICS(double energy,
993 double factual_minimal_threshold) const {
994 mfunname("double AtomPhotoAbsCS::get_TICS(double energy, double "
995 "factual_minimal_threshold) const");
996
997 if (factual_minimal_threshold <= energy) {
998 // All what is absorbed, should ionize
999 return get_ACS(energy);
1000 }
1001 return 0.0;
1002}
1003
1005 double energy1, double energy2, double factual_minimal_threshold) const {
1006 mfunname("double AtomPhotoAbsCS::get_integral_TICS(double energy1, double "
1007 "energy2, double factual_minimal_threshold) const");
1008
1009 if (factual_minimal_threshold <= energy2) {
1010 // All what is absorbed, should ionize
1011 if (energy1 < factual_minimal_threshold) {
1012 energy1 = factual_minimal_threshold;
1013 }
1014 return get_integral_ACS(energy1, energy2);
1015 }
1016 return 0.0;
1017}
1018
1019double AtomPhotoAbsCS::get_TICS(int nshell, double energy,
1020 double factual_minimal_threshold) const {
1021 mfunname("double AtomPhotoAbsCS::get_TICS(int nshell, double energy, double "
1022 "factual_minimal_threshold) const");
1023
1024 if (s_ignore_shell[nshell] == 0) {
1025 if (factual_minimal_threshold <= energy) {
1026 // All what is absorbed, should ionize
1027 return get_integral_ACS(nshell, energy);
1028 }
1029 }
1030 return 0.0;
1031}
1032
1034 int nshell, double energy1, double energy2,
1035 double factual_minimal_threshold) const {
1036 mfunname("double AtomPhotoAbsCS::get_integral_TICS(int nshell, double "
1037 "energy1, double energy2, double factual_minimal_threshold) const");
1038
1039 if (s_ignore_shell[nshell] == 0) {
1040 if (factual_minimal_threshold <= energy1) {
1041 // All what is absorbed, should ionize
1042 return get_integral_ACS(nshell, energy1, energy2);
1043 }
1044 if (factual_minimal_threshold >= energy2) return 0.0;
1045 return get_integral_ACS(nshell, factual_minimal_threshold, energy2);
1046 }
1047 return 0.0;
1048}
1049
1051 mfunname("void AtomPhotoAbsCS::remove_shell(int nshell)");
1052 check_econd21(nshell, < 0 ||, >= qshell, mcerr);
1053 s_ignore_shell[nshell] = 1;
1054}
1055
1057 mfunname("void AtomPhotoAbsCS::restore_shell(int nshell)");
1058 check_econd21(nshell, < 0 ||, >= qshell, mcerr);
1059 s_ignore_shell[nshell] = 0;
1060}
1061
1062void AtomPhotoAbsCS::print(std::ostream& file, int l) const {
1063 mfunnamep("void AtomPhotoAbsCS::print(std::ostream& file, int l) const");
1064 if (l > 0) {
1065 Ifile << "AtomPhotoAbsCS(l=" << l << "): name=" << name << " Z = " << Z
1066 << " qshell = " << qshell << std::endl;
1067 Iprintn(mcout, asp.get_qel());
1068 long q = asp.get_qel();
1069 if (q == 0) {
1070 q = s_ignore_shell.get_qel();
1071 indn.n += 2;
1072 for (long n = 0; n < q; ++n) {
1073 Ifile << "n=" << n << " s_ignore_shell[n] = " << s_ignore_shell[n]
1074 << '\n';
1075 }
1076 indn.n -= 2;
1077 } else {
1078 check_econd12(asp.get_qel(), !=, s_ignore_shell.get_qel(), mcerr);
1079 indn.n += 2;
1080 for (long n = 0; n < q; ++n) {
1081 Ifile << "n=" << n << " s_ignore_shell[n] = " << s_ignore_shell[n]
1082 << '\n';
1083 asp[n].print(mcout, l);
1084 }
1085 indn.n -= 2;
1086 }
1087 }
1088}
1089
1090std::ostream& operator<<(std::ostream& file, const AtomPhotoAbsCS& f) {
1091 f.print(file, 1);
1092 return file;
1093}
1094
1095double AtomPhotoAbsCS::get_I_min(void) const {
1096 mfunname("double AtomPhotoAbsCS::get_I_min(void) const");
1097 double st = DBL_MAX;
1098 for (int n = 0; n < qshell; ++n) {
1099 // currently the minimal shell is the last,
1100 // but to avoid this assumption we check all
1101 if (get_threshold(n) < st) st = get_threshold(n);
1102 }
1103 return st;
1104}
1105
1106void AtomPhotoAbsCS::get_escape_particles(int nshell, double energy,
1107 DynLinArr<double>& el_energy,
1108 DynLinArr<double>& ph_energy) const {
1109 mfunname("void AtomPhotoAbsCS::get_escape_particles(...)");
1110#ifdef DEBUG_PRINT_get_escape_particles
1111 mcout << "AtomPhotoAbsCS::get_escape_particles is started\n";
1112 // to find minimal shell
1113 Iprintn(mcout, nshell);
1114 Iprintn(mcout, energy);
1115#endif
1116 check_econd12(energy, <, 0.5 * get_threshold(nshell), mcerr);
1117 // In principle, the energy is allowed to be slightly less than threshold
1118 // due to unprecision of definition of point-wise cross sections.
1119 // To keep correct norm it is better not to ignore such events.
1120 // They usually can be treated quite well.
1121 // The factor 0.5 is put there just as arbitrary check for full stupidity.
1122
1123 // program is complicated; to make sure that there is no remains
1124 // from previous runs are allowed to pass, I clear arrays:
1125 el_energy = DynLinArr<double>(0);
1126 ph_energy = DynLinArr<double>(0);
1127
1128 int n_min = 0;
1129 double st = DBL_MAX;
1130 for (int n = 0; n < qshell; ++n) {
1131 // currently the minimal shell is the last,
1132 // but to avoid this assumption we check all
1133 if (get_threshold(n) < st) {
1134 n_min = n;
1135 st = get_threshold(n);
1136 }
1137 }
1138#ifdef DEBUG_PRINT_get_escape_particles
1139 Iprintn(mcout, n_min);
1140#endif
1141 if (nshell == n_min) {
1142 // generate delta-electron, here it will be the only one, since the shell
1143 // is the outmost valent one.
1144 double en;
1145 if (energy - get_threshold(n_min) >= 0.0) {
1146 en = energy - get_threshold(n_min);
1147 } else {
1148 en = 0.0;
1149 }
1150 el_energy = DynLinArr<double>(1, en);
1151 ph_energy = DynLinArr<double>(0);
1152 } else {
1153 // Energy of photo-electron
1154 double en = energy - get_threshold(nshell);
1155 double hdist = 0.0; // used to preserve the balance of energy
1156 // virtual gamma are generated by energy mesh
1157 // and their energy could be little less than the shell energy.
1158 // To avoid generation of electrons with negative energy
1159 // their energy is corrected. The value of correction is hdist.
1160 // To preserve energy this hdist is then distributed over
1161 // the other secondary products if they exist.
1162 if (en < 0.0) {
1163 hdist = -en;
1164 en = 0.0;
1165 }
1166 int is = 0;
1167 DynLinArr<double> felectron_energy;
1168 DynLinArr<double> fphoton_energy;
1169#ifdef DEBUG_PRINT_get_escape_particles
1170 Iprint2n(mcout, asp.get_qel(), get_qshell());
1171#endif
1172#ifndef DEBUG_ignore_non_standard_channels
1173 if (asp.get_qel() == get_qshell()) {
1174 // works only in this case?
1175 is = asp[nshell].get_channel(felectron_energy, fphoton_energy);
1176 // Here zero can be if the shell is not included in database
1177 // or if not standard channel is not choosen by random way.
1178 // In both cases the standard way should be invoked.
1179 }
1180#endif
1181 int main_n = get_main_shell_number(nshell);
1182#ifdef DEBUG_PRINT_get_escape_particles
1183 Iprint2n(mcout, nshell, main_n);
1184 Iprintn(mcout, is);
1185 Iprint(mcout, felectron_energy);
1186 Iprint(mcout, fphoton_energy);
1187#endif
1188 //check_econd22( is , == 0 && , main_n , < 0 , mcerr); // since
1189 //// there will be neither default channel nor the explicit one
1190 //if(is == 0 && main_n > 0)
1191 if (is == 0) {
1192 // generate default channel
1193 if (main_n > 0) {
1194 // first find the principal quantum number of the deepest shell
1195 int main_n_largest = 0;
1196 for (int n = 0; n < qshell; ++n) {
1197 int t = get_main_shell_number(n);
1198 if (t > main_n_largest) main_n_largest = t;
1199 }
1200#ifdef DEBUG_PRINT_get_escape_particles
1201 Iprintn(mcout, main_n_largest);
1202#endif
1203 if (main_n_largest - main_n >= 2) {
1204 // At least K, l, M shells exists
1205 // In this case we use more advanced scheme
1206 // Look for shell with larger main number and with less energy
1207 int n_choosen = -1;
1208 double thr = DBL_MAX; // this will be the least threshold
1209 // among the shells with next principal number
1210 for (int n = 0; n < qshell; ++n) {
1211 // currently the minimal shell is the last,
1212 // but to avoid this assumption we check all
1213 int main_n_t = get_main_shell_number(n);
1214 if (main_n_t > 0 && main_n_t == main_n + 1) {
1215 if (thr > get_threshold(n)) {
1216 n_choosen = n;
1217 thr = get_threshold(n);
1218 }
1219 }
1220 }
1221#ifdef DEBUG_PRINT_get_escape_particles
1222 Iprint2n(mcout, n_choosen, thr);
1223#endif
1224 check_econd11(n_choosen, < 0, mcerr);
1225 double e_temp = 0.0;
1226 if ((e_temp = get_threshold(nshell) - hdist -
1227 2.0 * get_threshold(n_choosen))
1228 //if(get_threshold(nshell) - hdist - 2.0*get_threshold(n_choosen)
1229 > 0.0) {
1230 el_energy = DynLinArr<double>(4);
1231 el_energy[0] = en; // photo-electron
1232 el_energy[1] = e_temp;
1233 //el_energy[1] = get_threshold(nshell) - hdist -
1234 // 2.0*get_threshold(n_choosen);
1235 // first Auger from choosen shell
1236 // Then filling two vacancies at the next (choosen) shell
1237 // from the outermost one
1238 if ((e_temp = get_threshold(n_choosen) -
1239 2.0 * get_threshold(n_min)) > 0.0) {
1240 el_energy[2] = e_temp;
1241 //el_energy[2] =
1242 // get_threshold(n_choosen) - 2.0 * get_threshold(n_min);
1243 el_energy[3] = el_energy[2];
1244 check_econd11(el_energy[2], < 0.0, mcerr);
1245 } else { // ignore two last Auger electrons
1246 el_energy.put_qel(2);
1247 }
1248 } else if ((e_temp = get_threshold(nshell) - hdist -
1249 get_threshold(n_choosen) - get_threshold(n_min)) > 0.0) {
1250 el_energy = DynLinArr<double>(3);
1251 el_energy[0] = en; // photo-electron
1252 el_energy[1] = e_temp;
1253 //el_energy[1] = get_threshold(nshell) - hdist -
1254 // get_threshold(n_choosen) - get_threshold(n_min);
1255 // filling initially ionized level from choosen
1256 // and emmitance of Auger from outermost.
1257 check_econd11(el_energy[1], < 0.0, mcerr);
1258 if ((e_temp = get_threshold(n_choosen) -
1259 2.0 * get_threshold(n_min)) > 0.0) {
1260 el_energy[2] = e_temp;
1261 //el_energy[2] = get_threshold(n_choosen) - 2.0 *
1262 //get_threshold(n_min);
1263 } else { // ignore this last Auger electron
1264 el_energy.put_qel(2);
1265 }
1266 }
1267 ph_energy = DynLinArr<double>(0);
1268 } else // for if(main_n_largest - main_n >= 2)
1269 {
1270 //generate Auger from the outermost shell
1271 double en1 =
1272 get_threshold(nshell) - hdist - 2.0 * get_threshold(n_min);
1273 if (en1 >= 0.0) {
1274 el_energy = DynLinArr<double>(2);
1275 el_energy[0] = en;
1276 el_energy[1] = en1;
1277 } else {
1278 el_energy = DynLinArr<double>(1);
1279 el_energy[0] = en;
1280 }
1281 ph_energy = DynLinArr<double>(0);
1282 }
1283 } else // for if(main_n > 0)
1284 { // principal numbers are not available, then we generate Auger
1285 // to outmost shell
1286 double en1 = get_threshold(nshell) - hdist - 2.0 * get_threshold(n_min);
1287 if (en1 >= 0.0) {
1288 el_energy = DynLinArr<double>(2);
1289 el_energy[0] = en;
1290 el_energy[1] = en1;
1291 } else {
1292 el_energy = DynLinArr<double>(1);
1293 el_energy[0] = en;
1294 }
1295 ph_energy = DynLinArr<double>(0);
1296 }
1297 } else // for if(is == 0), generate explicit channel
1298 {
1299 // Generate photo-electron and
1300 // just copy all what is proposed by get_channel
1301 // with corrections by hdist.
1302
1303 el_energy = DynLinArr<double>(1 + felectron_energy.get_qel());
1304 el_energy[0] = en;
1305 long q = felectron_energy.get_qel();
1306 for (long n = 0; n < q; ++n) {
1307 check_econd21(felectron_energy[n], < 0 ||, > get_threshold(nshell),
1308 mcerr);
1309 el_energy[1 + n] = felectron_energy[n] - hdist;
1310 if (el_energy[1 + n] < 0) {
1311 hdist = -el_energy[1 + n];
1312 el_energy[1 + n] = 0.0;
1313 } else {
1314 hdist = 0.0;
1315 }
1316 }
1317 ph_energy = DynLinArr<double>(fphoton_energy.get_qel());
1318 q = fphoton_energy.get_qel();
1319 for (long n = 0; n < q; ++n) {
1320 check_econd21(fphoton_energy[n], < 0 ||, > get_threshold(nshell),
1321 mcerr);
1322 ph_energy[n] = fphoton_energy[n] - hdist;
1323 if (ph_energy[n] < 0) {
1324 hdist = -ph_energy[n];
1325 ph_energy[n] = 0.0;
1326 } else {
1327 hdist = 0.0;
1328 }
1329 }
1330 }
1331 }
1332#ifdef DEBUG_PRINT_get_escape_particles
1333 Iprintn(mcout, ph_energy);
1334 Iprintn(mcout, el_energy);
1335 mcout << "AtomPhotoAbsCS::get_escape_particles: exitting\n";
1336#endif
1337}
1338
1340 mfunnamep("AtomicSecondaryProducts* AtomPhotoAbsCS::get_asp(int nshell)");
1341 check_econd21(nshell, < 0 ||, >= qshell, mcerr);
1342 return &(asp[nshell]);
1343}
1344
1346
1348 : file_name(ffile_name) {
1349 mfunnamep("SimpleAtomPhotoAbsCS::SimpleAtomPhotoAbsCS(int fZ, const String& "
1350 "ffile_name)");
1351 check_econd11(fZ, < 1, mcerr);
1352#ifdef USE_STLSTRING
1353 std::ifstream file(file_name.c_str());
1354#else
1355 std::ifstream file(file_name);
1356#endif
1357 if (!file) {
1358 funnw.ehdr(mcerr);
1359 mcerr << "cannot open file " << file_name << std::endl;
1360 spexit(mcerr);
1361 }
1362 while (findmark(file, "#") == 1) {
1363 file >> Z;
1364 if (Z == fZ) {
1365 file >> qshell;
1366 check_econd21(qshell, < 1 ||, > 10000, mcerr);
1368 file >> name;
1372 int sZshell = 0;
1373 for (int nshell = 0; nshell < qshell; ++nshell) {
1374 double thr = 0.0;
1375 int Zshell = 0;
1376 //double fl;
1377 String shell_name;
1378 file >> thr;
1379 check_econd11(thr, <= 0.0, mcerr);
1380 file >> Zshell;
1381 check_econd11(Zshell, <= 0, mcerr);
1382 sZshell += Zshell;
1383 file >> fl[nshell];
1384 findmark(file, "!");
1385 file >> shell_name;
1386 //PhotoAbsCS* ap = new PhenoPhotoAbsCS(shell_name, Zshell, thr);
1387 //acs[nshell].pass( ap );
1388 acs[nshell].pass(new PhenoPhotoAbsCS(shell_name, Zshell, thr * 1.0e-6));
1389 }
1390 check_econd12(sZshell, !=, Z, mcerr);
1391
1392 int n_min = 0;
1393 double st = DBL_MAX;
1394 for (int nshell = 0; nshell < qshell; ++nshell) {
1395 // currently the minimal shell is the last,
1396 // but to avoid this assumption we check all
1397 if (get_threshold(nshell) < st) n_min = nshell;
1398 }
1399 for (int nshell = 0; nshell < qshell; ++nshell) {
1400 if (fl[nshell] > 0) {
1401 check_econd12(nshell, ==, n_min, mcerr);
1402 DynLinArr<double> felectron_energy(0);
1403 DynLinArr<double> fphoton_energy(1);
1404 fphoton_energy[0] = get_threshold(nshell) - get_threshold(n_min);
1405 asp[nshell].add_channel(fl[nshell], felectron_energy, fphoton_energy);
1406 }
1407 }
1408 return;
1409 }
1410 }
1411 funnw.ehdr(mcerr);
1412 mcerr << "there is no element Z=" << fZ << " in file " << file_name
1413 << std::endl;
1414 spexit(mcerr);
1415}
1416
1417SimpleAtomPhotoAbsCS::SimpleAtomPhotoAbsCS(int fZ, const PhotoAbsCS& facs) {
1418 mfunname("SimpleAtomPhotoAbsCS::SimpleAtomPhotoAbsCS(int fZ, const "
1419 "PhotoAbsCS& facs)");
1420 check_econd11(fZ, <= 0, mcerr);
1421 check_econd12(fZ, !=, facs.get_Z(), mcerr);
1422 Z = fZ;
1423 qshell = 1;
1425 name = facs.get_name();
1426 acs.put_qel(1);
1427 acs[0].put(&facs);
1428}
1429
1430double SimpleAtomPhotoAbsCS::get_threshold(int nshell) const {
1431 mfunname("double SimpleAtomPhotoAbsCS::get_threshold(int nshell) const");
1432 check_econd21(nshell, < 0 ||, > qshell, mcerr);
1433 return acs[nshell]->get_threshold();
1434}
1435
1436double SimpleAtomPhotoAbsCS::get_ACS(double energy) const {
1437 mfunname("double SimpleAtomPhotoAbsCS::get_ACS(double energy) const");
1438 double s = 0.0;
1439 for (int n = 0; n < qshell; ++n) {
1440 if (s_ignore_shell[n] == 0) {
1441 s += acs[n]->get_CS(energy);
1442 }
1443 }
1444 return s;
1445}
1447 double energy2) const {
1448 mfunnamep(
1449 "double SimpleAtomPhotoAbsCS::get_integral_ACS(double energy)const");
1450 double s = 0.0;
1451 for (int n = 0; n < qshell; ++n) {
1452 if (s_ignore_shell[n] == 0) {
1453 double t;
1454 s += (t = acs[n]->get_integral_CS(energy1, energy2));
1455 if (t < 0) {
1456 funnw.ehdr(mcout);
1457 mcout << "t < 0\n";
1458 Iprintn(mcout, t);
1459 print(mcout, 4);
1460 spexit(mcout);
1461 }
1462 }
1463 //check_econd11a(t , < 0 , "n="<<n<<'\n', mcerr);
1464 }
1465 return s;
1466}
1467
1468double SimpleAtomPhotoAbsCS::get_ACS(int nshell, double energy) const {
1469 mfunname("double SimpleAtomPhotoAbsCS::get_ACS(int nshell, double energy)");
1470 check_econd21(nshell, < 0 ||, > qshell, mcerr);
1471 if (s_ignore_shell[nshell] == 0) {
1472 return acs[nshell]->get_CS(energy);
1473 }
1474 return 0.0;
1475}
1476
1477double SimpleAtomPhotoAbsCS::get_integral_ACS(int nshell, double energy1,
1478 double energy2) const {
1479 mfunname("double SimpleAtomPhotoAbsCS::get_integral_ACS(int nshell, double "
1480 "energy1, double energy2)");
1481 check_econd21(nshell, < 0 ||, > qshell, mcerr);
1482 if (s_ignore_shell[nshell] == 0) {
1483 return acs[nshell]->get_integral_CS(energy1, energy2);
1484 }
1485 return 0.0;
1486}
1487
1488double SimpleAtomPhotoAbsCS::get_ICS(double energy) const {
1489 mfunname("double SimpleAtomPhotoAbsCS::get_ICS(double energy) const");
1490 double s = 0.0;
1491 for (int n = 0; n < qshell; ++n) {
1492 if (s_ignore_shell[n] == 0) {
1493 s += acs[n]->get_CS(energy);
1494 }
1495 }
1496 return s;
1497}
1498
1500 double energy2) const {
1501 mfunname("double SimpleAtomPhotoAbsCS::get_integral_ICS(double energy1, "
1502 "double energy2) const");
1503 double s = 0.0;
1504 for (int n = 0; n < qshell; ++n) {
1505 if (s_ignore_shell[n] == 0) {
1506 s += acs[n]->get_integral_CS(energy1, energy2);
1507 }
1508 }
1509 return s;
1510}
1511
1512double SimpleAtomPhotoAbsCS::get_ICS(int nshell, double energy) const {
1513 mfunname("double SimpleAtomPhotoAbsCS::get_ICS(int nshell, double energy)");
1514 check_econd21(nshell, < 0 ||, > qshell, mcerr);
1515 if (s_ignore_shell[nshell] == 0) {
1516 return acs[nshell]->get_CS(energy);
1517 }
1518 return 0.0;
1519}
1520
1521double SimpleAtomPhotoAbsCS::get_integral_ICS(int nshell, double energy1,
1522 double energy2) const {
1523 mfunname("double SimpleAtomPhotoAbsCS::get_integral_ICS(int nshell, double "
1524 "energy1, double energy2)");
1525 check_econd21(nshell, < 0 ||, > qshell, mcerr);
1526 if (s_ignore_shell[nshell] == 0) {
1527 return acs[nshell]->get_integral_CS(energy1, energy2);
1528 }
1529 return 0.0;
1530}
1531
1532void SimpleAtomPhotoAbsCS::print(std::ostream& file, int l) const {
1533 if (l > 0) {
1534 Ifile << "SimpleAtomPhotoAbsCS(l=" << l << "): name=" << name
1535 << " Z = " << Z << " qshell = " << qshell
1536 << " file_name=" << file_name << std::endl;
1537 l--;
1538 if (l > 0) {
1539 indn.n += 2;
1540 for (int n = 0; n < qshell; ++n) {
1541 Ifile << "nshell=" << n << std::endl;
1542 acs[n].print(file, l);
1543 }
1544 AtomPhotoAbsCS::print(file, l);
1545 indn.n -= 2;
1546 }
1547 }
1548}
1549
1551 mfunname("int SimpleAtomPhotoAbsCS::get_main_shell_number(int nshell) const");
1552 String shell_name = acs[nshell]->get_name();
1553#ifdef STRSTREAM_AVAILABLE
1554#ifdef USE_STLSTRING
1555 istrstream sfile(shell_name.c_str());
1556#else
1557 istrstream sfile(shell_name);
1558#endif
1559#else
1560 std::istringstream sfile(shell_name.c_str());
1561#endif
1562 int i = -100;
1563 sfile >> i;
1564 //Iprintn(mcout, i);
1565 //check_econd(i < 1 || i > 50 , "i="<<i<<" shell_name="<<shell_name<<'\n' ,
1566 // mcerr);
1567 if (i < 1 || i > 50) {
1568 i = -1;
1569 }
1570 return i;
1571}
1572
1573//----------------------------------------------------------------------
1574
1575ExAtomPhotoAbsCS::ExAtomPhotoAbsCS(int fZ, const String& fthreshold_file_name,
1576 const String& fsimple_table_file_name,
1577 const String& fname,
1578 double fminimal_threshold)
1579 : threshold_file_name(fthreshold_file_name),
1580 simple_table_file_name(fsimple_table_file_name),
1581 BT_file_name("none"),
1582 minimal_threshold(fminimal_threshold) {
1583 mfunnamep("SimpleAtomPhotoAbsCS::SimpleAtomPhotoAbsCS(int fZ, const String& "
1584 "fname, , const String& fthreshold_file_name, const String& "
1585 "fsimple_table_file_name, double fminimal_threshold)");
1586 //Imcout<<"ExAtomPhotoAbsCS::ExAtomPhotoAbsCS is run for fZ="<<fZ<<std::endl;
1587 check_econd11(fZ, < 1, mcerr);
1588#ifdef USE_STLSTRING
1589 std::ifstream threshold_file(threshold_file_name.c_str());
1590#else
1591 std::ifstream threshold_file(threshold_file_name);
1592#endif
1593 if (!threshold_file) {
1594 funnw.ehdr(mcerr);
1595 mcerr << "cannot open file " << threshold_file_name << std::endl;
1596 spexit(mcerr);
1597 }
1598 DynLinArr<double> thr(0, 0.0);
1599 DynLinArr<int> Zshell(0, 0);
1600 DynLinArr<double> fl(0, 0.0);
1601 DynLinArr<String> shell_name(0);
1602 while (findmark(threshold_file, "#") == 1) {
1603 threshold_file >> Z;
1604 if (Z == fZ) {
1605 threshold_file >> qshell;
1606 check_econd21(qshell, < 1 ||, > 10000, mcerr);
1608 //Iprintn(mcout, qshell);
1609 thr = DynLinArr<double>(qshell, 0.0);
1610 Zshell = DynLinArr<int>(qshell, 0);
1611 fl = DynLinArr<double>(qshell, 0.0);
1612 shell_name = DynLinArr<String>(qshell);
1615 String temp_name;
1616 threshold_file >> temp_name;
1617 if (fname == "none")
1618 name = temp_name;
1619 else
1620 name = fname;
1621 int nshell;
1622 int sZshell = 0;
1623 for (nshell = 0; nshell < qshell; nshell++) {
1624 threshold_file >> thr[nshell];
1625 check_econd11(thr[nshell], <= 0.0, mcerr);
1626 thr[nshell] *= 1.0e-6;
1627 threshold_file >> Zshell[nshell];
1628 check_econd11(Zshell[nshell], <= 0, mcerr);
1629 sZshell += Zshell[nshell];
1630 threshold_file >> fl[nshell];
1631 findmark(threshold_file, "!");
1632 threshold_file >> shell_name[nshell];
1633 //PhotoAbsCS* ap = new PhenoPhotoAbsCS(shell_name, Zshell, thr);
1634 //acs[nshell].pass( ap );
1635 }
1636 check_econd12(sZshell, !=, Z, mcerr);
1637 int n_min = 0;
1638 double st = DBL_MAX;
1639 for (nshell = 0; nshell < qshell;
1640 nshell++) { // currently the minimal shell is the last,
1641 // but to avoid this assumption
1642 // we check all
1643 if (thr[nshell] < st) {
1644 n_min = nshell;
1645 st = thr[nshell];
1646 }
1647 }
1648 for (nshell = 0; nshell < qshell; nshell++) {
1649 if (fl[nshell] > 0) {
1650 check_econd12(nshell, ==, n_min, mcerr);
1651 DynLinArr<double> felectron_energy(0);
1652 DynLinArr<double> fphoton_energy(1);
1653 fphoton_energy[0] = thr[nshell] - thr[n_min];
1654 asp[nshell].add_channel(fl[nshell], felectron_energy, fphoton_energy);
1655 }
1656 }
1657 goto mark1;
1658 }
1659 }
1660 funnw.ehdr(mcerr);
1661 mcerr << "there is no element Z=" << fZ << " in file " << threshold_file_name
1662 << std::endl;
1663 spexit(mcerr);
1664mark1:
1665 // Here it reads the PACS as an one shell curve:
1666 SimpleTablePhotoAbsCS stpacs(name, Z, 0.0, fsimple_table_file_name);
1667 const DynLinArr<double>& ener = stpacs.get_arr_ener();
1668 const DynLinArr<double>& CS = stpacs.get_arr_CS();
1669 //Iprintn(mcout, ener.get_qel());
1670 DynLinArr<double> left_CS = CS; // used in sequencial algorithm
1671 DynLinArr<DynLinArr<double> > SCS(qshell); // here cs is saved
1672 int n;
1673 for (n = 0; n < qshell; n++)
1674 SCS[n] = DynLinArr<double>(ener.get_qel(), 0.0);
1675 /*
1676 if(ener[0] < thr[qshell-1])
1677 {
1678 funnw.ehdr(cerr);
1679 mcerr<<"ener[0] < thr[qshell-1]\n";
1680 mcerr<<"This possibility is not bad, but it is not implemented \n"
1681 <<"in the program yet, sorry\n";
1682 spexit(mcerr);
1683 }
1684 */
1685 int nct = qshell - 1; //"current" threshold index
1686 long nce = 0; // "current" energy index
1687 // Let uss ignore values below the lowest threshold
1688 // It is not clear whether it is right, perhaps this is one of possible ways
1689 if (ener[0] < thr[qshell - 1]) {
1690 long ne;
1691 for (ne = 0; ne < ener.get_qel() && (ener[ne] < thr[qshell - 1] ||
1692 (ener[ne] >= thr[qshell - 1] &&
1693 ne > 1 && CS[ne - 1] <= CS[ne - 2]));
1694 ne++) {
1695 if (ne > 0) left_CS[ne - 1] = 0.0;
1696 nce = ne;
1697 }
1698 }
1699 //Iprintn(mcout, nce);
1700 //mcout<<"ener[nce]="<<ener[nce]<<" CS[nce]="<<CS[nce]
1701 // <<" CS[nce+1]="<<CS[nce+1]<<std::endl;
1702 int s_more;
1703 int nt2 = 0; // < nt1
1704 int s_spes = 0;
1705 do // Actually this is a loop by the group of thresholds
1706 {
1707 //mcout<<"nct="<<nct<<" nce="<<nce<<std::endl;
1708 // Find all thresholds which are less then the current energy
1709 int nt;
1710 s_more = 0; // sign that there are thresholds more than
1711 // the current energy
1712 for (nt = nct; nt >= 0; nt--) {
1713 if (s_spes == 0) {
1714 if (thr[nt] > ener[nce]) {
1715 s_more = 1;
1716 break;
1717 }
1718 } else {
1719 if (thr[nt] > ener[nce + 1]) {
1720 s_more = 1;
1721 break;
1722 }
1723 }
1724 }
1725 // nt is now index of the next threshold or -1, if the thresholds are
1726 // made up.
1727 int nt1 = nct;
1728 int nce_next = ener.get_qel();
1729 nt2 = nt + 1;
1730 //mcout<<"nt="<<nt<<" nt1="<<nt1<<" nt2="<<nt2<<" s_more="<<s_more<<'\n';
1731 if (s_more == 1) {
1732 //if(nt >= 0) // so if there are other larger thresholds,
1733 //{ // we should check how far we can pass at this step
1734 long ne;
1735 for (ne = nce; ne < ener.get_qel();
1736 ne++) { // finding energy larger than the next threshold
1737 if (thr[nt] <= ener[ne]) {
1738 nce_next = ne;
1739 s_spes = 0;
1740 break;
1741 } else { // At the following condition energy could be less then
1742 // threshold,
1743 // but this point will anyway be assotiated with the next shell
1744 // corresponding to this threshold.
1745 // This is related to not precise measurement of cross section
1746 // and not precise knowledge of shell energies.
1747 // Occurence of this condition is marked by s_spes = 1.
1748 // At the next passing of this loop the thresholds are compared with
1749 // the next energy.
1750 if (ne > 1 && ne < ener.get_qel() - 1 && ne > nce + 2 &&
1751 thr[nt] <= ener[ne + 1] &&
1752 (thr[nt] - ener[ne]) / (ener[ne + 1] - ener[ne]) < 0.1 &&
1753 CS[ne] > CS[ne - 1]) {
1754 //mcout<<"special condition is satisf.\n";
1755 nce_next = ne;
1756 s_spes = 1;
1757 break;
1758 }
1759 }
1760 }
1761 if (ne == ener.get_qel()) // threshold is larger then energy mesh
1762 s_more = 0; // to finish the loop
1763 }
1764 //Iprintn(mcout, nce_next);
1765 //Iprintn(mcout, ener[nce_next-1]);
1766 int qt = nt1 - nt2 + 1; // quantity of the new thresholds
1767 //Iprintn(mcout, qt);
1768 DynLinArr<double> w(qt); // weights accouring to charges
1769 int s = 0; // sum of Z
1770 for (nt = 0; nt < qt; nt++) {
1771 int nshell = nct - nt;
1772 s += Zshell[nshell];
1773 }
1774 for (nt = 0; nt < qt; nt++) {
1775 int nshell = nct - nt;
1776 w[nt] = double(Zshell[nshell]) / s;
1777 }
1778 double save_left_CS = left_CS[nce_next - 1];
1779 long ne;
1780 for (ne = 0; ne < nce_next; ne++) {
1781 for (nt = 0; nt < qt; nt++) {
1782 int nshell = nct - nt;
1783 SCS[nshell][ne] = left_CS[ne] * w[nt];
1784 }
1785 left_CS[ne] = 0.0;
1786 }
1787 for (ne = nce_next; ne < ener.get_qel(); ne++) {
1788 double extrap_CS =
1789 save_left_CS * pow(ener[nce_next - 1], 2.75) / pow(ener[ne], 2.75);
1790 if (extrap_CS > left_CS[ne]) extrap_CS = left_CS[ne];
1791 for (nt = 0; nt < qt; nt++) {
1792 int nshell = nct - nt;
1793 SCS[nshell][ne] = extrap_CS * w[nt];
1794 }
1795 left_CS[ne] -= extrap_CS;
1796 }
1797 nce = nce_next;
1798 nct = nt2 - 1;
1799 } while (s_more != 0);
1800 // now nt2 will be index of last filled shell
1801 // Now to fill the shells which are absent in the input table.
1802 // They will be initialized phenomenologically, basing on the sum rule:
1803 int ns;
1804 //Iprintn(mcout, nt2);
1805 for (ns = 0; ns < nt2; ns++) {
1806 acs[ns].pass(new PhenoPhotoAbsCS(shell_name[ns], Zshell[ns], thr[ns]));
1807 }
1808 // Initialization of input shells:
1809 for (ns = nt2; ns < qshell; ns++) {
1810 //acs[ns].pass( new SimpleTablePhotoAbsCS(shell_name[ns], Zshell[ns],
1811 // thr[ns],
1812 // ener, SCS[ns]) );
1814 shell_name[ns], Zshell[ns], thr[ns], ener, SCS[ns]);
1815 adr->remove_leading_zeros();
1816 acs[ns].pass(adr);
1817 }
1819 exener[0] = exener[1] = 0.0;
1820 double integ = get_integral_ACS(0.0, DBL_MAX);
1821 //Iprintn(mcout, integ);
1822 integ_abs_before_corr = integ;
1823 double pred_integ = Thomas_sum_rule_const_Mb * Z;
1824 //Iprintn(mcout, pred_integ);
1825 if (pred_integ > integ) {
1827 // add excitation
1828 exener[0] =
1829 low_boundary_of_excitations * acs[qshell - 1]->get_threshold();
1830 exener[1] = 1.0 * acs[qshell - 1]->get_threshold();
1831 height_of_excitation = (pred_integ - integ) / (exener[1] - exener[0]);
1832 if (minimal_threshold > 0.0) {
1833 if (minimal_threshold > acs[qshell - 1]->get_threshold()) {
1834 // currently the minimal shell is the last one
1835 exener[0] += minimal_threshold - acs[qshell - 1]->get_threshold();
1836 exener[1] += minimal_threshold - acs[qshell - 1]->get_threshold();
1837 }
1838 }
1839 }
1840 /*
1841 else
1842 {
1843 height_of_excitation = 0.0;
1844 exener[0] = exener[1] = 0.0;
1845 }
1846 */
1847 } else if (pred_integ < integ) {
1849 const double fact = pred_integ / integ;
1850 for (int nshell = 0; nshell < qshell; ++nshell) {
1851 acs[nshell]->scale(fact);
1852 }
1853 }
1854 }
1855
1856 integ_abs_after_corr = get_integral_ACS(0.0, DBL_MAX);
1858}
1859
1861 const String& fBT_file_name, int id,
1862 double fminimal_threshold)
1863 : threshold_file_name("none"),
1864 simple_table_file_name("none"),
1865 BT_file_name(fBT_file_name),
1866 minimal_threshold(fminimal_threshold) {
1867 mfunnamep("ExAtomPhotoAbsCS::ExAtomPhotoAbsCS(int fZ, const String& fname, "
1868 "const String& fBT_file_name, int id, double fminimal_threshold)");
1869 check_econd11(fZ, < 1, mcerr);
1870 check_econd21(id, < 1 ||, > 2, mcerr);
1871
1872 name = fname;
1873#ifdef USE_STLSTRING
1874 std::ifstream BT_file(BT_file_name.c_str());
1875#else
1876 std::ifstream BT_file(BT_file_name);
1877#endif
1878 if (!BT_file) {
1879 funnw.ehdr(mcerr);
1880 mcerr << "cannot open file " << BT_file_name << std::endl;
1881 spexit(mcerr);
1882 }
1883 DynLinArr<double> thresh(0, 0.0);
1884 DynLinArr<double> fl(0, 0.0);
1885 Z = fZ;
1886 int i = findmark(BT_file, "NUCLEAR CHARGE =");
1887 check_econd11a(i, != 1, "wrong file format", mcerr);
1888 int Z_from_file;
1889 BT_file >> Z_from_file;
1890 check_econd12(Z_from_file, !=, Z, mcerr);
1891 qshell = 0;
1892 while ((i = findmark(BT_file, "Z =")) == 1) {
1893 BT_file >> i;
1894 check_econd11(i, != Z, mcerr);
1895 String shellname;
1896 BT_file >> shellname;
1897 //Iprintn(mcout, shellname);
1898 i = findmark(BT_file, "$");
1899 check_econd11(i, != 1, mcerr);
1900 long qen;
1901 BT_file >> qen;
1902 check_econd11(qen, <= 0, mcerr);
1903 DynLinArr<double> fener(qen, 0.0);
1904 DynLinArr<double> fcs(qen, 0.0);
1905 double thr = 0.0;
1906 BT_file >> thr;
1907 check_econd11(thr, <= 0, mcerr);
1908 thr *= 1.0e-3; // pass from keV to MeV
1909 if (id == 2) {
1910 thresh.increment();
1911 fl.increment();
1912 thresh[qshell] = thr;
1913 BT_file >> fl[qshell];
1914 check_econd21(fl[qshell], < 0.0 ||, > 1.0, mcerr);
1915 //Iprintn(mcout, fl[qshell]);
1916 }
1917 long nen;
1918 for (nen = 0; nen < qen; nen++) {
1919 BT_file >> fener[nen] >> fcs[nen];
1920 check_econd11(fener[nen], <= 0.0, mcerr);
1921 check_econd11(fcs[nen], < 0.0, mcerr);
1922 fener[nen] *= 1.0e-3; // pass from keV to MeV
1923 }
1924 qshell++;
1926 acs[qshell - 1]
1927 .pass(new SimpleTablePhotoAbsCS(shellname, 0, // unknown here
1928 thr, fener, fcs));
1929 }
1930 if (id == 2) {
1931 // a copy of similar thing from subroutine above
1932 int n_min = 0;
1933 double st = DBL_MAX;
1934 for (int nshell = 0; nshell < qshell; ++nshell) {
1935 // currently the minimal shell is the last,
1936 // but to avoid this assumption we check all
1937 if (thresh[nshell] < st) {
1938 n_min = nshell;
1939 st = thresh[nshell];
1940 }
1941 }
1943 for (int nshell = 0; nshell < qshell; ++nshell) {
1944 if (fl[nshell] > 0) {
1945 check_econd12(nshell, ==, n_min, mcerr);
1946 DynLinArr<double> felectron_energy(0);
1947 DynLinArr<double> fphoton_energy(1);
1948 fphoton_energy[0] = thresh[nshell] - thresh[n_min];
1949 asp[nshell].add_channel(fl[nshell], felectron_energy, fphoton_energy);
1950 }
1951 }
1952 }
1953
1954 check_econd11(qshell, <= 0, mcerr);
1957 exener[0] = exener[1] = 0.0;
1958 double integ = get_integral_ACS(0.0, DBL_MAX);
1959 //Iprintn(mcout, integ);
1960 integ_abs_before_corr = integ;
1961 double pred_integ = Thomas_sum_rule_const_Mb * Z;
1962 //Iprintn(mcout, pred_integ);
1963 if (pred_integ > integ) {
1965 // add excitation
1966 exener[0] =
1967 low_boundary_of_excitations * acs[qshell - 1]->get_threshold();
1968 exener[1] = 1.0 * acs[qshell - 1]->get_threshold();
1969 height_of_excitation = (pred_integ - integ) / (exener[1] - exener[0]);
1970 if (minimal_threshold > 0.0) {
1971 if (minimal_threshold > acs[qshell - 1]->get_threshold()) {
1972 // currently the minimal shell is the last one
1973 exener[0] += minimal_threshold - acs[qshell - 1]->get_threshold();
1974 exener[1] += minimal_threshold - acs[qshell - 1]->get_threshold();
1975 }
1976 }
1977 }
1978 } else {
1980 const double fact = pred_integ / integ;
1981 for (int nshell = 0; nshell < qshell; ++nshell) {
1982 acs[nshell]->scale(fact);
1983 }
1984 }
1985 }
1986 integ_abs_after_corr = get_integral_ACS(0.0, DBL_MAX);
1988}
1989
1990#define READ_FILE_WITH_PRINCIPAL_NUMBERS
1991
1993 const String& fFitBT_file_name,
1994 int id, // to distinguish it from
1995 // constructor above
1996 int s_no_scale, double fminimal_threshold)
1997 : threshold_file_name("none"),
1998 simple_table_file_name("none"),
1999 BT_file_name(fFitBT_file_name),
2000 minimal_threshold(fminimal_threshold) {
2001 mfunnamep(
2002 "ExAtomPhotoAbsCS::ExAtomPhotoAbsCS(int fZ, const String& fname, const "
2003 "String& fFitBT_file_name, int id, int id1, double fminimal_threshold)");
2004 check_econd11(fZ, < 1, mcerr);
2005 check_econd21(id, < 1 ||, > 2, mcerr);
2006 Z = fZ;
2007 name = fname;
2008#ifdef USE_STLSTRING
2009 std::ifstream BT_file(fFitBT_file_name.c_str());
2010#else
2011 std::ifstream BT_file(fFitBT_file_name);
2012#endif
2013 if (!BT_file) {
2014 funnw.ehdr(mcerr);
2015 mcerr << "cannot open file " << BT_file_name << std::endl;
2016 spexit(mcerr);
2017 }
2018 DynLinArr<double> thresh(0, 0.0);
2019 DynLinArr<double> fl(0, 0.0);
2020
2021 int i = 0;
2022 while ((i = findmark(BT_file, "$")) == 1) {
2023 long iZ;
2024 BT_file >> iZ;
2025 //Iprintn(mcout, iZ);
2026 if (iZ == Z) {
2027 BT_file >> qshell;
2028 //Iprintn(mcout, qshell);
2029 check_econd11(qshell, <= 0, mcerr);
2030 check_econd11(qshell, > 1000, mcerr);
2032 if (id == 2) {
2033 thresh.put_qel(qshell);
2034 fl.put_qel(qshell);
2035 }
2036 for (int nshell = 0; nshell < qshell; ++nshell) {
2037#ifdef READ_FILE_WITH_PRINCIPAL_NUMBERS
2038 int n_princ = 0;
2039#endif
2040 int l;
2041 double threshold;
2042 double E0;
2043 double yw;
2044 double ya;
2045 double P;
2046 double sigma;
2047 if (BT_file.eof()) {
2048 mcerr << "unexpected end of file " << BT_file_name << '\n';
2049 spexit(mcerr);
2050 }
2051 if (!BT_file.good()) {
2052 mcerr << "bad format of file " << BT_file_name << '\n';
2053 spexit(mcerr);
2054 }
2055#ifdef READ_FILE_WITH_PRINCIPAL_NUMBERS
2056 BT_file >> n_princ;
2057 if (!BT_file.good()) {
2058 mcerr << "bad format of file " << BT_file_name << '\n';
2059 spexit(mcerr);
2060 }
2061 check_econd21(n_princ, < 0 ||, > 10, mcerr);
2062#endif
2063 BT_file >> l >> threshold >> E0 >> sigma >> ya >> P >> yw;
2064 check_econd11(l, < 0, mcerr);
2065 check_econd11(l, > 20, mcerr);
2066 threshold *= 1.0e-6;
2067 E0 *= 1.0e-6;
2068
2069 check_econd11a(threshold, <= 2.0e-6,
2070 "n_princ=" << n_princ << " l=" << l << '\n', mcerr);
2071 check_econd11(E0, <= 0, mcerr);
2072 double flu = 0.0;
2073 if (id == 2) {
2074 if (BT_file.eof()) {
2075 mcerr << "unexpected end of file " << BT_file_name << '\n';
2076 spexit(mcerr);
2077 }
2078 if (!BT_file.good()) {
2079 mcerr << "bad format of file " << BT_file_name << '\n';
2080 spexit(mcerr);
2081 }
2082 BT_file >> flu;
2083 check_econd11(flu, < 0.0, mcerr);
2084 check_econd11(flu, > 1.0, mcerr);
2085 thresh[nshell] = threshold;
2086 fl[nshell] = flu;
2087 }
2088#ifdef READ_FILE_WITH_PRINCIPAL_NUMBERS
2089 String shellname(long_to_String(n_princ) + // necessary
2090 // for generation escape products
2091 String(" shell number ") + long_to_String(nshell));
2092#else
2093 String shellname(String("shell number ") + long_to_String(nshell));
2094#endif
2095 acs[nshell].pass(
2096 new SimpleTablePhotoAbsCS(shellname, 0, // unknown here
2097 threshold, l, E0, yw, ya, P, sigma));
2098 //Iprintn(mcout, nshell);
2099 //Iprint3n(mcout, l, threshold, E0);
2100 //Iprint4n(mcout, yw, ya, P, sigma);
2101 //acs[nshell]->print(mcout, 5);
2102 }
2103 goto mark1;
2104 }
2105 }
2106 funnw.ehdr(mcerr);
2107 mcerr << "there is no element Z=" << fZ << " in file " << fFitBT_file_name
2108 << std::endl;
2109 spexit(mcerr);
2110mark1:
2111 if (id == 2) {
2112 // a copy of similar thing from subroutine above
2113 int n_min = 0;
2114 double st = DBL_MAX;
2115 for (int nshell = 0; nshell < qshell; ++nshell) {
2116 // currently the minimal shell is the last,
2117 // but to avoid this assumption we check all
2118 if (thresh[nshell] < st) {
2119 n_min = nshell;
2120 st = thresh[nshell];
2121 }
2122 }
2124 for (int nshell = 0; nshell < qshell; ++nshell) {
2125 if (fl[nshell] > 0) {
2126 check_econd12(nshell, ==, n_min, mcerr);
2127 DynLinArr<double> felectron_energy(0);
2128 DynLinArr<double> fphoton_energy(1);
2129 fphoton_energy[0] = thresh[nshell] - thresh[n_min];
2130 asp[nshell].add_channel(fl[nshell], felectron_energy, fphoton_energy);
2131 }
2132 }
2133 }
2134
2135 check_econd11(qshell, <= 0, mcerr);
2138 exener[0] = exener[1] = 0.0;
2139 double integ = get_integral_ACS(0.0, DBL_MAX);
2140 //Iprintn(mcout, integ);
2141 integ_abs_before_corr = integ;
2142 double pred_integ = Thomas_sum_rule_const_Mb * Z;
2143 //Iprintn(mcout, pred_integ);
2144 if (pred_integ > integ) {
2146 // add excitation
2147 exener[0] =
2148 low_boundary_of_excitations * acs[qshell - 1]->get_threshold();
2149 exener[1] = 1.0 * acs[qshell - 1]->get_threshold();
2150 height_of_excitation = (pred_integ - integ) / (exener[1] - exener[0]);
2151 if (minimal_threshold > 0.0) {
2152 if (minimal_threshold > acs[qshell - 1]->get_threshold()) {
2153 // currently the minimal shell is the last one
2154 exener[0] += minimal_threshold - acs[qshell - 1]->get_threshold();
2155 exener[1] += minimal_threshold - acs[qshell - 1]->get_threshold();
2156 }
2157 }
2158 }
2159 } else {
2160 if (s_scale_to_normalize_if_more == 1 && s_no_scale == 0) {
2161 const double fact = pred_integ / integ;
2162 for (int nshell = 0; nshell < qshell; ++nshell) {
2163 acs[nshell]->scale(fact);
2164 }
2165 }
2166 }
2167 integ_abs_after_corr = get_integral_ACS(0.0, DBL_MAX);
2169}
2170
2172 int fZ, const String& fname, // just name of this atom
2173 const String& fFitBT_file_name,
2174 // const String& fthreshold_file_name, // minimal potential corresponding
2175 // to:
2176 const String& fsimple_table_file_name,
2177 double emax_repl, int id, // to distinguish it from constructor above
2178 double fminimal_threshold) {
2179 mfunname("ExAtomPhotoAbsCS::ExAtomPhotoAbsCS(...)");
2180 Z = fZ;
2181 name = fname;
2182 int s_no_scale = 1;
2183 *this =
2184 ExAtomPhotoAbsCS(fZ, fname, // just name of this atom
2185 fFitBT_file_name, id, s_no_scale, fminimal_threshold);
2186
2188 exener[0] = exener[1] = 0.0;
2189
2190 double thrmin = DBL_MAX;
2191 long nsmin = -1;
2192 // look for minimal shell (usually the last)
2193 for (long ns = 0; ns < qshell; ++ns) {
2194 //Iprint2n(mcout, ns, acs[ns]->get_threshold());
2195 if (thrmin > acs[ns]->get_threshold()) {
2196 nsmin = ns;
2197 thrmin = acs[ns]->get_threshold();
2198 }
2199 }
2200 //Iprint3n(mcout, nsmin, acs[nsmin]->get_threshold(), thrmin);
2201 check_econd11(nsmin, < 0, mcerr);
2202 check_econd11(nsmin, != qshell - 1, mcerr); // now it has to be by this way
2203 ActivePtr<PhotoAbsCS> facs = acs[nsmin]; // copying the valence shell
2204 PhotoAbsCS* apacs = facs.get();
2205 //PhotoAbsCS* apacs = acs[nsmin]->getver();
2206 SimpleTablePhotoAbsCS* first_shell =
2207 dynamic_cast<SimpleTablePhotoAbsCS*>(apacs);
2208 //mcout<<"first_shell:\n";
2209 //first_shell->print(mcout, 2);
2210 // Strange why the following does not work (because of something being not
2211 // implemented in commpiler)? - may be due to an error before, which is now
2212 // corrected.
2213 // assumes the first is the last
2214 // SimpleTablePhotoAbsCS* first_shell =
2215 // dynamic_cast<SimpleTablePhotoAbsCS*>(&(acs[nsmin]->getver()));
2216
2217 check_econd11(first_shell, == NULL, mcerr);
2218 /*
2219#ifdef USE_STLSTRING
2220 std::ifstream threshold_file(threshold_file_name.c_str());
2221#else
2222 std::ifstream threshold_file(threshold_file_name);
2223#endif
2224 if( !threshold_file )
2225 {
2226 funnw.ehdr(mcerr);
2227 mcerr<<"cannot open file "<<threshold_file_name<<std::endl;
2228 spexit(mcerr);
2229 }
2230 DynLinArr< double > thr(0, 0.0);
2231 DynLinArr< int > Zshell(0, 0);
2232 DynLinArr< double > fl(0, 0.0);
2233 DynLinArr< String > shell_name(0);
2234 int n_min=0;
2235 double thr_min = DBL_MAX;
2236 while (findmark(threshold_file, "#") == 1) {
2237 threshold_file >> Z;
2238 if (Z == fZ) {
2239 threshold_file >> qshell;
2240 check_econd21( qshell , < 1 || , > 10000 , mcerr);
2241 s_ignore_shell.put_qel(qshell, 0);
2242 //Iprintn(mcout, qshell);
2243 thr = DynLinArr< double >(qshell, 0.0);
2244 Zshell = DynLinArr< int >(qshell, 0);
2245 fl = DynLinArr< double >(qshell, 0.0);
2246 shell_name = DynLinArr< String >(qshell);
2247 String temp_name;
2248 threshold_file >> temp_name;
2249 if(fname == "none") name = temp_name;
2250 else name = fname;
2251 int sZshell = 0;
2252 for (int nshell = 0; nshell < qshell; nshell++) {
2253 threshold_file >> thr[nshell];
2254 check_econd11(thr[nshell] , <= 0.0, mcerr);
2255 thr[nshell] *= 1.0e-6;
2256 threshold_file >> Zshell[nshell];
2257 check_econd11(Zshell[nshell] , <= 0, mcerr);
2258 sZshell += Zshell[nshell];
2259 threshold_file >> fl[nshell];
2260 findmark(threshold_file, "!");
2261 threshold_file >> shell_name[nshell];
2262 //PhotoAbsCS* ap = new PhenoPhotoAbsCS(shell_name, Zshell, thr);
2263 //acs[nshell].pass( ap );
2264 }
2265 check_econd12( sZshell , != , Z , mcerr);
2266 for(nshell=0; nshell<qshell; nshell++) {
2267 // currently the minimal shell is the last,
2268 // but to avoid this assumption we check all
2269 if (thr[nshell] < thr_min) {
2270 n_min = nshell;
2271 thr_min = thr[nshell];
2272 }
2273 }
2274 goto mark1;
2275 }
2276 }
2277 funnw.ehdr(mcerr);
2278 mcerr<<"there is no element Z="<<fZ<<" in file
2279"<<threshold_file_name<<std::endl;
2280 spexit(mcerr);
2281 mark1:
2282 */
2283
2284 SimpleTablePhotoAbsCS stpacs(name, Z, 0.0, fsimple_table_file_name);
2285 stpacs.remove_leading_tiny(1.0e-10);
2286
2287 // Merging shells:
2288 acs[nsmin].pass(new SimpleTablePhotoAbsCS(*first_shell, stpacs, emax_repl));
2289
2292 exener[0] = exener[1] = 0.0;
2293 double integ = get_integral_ACS(0.0, DBL_MAX);
2294 //Iprintn(mcout, integ);
2295 integ_abs_before_corr = integ;
2296 double pred_integ = Thomas_sum_rule_const_Mb * Z;
2297 //Iprintn(mcout, pred_integ);
2298 if (pred_integ > integ) {
2300 // add excitation
2301 exener[0] =
2302 low_boundary_of_excitations * acs[qshell - 1]->get_threshold();
2303 exener[1] = 1.0 * acs[qshell - 1]->get_threshold();
2304 height_of_excitation = (pred_integ - integ) / (exener[1] - exener[0]);
2305 if (minimal_threshold > 0.0) {
2306 if (minimal_threshold > acs[qshell - 1]->get_threshold()) {
2307 // currently the minimal shell is the last one
2308 exener[0] += minimal_threshold - acs[qshell - 1]->get_threshold();
2309 exener[1] += minimal_threshold - acs[qshell - 1]->get_threshold();
2310 }
2311 }
2312 }
2313 } else {
2315 const double fact = pred_integ / integ;
2316 for (int nshell = 0; nshell < qshell; ++nshell) {
2317 acs[nshell]->scale(fact);
2318 }
2319 }
2320 }
2321 integ_abs_after_corr = get_integral_ACS(0.0, DBL_MAX);
2323}
2324
2325double ExAtomPhotoAbsCS::get_threshold(int nshell) const {
2326 mfunname("double ExAtomPhotoAbsCS::get_threshold(int nshell) const");
2327 check_econd21(nshell, < 0 ||, > qshell, mcerr);
2328 double r = acs[nshell]->get_threshold();
2329 if (minimal_threshold > 0.0) {
2331 }
2332 return r;
2333}
2334
2335double ExAtomPhotoAbsCS::get_ICS(double energy) const {
2336 mfunname("double ExAtomPhotoAbsCS::get_ACS(double energy) const");
2337 double s = 0.0;
2338 for (int n = 0; n < qshell; ++n) {
2339 if (s_ignore_shell[n] == 0) {
2340 double shift = 0.0;
2341 const double t = acs[n]->get_threshold();
2342 if (minimal_threshold > 0.0) {
2343 if (t < minimal_threshold) shift = minimal_threshold - t;
2344 }
2345 //Iprint3n(mcout, n, t, shift);
2346 s += acs[n]->get_CS(energy - shift);
2347 //Iprintn(mcout, s);
2348 }
2349 }
2350 return s;
2351}
2352
2354 double energy2) const {
2355 mfunname("double ExAtomPhotoAbsCS::get_integral_ICS(double energy)const");
2356 double s = 0.0;
2357 for (int n = 0; n < qshell; ++n) {
2358 if (s_ignore_shell[n] == 0) {
2359 double shift = 0.0;
2360 const double t = acs[n]->get_threshold();
2361 if (minimal_threshold > 0.0) {
2362 if (t < minimal_threshold) shift = minimal_threshold - t;
2363 }
2364 s += acs[n]->get_integral_CS(energy1 - shift, energy2 - shift);
2365 }
2366 }
2367 return s;
2368}
2369
2370double ExAtomPhotoAbsCS::get_ICS(int nshell, double energy) const {
2371 mfunname("double ExAtomPhotoAbsCS::get_ICS(int nshell, double energy)");
2372 check_econd21(nshell, < 0 ||, > qshell, mcerr);
2373 if (s_ignore_shell[nshell] == 0) {
2374 double shift = 0.0;
2375 const double t = acs[nshell]->get_threshold();
2376 if (minimal_threshold > 0.0) {
2377 if (t < minimal_threshold) shift = minimal_threshold - t;
2378 }
2379 return acs[nshell]->get_CS(energy - shift);
2380 }
2381 return 0.0;
2382}
2383
2384double ExAtomPhotoAbsCS::get_integral_ICS(int nshell, double energy1,
2385 double energy2) const {
2386 mfunname("double ExAtomPhotoAbsCS::get_integral_ICS(int nshell, double "
2387 "energy1, double energy2)");
2388 check_econd21(nshell, < 0 ||, > qshell, mcerr);
2389 if (s_ignore_shell[nshell] == 0) {
2390 double shift = 0.0;
2391 const double t = acs[nshell]->get_threshold();
2392 if (minimal_threshold > 0.0) {
2393 if (t < minimal_threshold) shift = minimal_threshold - t;
2394 }
2395 return acs[nshell]->get_integral_CS(energy1 - shift, energy2 - shift);
2396 }
2397 return 0.0;
2398}
2399
2400double ExAtomPhotoAbsCS::get_ACS(double energy) const {
2401 mfunname("double ExAtomPhotoAbsCS::get_ACS(double energy) const");
2402 double s = 0.0;
2403 for (int n = 0; n < qshell; ++n) {
2404 if (s_ignore_shell[n] == 0) {
2405 double shift = 0.0;
2406 const double t = acs[n]->get_threshold();
2407 if (minimal_threshold > 0.0) {
2408 if (t < minimal_threshold) shift = minimal_threshold - t;
2409 }
2410 s += acs[n]->get_CS(energy - shift);
2411 }
2412 }
2413 if (energy >= exener[0] && energy <= exener[1]) {
2415 }
2416 return s;
2417}
2418
2420 double energy2) const {
2421 mfunname("double ExAtomPhotoAbsCS::get_integral_ACS(double energy1, double "
2422 "energy2) const");
2423 double s = 0.0;
2424 for (int n = 0; n < qshell; ++n) {
2425 if (s_ignore_shell[n] == 0) {
2426 double shift = 0.0;
2427 const double t = acs[n]->get_threshold();
2428 if (minimal_threshold > 0.0) {
2429 if (t < minimal_threshold) shift = minimal_threshold - t;
2430 }
2431 s += acs[n]->get_integral_CS(energy1 - shift, energy2 - shift);
2432 }
2433 }
2434 //Imcout<<"energy1="<<setw(13)<<energy1
2435 // <<"energy2="<<setw(13)<<energy2
2436 // <<" s="<<s<<'\n';
2437 double b[2];
2438 b[0] = find_max(exener[0], energy1);
2439 b[1] = find_min(exener[1], energy2);
2440 //Iprint3n(mcout, b[0], b[1], height_of_excitation);
2441 if (b[1] >= b[0]) {
2442 s += height_of_excitation * (b[1] - b[0]);
2443 }
2444 //Iprintn(mcout, s);
2445 return s;
2446}
2447
2448double ExAtomPhotoAbsCS::get_ACS(int nshell, double energy) const {
2449 mfunname("double ExAtomPhotoAbsCS::get_ACS(int nshell, double energy)");
2450 check_econd21(nshell, < 0 ||, > qshell, mcerr);
2451 if (s_ignore_shell[nshell] == 0) {
2452 double shift = 0.0;
2453 const double t = acs[nshell]->get_threshold();
2454 if (minimal_threshold > 0.0) {
2455 if (t < minimal_threshold) shift = minimal_threshold - t;
2456 }
2457 double s = acs[nshell]->get_CS(energy - shift);
2458 if (nshell == qshell - 1 && energy >= exener[0] && energy <= exener[1]) {
2460 }
2461 return s;
2462 }
2463 return 0.0;
2464}
2465
2466double ExAtomPhotoAbsCS::get_integral_ACS(int nshell, double energy1,
2467 double energy2) const {
2468 mfunname("double ExAtomPhotoAbsCS::get_integral_ACS(int nshell, double "
2469 "energy1, double energy2)");
2470 check_econd21(nshell, < 0 ||, > qshell, mcerr);
2471 if (s_ignore_shell[nshell] == 0) {
2472 double shift = 0.0;
2473 const double t = acs[nshell]->get_threshold();
2474 if (minimal_threshold > 0.0) {
2475 if (t < minimal_threshold) shift = minimal_threshold - t;
2476 }
2477 double s = acs[nshell]->get_integral_CS(energy1 - shift, energy2 - shift);
2478 //Imcout<<"energy1="<<setw(13)<<energy1<<" s="<<s<<'\n';
2479 if (nshell == qshell - 1) {
2480 double b[2];
2481 b[0] = find_max(exener[0], energy1);
2482 b[1] = find_min(exener[1], energy2);
2483 if (b[1] >= b[0]) {
2484 s += height_of_excitation * (b[1] - b[0]);
2485 }
2486 }
2487 return s;
2488 }
2489 return 0.0;
2490}
2491
2492void ExAtomPhotoAbsCS::print(std::ostream& file, int l) const {
2493 if (l > 0) {
2494 Ifile << "ExAtomPhotoAbsCS(l=" << l << "): name=" << name << " Z = " << Z
2495 << " qshell = " << qshell << std::endl;
2496 indn.n += 2;
2497 Ifile << "threshold_file_name=" << threshold_file_name << '\n';
2498 Ifile << "simple_table_file_name=" << simple_table_file_name << '\n';
2499 Ifile << "BT_file_name=" << BT_file_name << std::endl;
2500 Ifile << "Thomas_sum_rule_const_Mb * Z = " << Thomas_sum_rule_const_Mb* Z
2501 << '\n';
2502 Ifile << "integ_abs_before_corr = " << integ_abs_before_corr << '\n';
2503 Ifile << "integ_abs_after_corr = " << integ_abs_after_corr << '\n';
2504 Ifile << "integ_ioniz_after_corr = " << integ_ioniz_after_corr
2505 << '\n';
2506 Ifile << "height_of_excitation=" << height_of_excitation
2507 << " exener=" << exener[0] << ' ' << exener[1] << '\n';
2509 Ifile << "integrals by shells:\n";
2510 Ifile << "nshell, int(acs), int(ics)\n";
2511 for (long n = 0; n < qshell; n++) {
2512 double ainteg = get_integral_ACS(n, 0.0, DBL_MAX);
2513 double iinteg = get_integral_ICS(n, 0.0, DBL_MAX);
2514 Ifile << n << " " << ainteg << " " << iinteg << '\n';
2515 }
2516
2517 if (l > 1) {
2518 l--;
2519 indn.n += 2;
2520 for (long n = 0; n < qshell; ++n) {
2521 Ifile << "nshell=" << n << std::endl;
2522 acs[n].print(file, l);
2523 }
2524 AtomPhotoAbsCS::print(file, l);
2525 indn.n -= 2;
2526 }
2527 indn.n -= 2;
2528 }
2529}
2530
2532 mfunname("int ExAtomPhotoAbsCS::get_main_shell_number(int nshell) const");
2533 String shell_name = acs[nshell]->get_name();
2534#ifdef STRSTREAM_AVAILABLE
2535#ifdef USE_STLSTRING
2536 istrstream sfile(shell_name.c_str());
2537#else
2538 istrstream sfile(shell_name);
2539#endif
2540#else
2541 std::istringstream sfile(shell_name.c_str());
2542#endif
2543 int i = -1;
2544 sfile >> i;
2545 //Iprintn(mcout, i);
2546 //check_econd(i < 1 || i > 50 , "i="<<i<<" shell_name="<<shell_name<<'\n' ,
2547 // mcerr);
2548 if (i < 1 || i > 50) {
2549 i = -1;
2550 }
2551 return i;
2552}
2553
2555 double fstep,
2556 long fmax_q_step) {
2557 mfunname("void ExAtomPhotoAbsCS::replace_shells_by_overage(...)");
2558 for (long n = 0; n < qshell; n++) {
2559 //Iprintn(mcout, n);
2560 //mcout<<"----------------before replacement:\n";
2561 //acs[n]->print(mcout, 10);
2562 PhotoAbsCS* a =
2563 new OveragePhotoAbsCS(acs[n].getver(), fwidth, fstep, fmax_q_step);
2564 //mcout<<"----------------after replacement:\n";
2565 //acs[n]->print(mcout, 10);
2566 acs[n].pass(a);
2567 }
2568}
2569
2570//---------------------------------------------------------
2571
2573 double fW, double fF)
2574 : W(fW), F(fF) {
2575 qatom = fqatom;
2576 qatom_ps.put_qel(1, qatom);
2577 atom.put_qel(1, PassivePtr<const AtomPhotoAbsCS>(&fatom));
2578 if (W == 0.0) {
2579 W = coef_I_to_W * atom[0]->get_I_min();
2580 }
2581}
2582
2584 const AtomPhotoAbsCS& fatom2, int fqatom_ps2,
2585 double fW, double fF)
2586 : W(fW), F(fF) {
2587 qatom = fqatom_ps1 + fqatom_ps2;
2588 qatom_ps.put_qel(2);
2589 atom.put_qel(2);
2590 qatom_ps[0] = fqatom_ps1;
2591 qatom_ps[1] = fqatom_ps2;
2592 atom[0] = fatom1;
2593 atom[1] = fatom2;
2594 if (W == 0.0) {
2595#ifdef CALC_W_USING_CHARGES
2596 W = coef_I_to_W * (qatom_ps[0] * atom[0]->get_Z() * atom[0]->get_I_min() +
2597 qatom_ps[1] * atom[1]->get_Z() * atom[1]->get_I_min()) /
2598 (qatom_ps[0] * atom[0]->get_Z() + qatom_ps[1] * atom[1]->get_Z());
2599#else
2600 W = coef_I_to_W * (qatom_ps[0] * atom[0]->get_I_min() +
2601 qatom_ps[1] * atom[1]->get_I_min()) / qatom;
2602#endif
2603 }
2604}
2605
2607 const AtomPhotoAbsCS& fatom2, int fqatom_ps2,
2608 const AtomPhotoAbsCS& fatom3, int fqatom_ps3,
2609 double fW, double fF)
2610 : W(fW), F(fF) {
2611 qatom = fqatom_ps1 + fqatom_ps2 + fqatom_ps3;
2612 qatom_ps.put_qel(3);
2613 atom.put_qel(3);
2614 qatom_ps[0] = fqatom_ps1;
2615 qatom_ps[1] = fqatom_ps2;
2616 qatom_ps[2] = fqatom_ps3;
2617 atom[0] = fatom1;
2618 atom[1] = fatom2;
2619 atom[2] = fatom3;
2620 if (W == 0.0) {
2621#ifdef CALC_W_USING_CHARGES
2622 W = coef_I_to_W * (qatom_ps[0] * atom[0]->get_Z() * atom[0]->get_I_min() +
2623 qatom_ps[1] * atom[1]->get_Z() * atom[1]->get_I_min() +
2624 qatom_ps[2] * atom[2]->get_Z() * atom[2]->get_I_min()) /
2625 (qatom_ps[0] * atom[0]->get_Z() + qatom_ps[1] * atom[1]->get_Z() +
2626 qatom_ps[2] * atom[2]->get_Z());
2627#else
2628 W = coef_I_to_W * (qatom_ps[0] * atom[0]->get_I_min() +
2629 qatom_ps[1] * atom[1]->get_I_min() +
2630 qatom_ps[2] * atom[2]->get_I_min()) / qatom;
2631#endif
2632 }
2633}
2634
2635double MolecPhotoAbsCS::get_ACS(double energy) const {
2636 mfunname("double MolecPhotoAbsCS::get_ACS(double energy) const");
2637 const long q = qatom_ps.get_qel();
2638 double s = 0.0;
2639 for (long n = 0; n < q; n++) {
2640 s += qatom_ps[n] * atom[n]->get_ACS(energy);
2641 }
2642 return s;
2643}
2644
2645double MolecPhotoAbsCS::get_integral_ACS(double energy1, double energy2) const {
2646 mfunname("double MolecPhotoAbsCS::get_integral_ACS(double energy1, double "
2647 "energy2) const");
2648 const long q = qatom_ps.get_qel();
2649 double s = 0.0;
2650 for (long n = 0; n < q; n++) {
2651 s += qatom_ps[n] * atom[n]->get_integral_ACS(energy1, energy2);
2652 }
2653 return s;
2654}
2655
2656double MolecPhotoAbsCS::get_ICS(double energy) const {
2657 mfunname("double MolecPhotoAbsCS::get_ICS(double energy) const");
2658 const long q = qatom_ps.get_qel();
2659 double s = 0.0;
2660 for (long n = 0; n < q; n++) {
2661 s += qatom_ps[n] * atom[n]->get_ICS(energy);
2662 }
2663 return s;
2664}
2665
2666double MolecPhotoAbsCS::get_integral_ICS(double energy1, double energy2) const {
2667 mfunname("double MolecPhotoAbsCS::get_integral_ICS(double energy1, double "
2668 "energy2) const");
2669 const long q = qatom_ps.get_qel();
2670 double s = 0.0;
2671 for (long n = 0; n < q; n++) {
2672 s += qatom_ps[n] * atom[n]->get_integral_ICS(energy1, energy2);
2673 }
2674 return s;
2675}
2676
2678 mfunname("int MolecPhotoAbsCS::get_total_Z() const");
2679 int s = 0;
2680 const long q = qatom_ps.get_qel();
2681 for (long n = 0; n < q; n++) {
2682 s += atom[n]->get_Z();
2683 }
2684 return s;
2685}
2686
2687void MolecPhotoAbsCS::print(std::ostream& file, int l) const {
2688 Ifile << "MolecPhotoAbsCS (l=" << l << "):\n";
2689 Iprintn(file, qatom);
2690 Iprintn(file, W);
2691 Iprintn(file, F);
2692 const long q = qatom_ps.get_qel();
2693 Ifile << "number of sorts of atoms is " << q << '\n';
2694 indn.n += 2;
2695 for (long n = 0; n < q; n++) {
2696 Ifile << "n=" << n << " qatom_ps[n]=" << qatom_ps[n] << " atom:\n";
2697 atom[n]->print(file, l);
2698 }
2699 indn.n -= 2;
2700}
2701
2702std::ostream& operator<<(std::ostream& file, const MolecPhotoAbsCS& f) {
2703 f.print(file, 1);
2704 return file;
2705}
2706
2707/*
2708PhotoAbsorptionCS::PhotoAbsorptionCS(void):sqh(0) {};
2709
2710PhotoAbsorptionCS::PhotoAbsorptionCS(const String& ffile_name, PassivePtr<
2711EnergyMesh > fenergy_mesh) {
2712 mfunname("PhotoAbsorptionCS::PhotoAbsorptionCS(const String& ffile_name,
2713PassivePtr< EnergyMesh > fenergy_mesh)");
2714
2715 std::ifstream file(ffile_name);
2716 if (!file) {
2717 mcerr<<" can not open file "<<ffile_name<<'\n';
2718 spexit(mcerr);
2719 }
2720 file>>z;
2721 check_econd12(z , < 1 || , >= max_poss_atom_z , mcerr);
2722 file>>qsh;
2723 check_econd12(qsh , < 1 || , >= z , mcerr);
2724 String atname;
2725 file >> atname;
2726 shell_energy = DinLinArr< double >(qsh, 0.0);
2727 fluorescence_yield = DinLinArr< double >(qsh, 0.0);
2728 z_shell = DinLinArr< int >(qsh, 0);
2729
2730 for (int nsh = qsh-1; nsh >= 0; nsh--) {
2731 file >> shell_energy[nsh] >> z_shell[nsh] >> fluorescence_yield[nsh];
2732 check_econd11(shell_energy[nsh] , > 0.0 , mcerr);
2733 check_econd12(z_shell[nsh] , < 1 || , >= z , mcerr);
2734 check_econd12(fluorescence_yield[nsh] , < 0.0 || , > 1.0 , mcerr);
2735 int ci; // now pass comments
2736 while( !((ci=getc(fli))=='\n') ) {
2737 check_econd11( ci , == EOF , mcerr);
2738 }
2739 }
2740 long qen;
2741 file >> qen;
2742 check_econd11(qen , < 1 , mcerr);
2743 DinLinArr< double >en(qen, 0.0);
2744 total_cs = DinLinArr< double >(qen, 0.0);
2745 shell_cs = DinArr< double >(qsh, qen, 0.0);
2746 for (long nen = 0; nen<qen; ++nen) {
2747 file>>en[nen];
2748 if(nen > 0) {
2749 check_econd11(en[nen] , < en[nen-1] , mcerr);
2750 }
2751 file>>total_cs[nen];
2752 check_econd11(total_cs[nen] , < 0.0 , mcerr);
2753 double s = 0;
2754 for (nsh = 0; nsh < qsh; nsh++) {
2755 file>>shell_cs[nsh][nen];
2756 check_econd11(shell_cs[nen] , < 0.0 , mcerr);
2757 s+=shell_cs[nsh][nen];
2758 }
2759 s = fabs((total_cs[nen] - s)/(total_cs[nen] + s));
2760 check_econd11(s , > 0.001 , mcerr); // precision of printing
2761 }
2762 *fenergy_mesh = EnergyMesh(en);
2763}
2764*/
2765
2766}
@ do_clone
Definition: AbsPtr.h:516
DoubleAc pow(const DoubleAc &f, double p)
Definition: DoubleAc.cpp:336
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:313
DoubleAc find_min(const DoubleAc &a, const DoubleAc &b)
Definition: DoubleAc.h:627
DoubleAc find_max(const DoubleAc &a, const DoubleAc &b)
Definition: DoubleAc.h:655
#define check_econd21(a, sign1_b1_sign0, sign2_b2, stream)
Definition: FunNameStack.h:428
#define check_econd11(a, signb, stream)
Definition: FunNameStack.h:366
#define check_econd11a(a, signb, add, stream)
Definition: FunNameStack.h:395
#define mfunnamep(string)
Definition: FunNameStack.h:77
#define spexit(stream)
Definition: FunNameStack.h:536
#define check_econd12(a, sign, b, stream)
Definition: FunNameStack.h:380
#define mfunname(string)
Definition: FunNameStack.h:67
String long_to_String(long n)
Definition: String.h:102
std::string String
Definition: String.h:75
Definition: BlkArr.h:80
void put_qel(long fqel)
Definition: BlkArr.h:354
void increment(const T *val=NULL)
Definition: AbsArr.h:440
long get_qel(void) const
Definition: AbsArr.h:420
void put_qel(long fqel)
Definition: AbsArr.h:774
void pass(long fqel, T *fel)
Definition: AbsArr.h:265
virtual void remove_shell(int nshell)
virtual double get_ACS(double energy) const =0
virtual int get_main_shell_number(int nshell) const =0
int get_qshell() const
Definition: PhotoAbsCS.h:325
virtual void print(std::ostream &file, int l) const
virtual double get_TICS(double energy, double factual_minimal_threshold) const
Definition: PhotoAbsCS.cpp:992
virtual double get_threshold(int nshell) const =0
DynLinArr< int > s_ignore_shell
Definition: PhotoAbsCS.h:388
virtual double get_integral_TICS(double energy1, double energy2, double factual_minimal_threshold) const
virtual void get_escape_particles(int nshell, double energy, DynLinArr< double > &el_energy, DynLinArr< double > &ph_energy) const
DynLinArr< AtomicSecondaryProducts > asp
Definition: PhotoAbsCS.h:394
virtual void restore_shell(int nshell)
virtual double get_I_min(void) const
AtomicSecondaryProducts * get_asp(int nshell)
virtual double get_integral_ACS(double energy1, double energy2) const =0
int get_channel(DynLinArr< double > &felectron_energy, DynLinArr< double > &fphoton_energy) const
Definition: PhotoAbsCS.cpp:914
DynLinArr< DynLinArr< double > > photon_energy
Definition: PhotoAbsCS.h:319
DynLinArr< double > channel_prob_dens
Definition: PhotoAbsCS.h:316
virtual void print(std::ostream &file, int l) const
Definition: PhotoAbsCS.cpp:958
DynLinArr< DynLinArr< double > > electron_energy
Definition: PhotoAbsCS.h:318
void add_channel(double fchannel_prob_dens, const DynLinArr< double > &felectron_energy, const DynLinArr< double > &fphoton_energy, int s_all_rest=0)
Definition: PhotoAbsCS.cpp:875
virtual double get_threshold(int nshell) const
void replace_shells_by_overage(double fwidth, double fstep, long fmax_q_step)
virtual double get_integral_ACS(double energy1, double energy2) const
virtual double get_ACS(double energy) const
virtual int get_main_shell_number(int nshell) const
DynLinArr< ActivePtr< PhotoAbsCS > > acs
Definition: PhotoAbsCS.h:517
virtual void print(std::ostream &file, int l) const
virtual double get_ICS(double energy) const
virtual double get_integral_ICS(double energy1, double energy2) const
virtual void print(std::ostream &file, int l) const
Definition: PhotoAbsCS.cpp:406
virtual double get_CS(double energy) const
Definition: PhotoAbsCS.cpp:371
virtual double get_integral_CS(double energy1, double energy2) const
Definition: PhotoAbsCS.cpp:384
virtual void scale(double fact)
Definition: PhotoAbsCS.cpp:404
virtual double get_ACS(double energy) const
virtual double get_ICS(double energy) const
virtual double get_integral_ICS(double energy1, double energy2) const
virtual double get_integral_ACS(double energy1, double energy2) const
virtual void print(std::ostream &file, int l) const
virtual void print(std::ostream &file, int l) const
Definition: PhotoAbsCS.cpp:359
virtual double get_integral_CS(double energy1, double energy2) const
Definition: PhotoAbsCS.cpp:326
virtual void scale(double fact)
Definition: PhotoAbsCS.cpp:354
virtual double get_CS(double energy) const
Definition: PhotoAbsCS.cpp:308
virtual double get_integral_CS(double energy1, double energy2) const
Definition: PhotoAbsCS.cpp:757
virtual double get_CS(double energy) const
Definition: PhotoAbsCS.cpp:752
virtual void print(std::ostream &file, int l) const
Definition: PhotoAbsCS.cpp:776
virtual void scale(double fact)
Definition: PhotoAbsCS.cpp:771
virtual double get_ACS(double energy) const
virtual void print(std::ostream &file, int l) const
virtual int get_main_shell_number(int nshell) const
virtual double get_integral_ACS(double energy1, double energy2) const
virtual double get_integral_ICS(double energy1, double energy2) const
DynLinArr< ActivePtr< PhotoAbsCS > > acs
Definition: PhotoAbsCS.h:430
virtual double get_threshold(int nshell) const
virtual double get_ICS(double energy) const
virtual double get_integral_CS(double energy1, double energy2) const
Definition: PhotoAbsCS.cpp:672
void remove_leading_tiny(double level)
Definition: PhotoAbsCS.cpp:565
virtual double get_CS(double energy) const
Definition: PhotoAbsCS.cpp:585
virtual void print(std::ostream &file, int l) const
Definition: PhotoAbsCS.cpp:723
const DynLinArr< double > & get_arr_CS() const
Definition: PhotoAbsCS.h:211
const DynLinArr< double > & get_arr_ener() const
Definition: PhotoAbsCS.h:210
virtual void scale(double fact)
Definition: PhotoAbsCS.cpp:715
int findmark(std::istream &file, const char *s)
Definition: findmark.cpp:18
Definition: BGMesh.cpp:3
DynLinArr< double > make_log_mesh_ec(double emin, double emax, long q)
Definition: EnergyMesh.cpp:146
double my_integr_fun(double xp1, double yp1, double xp2, double yp2, double xmin, double, double x1, double x2)
Definition: PhotoAbsCS.cpp:49
const int s_scale_to_normalize_if_more
Definition: PhotoAbsCS.h:440
const double low_boundary_of_excitations
Definition: PhotoAbsCS.h:442
std::ostream & operator<<(std::ostream &file, const BGMesh &bgm)
Definition: BGMesh.cpp:22
const double coef_I_to_W
Definition: PhotoAbsCS.h:552
const double Thomas_sum_rule_const_Mb
Definition: PhotoAbsCS.h:45
double glin_integ_ar(DynLinArr< double > x, DynLinArr< double > y, long q, double x1, double x2, double threshold)
Definition: PhotoAbsCS.cpp:74
const int s_add_excitations_to_normalize
Definition: PhotoAbsCS.h:433
int sign_nonlinear_interpolation(double e1, double cs1, double e2, double cs2, double threshold)
Definition: PhotoAbsCS.cpp:22
double my_val_fun(double xp1, double yp1, double xp2, double yp2, double xmin, double, double x)
Definition: PhotoAbsCS.cpp:60
indentation indn
Definition: prstream.cpp:13
#define mcout
Definition: prstream.h:133
#define Ifile
Definition: prstream.h:207
#define Iprint(file, name)
Definition: prstream.h:209
#define mcerr
Definition: prstream.h:135
#define Iprintn(file, name)
Definition: prstream.h:216
#define Iprint2n(file, name1, name2)
Definition: prstream.h:236
ffloat SRANLUX(void)
Definition: ranluxint.h:262
T t_integ_generic_point_ar(const M &mesh, const D &y, T(*fun)(T xp1, T yp1, T xp2, T yp2, T xmin, T xmax, T x1, T x2), T x1, T x2, int s_extrap_left, T left_bond, int s_extrap_right, T right_bond)
Definition: tline.h:2786
T t_value_generic_point_ar(const M &mesh, const D &y, T(*funval)(T xp1, T yp1, T xp2, T yp2, T xmin, T xmax, T x), T x, int s_extrap_left, T left_bond, int s_extrap_right, T right_bond)
Definition: tline.h:2572