diff --git a/HEN_HOUSE/data/photo_shellwise.data b/HEN_HOUSE/data/photo_shellwise.data new file mode 100644 index 000000000..1d35d5083 Binary files /dev/null and b/HEN_HOUSE/data/photo_shellwise.data differ diff --git a/HEN_HOUSE/egs++/egs_advanced_application.cpp b/HEN_HOUSE/egs++/egs_advanced_application.cpp index 6050bdc07..206ccaef9 100644 --- a/HEN_HOUSE/egs++/egs_advanced_application.cpp +++ b/HEN_HOUSE/egs++/egs_advanced_application.cpp @@ -277,6 +277,9 @@ class EGS_LOCAL EGS_TransportProperty { options.push_back(opt); type = 2; }; + void setOption(const int &index, const char *opt) { + options[index] = opt; + }; void getInput(EGS_Input *input) { if (type == 1) { EGS_Float aux; @@ -496,9 +499,13 @@ int EGS_AdvancedApplication::helpInit(EGS_Input *transportp, bool do_hatch) { rayl.addOption("custom"); EGS_TransportProperty ff_med("ff media names",24,MXMED,&ff_media); EGS_TransportProperty ff_files("ff file names",128,MXMED, &ff_names); + EGS_TransportProperty relax("Atomic relaxations",&the_xoptions->iedgfl); relax.addOption("Off"); relax.addOption("On"); + relax.addOption("eadl"); + relax.addOption("simple"); + EGS_TransportProperty iphter("Photoelectron angular sampling", &the_xoptions->iphter); iphter.addOption("Off"); @@ -562,6 +569,16 @@ int EGS_AdvancedApplication::helpInit(EGS_Input *transportp, bool do_hatch) { setRayleighData(ff_media,ff_names); } relax.getInput(transportp); + if (the_xoptions->iedgfl == 0 || the_xoptions->iedgfl == 3) { + the_xoptions->eadl_relax = 0; + } + else { + the_xoptions->eadl_relax = 1; + // Now, 'On' means EADL relaxation + if (the_xoptions->iedgfl == 1) { + relax.setOption(1,"eadl"); + } + } iphter.getInput(transportp); spin.getInput(transportp); trip.getInput(transportp); @@ -576,6 +593,32 @@ int EGS_AdvancedApplication::helpInit(EGS_Input *transportp, bool do_hatch) { the_etcontrol->transport_algorithm = 0; } bca.getInput(transportp); + + if (egsEquivStr(string("mcdf-xcom "), + string(the_media->pxsec).substr(0,pxsec.size()))) { + the_xoptions->mcdf_pe_xsections = 1; + string("xcom ").copy(the_media->pxsec,16,0); + } + else if (egsEquivStr(string("mcdf-epdl "), + string(the_media->pxsec).substr(0,pxsec.size()))) { + the_xoptions->mcdf_pe_xsections = 1; + string("epdl ").copy(the_media->pxsec,16,0); + } + else { + the_xoptions->mcdf_pe_xsections = 0; + } + + // iedgfl == 3 implies eadl_relax == 0 + if (the_xoptions->iedgfl == 3 && + the_xoptions->mcdf_pe_xsections == 1) { + egsWarning("\n**** Warning:" + "\n Simplified atomic relaxation not allowed" + "\n with shellwise PE cross sections. Resetting" + "\n to detailed EADL atomic relaxation!!!\n\n"); + the_xoptions->eadl_relax = 1; + the_xoptions->iedgfl = 2; + relax.setOption(1,"eadl"); + } } if (do_hatch) { @@ -668,6 +711,19 @@ int EGS_AdvancedApplication::helpInit(EGS_Input *transportp, bool do_hatch) { egsInformation("\n\nTransport parameter and cross section options:\n" "==============================================\n"); int nc = 50; + // Recover initial input string to reflect what's actually being used + // the_media->pxsec was already used in egsHatch() to read photon data + if (the_xoptions->mcdf_pe_xsections && + egsEquivStr(string("xcom "), + string(the_media->pxsec).substr(0,pxsec.size()))) { + string("mcdf-xcom ").copy(the_media->pxsec,16,0); + } + else if (the_xoptions->mcdf_pe_xsections && + egsEquivStr(string("epdl "), + string(the_media->pxsec).substr(0,pxsec.size()))) { + string("mcdf-epdl ").copy(the_media->pxsec,16,0); + } + if (!isspace(the_media->pxsec[0])) { pxsec.info(nc); } @@ -715,7 +771,6 @@ int EGS_AdvancedApplication::helpInit(EGS_Input *transportp, bool do_hatch) { if (efield.size()==3 || bfield.size()==3) { estepem.info(nc); } - egsInformation("==============================================\n\n"); delete [] ind; @@ -776,7 +831,7 @@ void EGS_AdvancedApplication::setEIIData(EGS_I32 len) { the_xoptions->eii_flag = 0; } else if (str_eii == on) { - strcpy(the_media->eiixsec,ik.c_str()); + ik.copy(the_media->eiixsec,16,0); } } diff --git a/HEN_HOUSE/egs++/egs_functions.cpp b/HEN_HOUSE/egs++/egs_functions.cpp index 214daf2a0..013652422 100644 --- a/HEN_HOUSE/egs++/egs_functions.cpp +++ b/HEN_HOUSE/egs++/egs_functions.cpp @@ -23,7 +23,7 @@ # # Author: Iwan Kawrakow, 2005 # -# Contributors: +# Contributors: Ernesto Mainegra-Hing # ############################################################################### */ @@ -331,3 +331,14 @@ bool egsIsAbsolutePath(const string &path) { #endif } +bool egsEquivStr(const string &a, const string &b) { + unsigned int sz = a.size(); + if (b.size() != sz) { + return false; + } + for (unsigned int i = 0; i < sz; ++i) + if (tolower(a[i]) != tolower(b[i])) { + return false; + } + return true; +} diff --git a/HEN_HOUSE/egs++/egs_functions.h b/HEN_HOUSE/egs++/egs_functions.h index 8bb666f5f..5a361298a 100644 --- a/HEN_HOUSE/egs++/egs_functions.h +++ b/HEN_HOUSE/egs++/egs_functions.h @@ -23,7 +23,7 @@ # # Author: Iwan Kawrakow, 2005 # -# Contributors: +# Contributors: Ernesto Mainegra-Hing # ############################################################################### */ @@ -225,5 +225,6 @@ int EGS_EXPORT egsGetEndian(); */ bool EGS_EXPORT egsIsAbsolutePath(const string &path); +bool EGS_EXPORT egsEquivStr(const string &a, const string &b); #endif diff --git a/HEN_HOUSE/gui/egs_inprz/include/tooltips.h b/HEN_HOUSE/gui/egs_inprz/include/tooltips.h index 252cfa021..15027ac7d 100644 --- a/HEN_HOUSE/gui/egs_inprz/include/tooltips.h +++ b/HEN_HOUSE/gui/egs_inprz/include/tooltips.h @@ -442,21 +442,25 @@ static const char* sources[] = { #define BREMS_ANG_SAMPLING "If KM: complete modified Koch-Motz 2BS used." \ "\nIf Simple: only leading term used." -#define BOUND_COMPTON "If On, Compton scattering treated in the Impuls Approximation," \ +#define BOUND_COMPTON "If On, Compton scattering treated with the RIA model," \ "\nelse treated with Klein-Nishina. Make sure to turn it on for" \ "\nlow energy applications. Default is On." #define PE_ANG_SAMPLING "If Off, photo-electrons get direction of `mother' photon," \ "\nelse Sauter's formula used (which is, striktly speaking," \ "\nvalid only for K-shell photo-absorption)." -#define RAYLEIGH_SCAT "Turns on/off coherent (Rayleigh) scattering." \ - "\nDefault is Off. Should be turned on for low energy" \ - "\napplications but needs PEGS4 data set with data." +#define RAYLEIGH_SCAT "Turns on/off coherent (Rayleigh) scattering." \ + "\nShould be turned on for low energy applications." \ + "\nDefault is On." #define RELAXATION_TIP "If On, shell vacancies relaxed via emission of fluorescent X-Rays," \ "\nAuger and Koster-Cronig electrons (element and shell sampled from" \ "\nappropriate cross sections). Make sure to turn this option on for" \ - "\nlow energy applications. Default is On." + "\nlow energy applications. Default is EADL (On)."\ + "\n On: Full EADL atomic relaxation"\ + "\n Off: No atomic relaxation"\ + "\n EADL: Full EADL atomic relaxation"\ + "\n simple: Atomic relaxations for K, L1-L3,, and shells" #define ESTEPE_ALGORITHM "Algorithm used to take into account lateral and longitudinal" \ "\ncorrelations in a condensed history step." @@ -466,9 +470,11 @@ static const char* sources[] = { "\nICRU radiative stopping powers)." #define PHOTON_XSECTION "Select photon cross section compilation." \ - "\n si: Default photon cross section (Storm&Israel)" \ - "\n xcom: Data taken directly from the XCOM library" \ + "\n si: Data taken from the compilation by Storm&Israel" \ + "\n xcom: Data taken directly from the XCOM library (default)" \ "\n epdl: Data taken directly from the EPDL97 library" \ + "\n mcdf-xcom: Renormalized photoelectric xsections, XCOM for the rest" \ + "\n mcdf-epdl: Renormalized photoelectric xsections, EPDL for the rest" \ "\nPEGS4: Data taken directly from the PEGS4 file" \ "\n\nThe prefix for any other photon cross setion compilation" \ "\nplaced in $HEN_HOUSE/data/ in the form prefix_*.data" \ diff --git a/HEN_HOUSE/gui/egs_inprz/src/inputRZForm.cpp b/HEN_HOUSE/gui/egs_inprz/src/inputRZForm.cpp index ff14762ac..9c9b5bb64 100644 --- a/HEN_HOUSE/gui/egs_inprz/src/inputRZForm.cpp +++ b/HEN_HOUSE/gui/egs_inprz/src/inputRZForm.cpp @@ -3815,6 +3815,9 @@ for(QStringList::Iterator it = ini.begin(); it != ini.end(); ++it ) { prefix += (*it); } } +/* Insert prefixes requesting use of either xcom or epdl library + * with renormalized photoelectric xsections */ +prefix +="mcdf-xcom"; prefix +="mcdf-epdl"; /* Insert prefix "pegs4" which requests use of the pegs4 data */ prefix +="PEGS4"; diff --git a/HEN_HOUSE/gui/egs_inprz/ui/inputRZ.ui b/HEN_HOUSE/gui/egs_inprz/ui/inputRZ.ui index 4cd2faaab..68067cfe5 100644 --- a/HEN_HOUSE/gui/egs_inprz/ui/inputRZ.ui +++ b/HEN_HOUSE/gui/egs_inprz/ui/inputRZ.ui @@ -371,7 +371,7 @@ QTabWidget::Rounded - 5 + 7 @@ -8228,6 +8228,16 @@ Will start in /pegs4/density_corrections/ subdirectory off the directory selecte Off in Regions + + + eadl + + + + + simple + + diff --git a/HEN_HOUSE/interface/egs_c_interface2.macros b/HEN_HOUSE/interface/egs_c_interface2.macros index 32bf0e79c..5a1a73151 100644 --- a/HEN_HOUSE/interface/egs_c_interface2.macros +++ b/HEN_HOUSE/interface/egs_c_interface2.macros @@ -248,11 +248,13 @@ REPLACE {;COMIN/USEFUL/;} WITH { REPLACE {;COMIN/X-OPTIONS/;} WITH { ;common/xsection_options/ ibrdst, iprdst, ibr_nist, spin_effects, ibcmp,iraylr,iedgfl,iphter,pair_nrc,itriplet, - radc_flag,eii_flag,iphotonuc; + radc_flag,eii_flag,iphotonuc, + eadl_relax, mcdf_pe_xsections; $INTEGER ibrdst, iprdst, ibr_nist,eii_flag,iphotonuc, ibcmp,iraylr,iedgfl,iphter,pair_nrc, radc_flag,itriplet; $LOGICAL spin_effects; + $LOGICAL eadl_relax, mcdf_pe_xsections; }; " Again, to remove geometry dependence make e_max_rr and i_do_rr sclars" @@ -293,7 +295,8 @@ REPLACE {$set-region-by-region-defaults;} WITH {; ecut = 0; pcut = 0; smaxir = $MAX-SMAX; ecut_new = 0; pcut_new = 0; smax_new = $MAX-SMAX; rhor = 1; rhor_new = 1; - ibcmp = $IBCMP-DEFAULT; iraylr = 0; iphotonuc=$IPHOTONUCR-DEFAULT; + ibcmp = $IBCMP-DEFAULT; + iraylr = $IRAYLR-DEFAULT; iphotonuc=$IPHOTONUCR-DEFAULT; iedgfl = $IEDGFL-DEFAULT; iphter = $IPHTER-DEFAULT; i_do_rr = 0; e_max_rr = 0; e_max_rr_new = 0; }; @@ -512,7 +515,7 @@ APPEND {;COMIN/X-OPTIONS/;} TO {$COMIN-PHOTO;}; APPEND {;COMIN/X-OPTIONS/;} TO {$COMIN-PHOTON;}; APPEND {;COMIN/X-OPTIONS/;} TO {$COMIN-EII-INIT;}; APPEND {;COMIN/X-OPTIONS/;} TO {$COMIN-INIT-TRIPLET;}; -APPEND {;COMIN/X-OPTIONS/;} TO {$COMIN-GET-TRANSPORTP;}; +"APPEND {;COMIN/X-OPTIONS/;} TO {$COMIN-GET-TRANSPORTP;}; already there" APPEND {;COMIN/X-OPTIONS/;} TO {$COMIN-MOLLER;}; APPEND {;COMIN/X-OPTIONS/;} TO {$COMIN-RADC-INIT;}; APPEND {;COMIN/X-OPTIONS/;} TO {$COMIN-INIT-NIST-BREMS;}; diff --git a/HEN_HOUSE/interface/egs_interface2.h b/HEN_HOUSE/interface/egs_interface2.h index 9300960bd..4ec9c306b 100644 --- a/HEN_HOUSE/interface/egs_interface2.h +++ b/HEN_HOUSE/interface/egs_interface2.h @@ -411,7 +411,7 @@ struct EGS_XOptions { Can be set in the input file using \verbatim - Pair angular samplin= Off or Simple or KM + Pair angular sampling= Off or Simple or KM \endverbatim */ EGS_I32 iprdst; @@ -466,7 +466,7 @@ struct EGS_XOptions { /*! Atomic relaxations flag. If set to 1 (the default), vacancies created in shells with binding energies above 1 keV will be relaxed - by fluorscent, Auger and Coster-Kronig transitions. Vacancies + by fluorescent, Auger and Coster-Kronig transitions. Vacancies can be currently created after photo-absorption, bound Compton scattering and electron inelastic scattering (if eii_flag is not 0). If set to 0, the binding energy will be given to the photo-electron @@ -546,6 +546,10 @@ struct EGS_XOptions { */ EGS_I32 iphotonuc; + EGS_I32 eadl_relax; + + EGS_I32 mcdf_pe_xsections; + }; /*! \brief A structure corresponding to the \c egs_vr common block diff --git a/HEN_HOUSE/omega/progs/gui/beamnrc/beamnrc_gui.tcl b/HEN_HOUSE/omega/progs/gui/beamnrc/beamnrc_gui.tcl index 494e8253b..8d971a41e 100755 --- a/HEN_HOUSE/omega/progs/gui/beamnrc/beamnrc_gui.tcl +++ b/HEN_HOUSE/omega/progs/gui/beamnrc/beamnrc_gui.tcl @@ -300,6 +300,10 @@ global numopts options } # Place PEGS4 option to end of the list set options(photon_xsections,0) $xsect + set options(photon_xsections,[expr $numopts(photon_xsections)-1]) "mcdf-xcom" + set numopts(photon_xsections) [expr $numopts(photon_xsections)+1] + set options(photon_xsections,[expr $numopts(photon_xsections)-1]) "mcdf-epdl" + set numopts(photon_xsections) [expr $numopts(photon_xsections)+1] set options(photon_xsections,[expr $numopts(photon_xsections)-1]) "PEGS4" } else { # If for some reason no photon data available, one can still use the diff --git a/HEN_HOUSE/omega/progs/gui/beamnrc/beamnrc_params.tcl b/HEN_HOUSE/omega/progs/gui/beamnrc/beamnrc_params.tcl index 015e1b619..5b5d1ac62 100644 --- a/HEN_HOUSE/omega/progs/gui/beamnrc/beamnrc_params.tcl +++ b/HEN_HOUSE/omega/progs/gui/beamnrc/beamnrc_params.tcl @@ -1351,15 +1351,19 @@ set names(photon_xsections) "Photon cross-sections" set help_text(photon_xsections) { Photon cross-sections (photon_xsections): -si (default), epdl, xcom, PEGS4, etc. \ +si, epdl, xcom (default), mcdf-xcom, mcdf-epdl, PEGS4, etc. \ Determines which photon cross-sections are used. -By default EGSnrc uses the Storm-Israel cross-sections (si_*.data files)\ +By default EGSnrc uses the XCOM cross-sections (xcom_*.data files)\ . Other photon cross section compilations provided with EGSnrc are the \ -EPDL (epdl_*.data files) and XCOM (xcom_*.data files) photon cross section\ +EPDL (epdl_*.data files) and Storm-Israel (si_*.data files) photon cross section\ compilations. For more details about these options, see \ the EGSnrc Manual (PIRS-701). +Selecting either mcdf-xcom or mcdf-epdl allows the use of renormalized\ +photoelectric cross sections by Sabbatucci and Salvat with either XCOM\ +or EDPL cross sections for the rest of the photon interactions. + Selecting PEGS4 forces the use of the photon cross section data in the\ PEGS4 file. This recently added option allows the use of PEGS4 files\ if the need arises. For instance, for comparison with previous calculations.\ @@ -1470,7 +1474,7 @@ uses actual bound Compton cross sections and does not\ reject any Compton interactions at run time. If the Norej\ option is selected, then the user has the option\ of specifying their own Compton cross sections.\ -Default is Off\ +Default is Norej\ Make sure to use bound Compton for low energy applications.\ Not necessary above, say, 1 MeV.\ @@ -1617,7 +1621,7 @@ Rayleigh scattering (IRAYLR): Off, On, On in regions, Off in regions, custom\ If On, turn on coherent (Rayleigh) scattering.\ -Default is Off. Should be turned on for low energy\ +Default is On. Should be turned on for low energy\ applications. Not set to On by default because\ On requires a sperial PEGS4 data set.\ @@ -1636,7 +1640,7 @@ set options(iraylr,1) "On" set options(iraylr,2) "On in regions" set options(iraylr,3) "Off in regions" set options(iraylr,4) "custom" -set values(iraylr) $options(iraylr,0) +set values(iraylr) $options(iraylr,1) set level(iraylr) 0 set num_rayl_custom 0 for {set i 1} {$i <= 200} {incr i} { @@ -1649,33 +1653,43 @@ set names(iedgfl) "Atomic relaxations" set help_text(iedgfl) { Atomic relaxations (IEDGFL): -Off, On, On in regions, Off in regions\ -The default is Off. The effect of using On is twofold:\ -- In photo-electric absorption events, the element\ - (if material is mixture) and the shell the photon\ +Off, On, On in regions, Off in regions, EADL, simple. \ +The default is EADL (On). The effect of using On is twofold:\ +- In photo-electric absorption, incoherent scattering with\ + bound electrons, and electron impact ionization events, \ + the element (if material is mixture) and the shell the particle\ is interacting with are sampled from the appropriate\ - cross seections\ -- Shell vacancies created in photo-absorption events\ + cross sections\ +- Shell vacancies created in these events\ are relaxed via emission of fluorescent X-Rays,\ Auger and Koster-Cronig electrons.\ Make sure to turn this option on for low energy\ applications.\ +Selecting On or EADl is equivalent. A full atomic relaxation cascade\ +is simulated using EADL transition probabilities. + +The option to use the original algorithm only accounting for M and\ +N shells in an average way is left for comparison and can be invoked\ +by selecting the option simple. + If you use "On in regions", you will be prompted to enter\ pairs of region numbers. Each pair defines a range of regions\ in which Atomic relaxations will be turned On. Everywhere\ -outside these regions Atomic relaxations will be turned Off.\ +outside these regions Atomic relaxations will be turned Off. If you use "Off in regions", you will be prompted to enter\ ranges of regions where Atomic relaxation is to be turned Off.\ Everywhere else it will be turned on. } -set numopts(iedgfl) 4 +set numopts(iedgfl) 6 set options(iedgfl,0) "Off" set options(iedgfl,1) "On" set options(iedgfl,2) "On in regions" set options(iedgfl,3) "Off in regions" -set values(iedgfl) $options(iedgfl,0) +set options(iedgfl,4) "EADL" +set options(iedgfl,5) "simple" +set values(iedgfl) $options(iedgfl,4) set level(iedgfl) 0 #################end of options for EGSnrc parameters#################### diff --git a/HEN_HOUSE/omega/progs/gui/dosxyznrc/dosxyznrc_gui.tcl b/HEN_HOUSE/omega/progs/gui/dosxyznrc/dosxyznrc_gui.tcl index 76c7d2d2a..0d88a2947 100755 --- a/HEN_HOUSE/omega/progs/gui/dosxyznrc/dosxyznrc_gui.tcl +++ b/HEN_HOUSE/omega/progs/gui/dosxyznrc/dosxyznrc_gui.tcl @@ -263,6 +263,10 @@ global numopts options } # Place PEGS4 option to end of the list set options(photon_xsections,0) $xsect + set options(photon_xsections,[expr $numopts(photon_xsections)-1]) "mcdf-xcom" + set numopts(photon_xsections) [expr $numopts(photon_xsections)+1] + set options(photon_xsections,[expr $numopts(photon_xsections)-1]) "mcdf-epdl" + set numopts(photon_xsections) [expr $numopts(photon_xsections)+1] set options(photon_xsections,[expr $numopts(photon_xsections)-1]) "PEGS4" } else { # If for some reason no photon data available, one can still use the diff --git a/HEN_HOUSE/omega/progs/gui/dosxyznrc/xyznrc_parameters.tcl b/HEN_HOUSE/omega/progs/gui/dosxyznrc/xyznrc_parameters.tcl index 67508490d..7dee506c2 100644 --- a/HEN_HOUSE/omega/progs/gui/dosxyznrc/xyznrc_parameters.tcl +++ b/HEN_HOUSE/omega/progs/gui/dosxyznrc/xyznrc_parameters.tcl @@ -951,7 +951,7 @@ set names(photon_xsections) "Photon cross-sections" set help_text(photon_xsections) { Photon cross-sections (photon_xsections): -si , epdl, xcom (default), PEGS4, etc. \ +si, epdl, xcom (default), mcdf-xcom, mcdf-epdl, PEGS4, etc. \ Determines which photon cross-sections are used. By default EGSnrc uses the NIST XCOM cross-sections (xcom_*.data files)\ @@ -960,6 +960,10 @@ EPDL (epdl_*.data files) and Strom & Israel (si_*.data files) photon cross secti compilations. For more details about these options, see \ the EGSnrc Manual (PIRS-701). +Selecting either mcdf-xcom or mcdf-epdl allows the use of renormalized\ +photoelectric cross sections by Sabbatucci and Salvat with either XCOM\ +or EDPL cross sections for the rest of the photon interactions. + Selecting PEGS4 forces the use of the photon cross section data in the\ PEGS4 file. This recently added option allows the use of PEGS4 files\ if the need arises. For instance, for comparison with previous calculations.\ @@ -1067,7 +1071,7 @@ uses actual bound Compton cross sections and does not\ reject any Compton interactions at run time. If the Norej\ option is selected, then the user has the option\ of specifying their own Compton cross sections.\ -Default is Off\ +Default is Norej\ Make sure to use bound Compton for low energy applications.\ Not necessary above, say, 1 MeV.\ @@ -1088,7 +1092,7 @@ set options(ibcmp,1) "On" set options(ibcmp,2) "On in regions" set options(ibcmp,3) "Off in regions" set options(ibcmp,4) "Norej" -set values(ibcmp) $options(ibcmp,0) +set values(ibcmp) $options(ibcmp,4) set level(ibcmp) 0 set comp_xsections {} @@ -1193,7 +1197,7 @@ Rayleigh scattering (IRAYLR): Off, On, On in regions, Off in regions, custom\ If On, turn on coherent (Rayleigh) scattering.\ -Default is Off. Should be turned on for low energy\ +Default is On. Should be turned on for low energy\ applications. Not set to On by default because\ On requires a sperial PEGS4 data set.\ @@ -1212,23 +1216,34 @@ set options(iraylr,1) "On" set options(iraylr,2) "On in regions" set options(iraylr,3) "Off in regions" set options(iraylr,4) "custom" +set values(iraylr) $options(iraylr,1) set iedgfl {} set names(iedgfl) "Atomic relaxations" set help_text(iedgfl) { Atomic relaxations (IEDGFL): -Off, On, On in regions, Off in regions\ -The default is Off. The effect of using On is twofold:\ -- In photo-electric absorption events, the element\ - (if material is mixture) and the shell the photon\ +Off, On, On in regions, Off in regions, EADL, simple. + +The default is EADL (On). The effect of using On is twofold: + +- In photo-electric absorption, incoherent scattering with\ + bound electrons, and electron impact ionization events, \ + the element (if material is mixture) and the shell the particle\ is interacting with are sampled from the appropriate\ - cross seections\ -- Shell vacancies created in photo-absorption events\ + cross sections +- Shell vacancies created in these events\ are relaxed via emission of fluorescent X-Rays,\ Auger and Koster-Cronig electrons.\ Make sure to turn this option on for low energy\ - applications.\ + applications. + +Selecting On or EADl is equivalent. A full atomic relaxation cascade\ +is simulated using EADL transition probabilities. + +The option to use the original algorithm only accounting for M and\ +N shells in an average way is left for comparison and can be invoked\ +by selecting the option simple. If you use "On in regions", you will be prompted to enter\ pairs of region numbers. Each pair defines a range of regions\ @@ -1239,11 +1254,14 @@ If you use "Off in regions", you will be prompted to enter\ ranges of regions where Atomic relaxation is to be turned Off.\ Everywhere else it will be turned on. } -set numopts(iedgfl) 4 +set numopts(iedgfl) 6 set options(iedgfl,0) "Off" set options(iedgfl,1) "On" set options(iedgfl,2) "On in regions" set options(iedgfl,3) "Off in regions" +set options(iedgfl,4) "EADL" +set options(iedgfl,5) "simple" +set values(iedgfl) $options(iedgfl,4) #################end of options for EGSnrc parameters#################### diff --git a/HEN_HOUSE/pegs4/density_corrections/compounds/water_ICRU90.density b/HEN_HOUSE/pegs4/density_corrections/compounds/water_ICRU90.density new file mode 100644 index 000000000..639717a5c --- /dev/null +++ b/HEN_HOUSE/pegs4/density_corrections/compounds/water_ICRU90.density @@ -0,0 +1,32 @@ +water_icru90 + 113 78.00000 1.000000 2 + 1 0.111894 8 0.888106 + 1.00E-03,0.000E+00 1.25E-03,0.000E+00 1.50E-03,0.000E+00 1.75E-03,0.000E+00 + 2.00E-03,0.000E+00 2.50E-03,0.000E+00 3.00E-03,0.000E+00 3.50E-03,0.000E+00 + 4.00E-03,0.000E+00 4.50E-03,0.000E+00 5.00E-03,0.000E+00 5.50E-03,0.000E+00 + 6.00E-03,0.000E+00 7.00E-03,0.000E+00 8.00E-03,0.000E+00 9.00E-03,0.000E+00 + 1.00E-02,0.000E+00 1.25E-02,0.000E+00 1.50E-02,0.000E+00 1.75E-02,0.000E+00 + 2.00E-02,0.000E+00 2.50E-02,0.000E+00 3.00E-02,0.000E+00 3.50E-02,0.000E+00 + 4.00E-02,0.000E+00 4.50E-02,0.000E+00 5.00E-02,0.000E+00 5.50E-02,0.000E+00 + 6.00E-02,0.000E+00 7.00E-02,0.000E+00 8.00E-02,0.000E+00 9.00E-02,0.000E+00 + 1.00E-01,0.000E+00 1.25E-01,0.000E+00 1.50E-01,0.000E+00 1.75E-01,0.000E+00 + 2.00E-01,0.000E+00 2.50E-01,0.000E+00 3.00E-01,0.000E+00 3.50E-01,0.000E+00 + 4.00E-01,0.000E+00 4.50E-01,0.000E+00 5.00E-01,0.000E+00 5.50E-01,8.304E-04 + 6.00E-01,1.533E-02 7.00E-01,5.387E-02 8.00E-01,1.011E-01 9.00E-01,1.537E-01 + 1.00E+00,2.094E-01 1.25E+00,3.545E-01 1.50E+00,4.993E-01 1.75E+00,6.389E-01 + 2.00E+00,7.716E-01 2.50E+00,1.015E+00 3.00E+00,1.233E+00 3.50E+00,1.427E+00 + 4.00E+00,1.602E+00 4.50E+00,1.762E+00 5.00E+00,1.908E+00 5.50E+00,2.042E+00 + 6.00E+00,2.167E+00 7.00E+00,2.392E+00 8.00E+00,2.590E+00 9.00E+00,2.768E+00 + 1.00E+01,2.930E+00 1.25E+01,3.278E+00 1.50E+01,3.568E+00 1.75E+01,3.819E+00 + 2.00E+01,4.040E+00 2.50E+01,4.419E+00 3.00E+01,4.736E+00 3.50E+01,5.011E+00 + 4.00E+01,5.254E+00 4.50E+01,5.471E+00 5.00E+01,5.667E+00 5.50E+01,5.847E+00 + 6.00E+01,6.012E+00 7.00E+01,6.307E+00 8.00E+01,6.565E+00 9.00E+01,6.794E+00 + 1.00E+02,7.000E+00 1.25E+02,7.438E+00 1.50E+02,7.798E+00 1.75E+02,8.104E+00 + 2.00E+02,8.369E+00 2.50E+02,8.812E+00 3.00E+02,9.175E+00 3.50E+02,9.483E+00 + 4.00E+02,9.749E+00 4.50E+02,9.984E+00 5.00E+02,1.019E+01 5.50E+02,1.038E+01 + 6.00E+02,1.056E+01 7.00E+02,1.087E+01 8.00E+02,1.113E+01 9.00E+02,1.137E+01 + 1.00E+03,1.158E+01 1.25E+03,1.203E+01 1.50E+03,1.239E+01 1.75E+03,1.270E+01 + 2.00E+03,1.296E+01 2.50E+03,1.341E+01 3.00E+03,1.378E+01 3.50E+03,1.408E+01 + 4.00E+03,1.435E+01 4.50E+03,1.459E+01 5.00E+03,1.480E+01 5.50E+03,1.499E+01 + 6.00E+03,1.516E+01 7.00E+03,1.547E+01 8.00E+03,1.574E+01 9.00E+03,1.597E+01 + 1.00E+04,1.618E+01 diff --git a/HEN_HOUSE/pegs4/density_corrections/elements/carbon_graphite_ICRU90_1.7g_cm3.density b/HEN_HOUSE/pegs4/density_corrections/elements/carbon_graphite_ICRU90_1.7g_cm3.density new file mode 100644 index 000000000..b8902312b --- /dev/null +++ b/HEN_HOUSE/pegs4/density_corrections/elements/carbon_graphite_ICRU90_1.7g_cm3.density @@ -0,0 +1,32 @@ +c2.265 + 113 81.00000 1.7000 1 + 6 1.000000 + 1.00E-03,2.261E-04 1.25E-03,2.835E-04 1.50E-03,3.413E-04 1.75E-03,3.994E-04 + 2.00E-03,4.579E-04 2.50E-03,5.759E-04 3.00E-03,6.953E-04 3.50E-03,8.161E-04 + 4.00E-03,9.383E-04 4.50E-03,1.062E-03 5.00E-03,1.187E-03 5.50E-03,1.313E-03 + 6.00E-03,1.441E-03 7.00E-03,1.700E-03 8.00E-03,1.965E-03 9.00E-03,2.235E-03 + 1.00E-02,2.511E-03 1.25E-02,3.223E-03 1.50E-02,3.967E-03 1.75E-02,4.742E-03 + 2.00E-02,5.549E-03 2.50E-02,7.251E-03 3.00E-02,9.070E-03 3.50E-02,1.100E-02 + 4.00E-02,1.304E-02 4.50E-02,1.518E-02 5.00E-02,1.741E-02 5.50E-02,1.974E-02 + 6.00E-02,2.217E-02 7.00E-02,2.727E-02 8.00E-02,3.269E-02 9.00E-02,3.841E-02 + 1.00E-01,4.440E-02 1.25E-01,6.044E-02 1.50E-01,7.778E-02 1.75E-01,9.618E-02 + 2.00E-01,1.154E-01 2.50E-01,1.559E-01 3.00E-01,1.982E-01 3.50E-01,2.415E-01 + 4.00E-01,2.854E-01 4.50E-01,3.293E-01 5.00E-01,3.731E-01 5.50E-01,4.166E-01 + 6.00E-01,4.595E-01 7.00E-01,5.434E-01 8.00E-01,6.245E-01 9.00E-01,7.025E-01 + 1.00E+00,7.774E-01 1.25E+00,9.519E-01 1.50E+00,1.110E+00 1.75E+00,1.253E+00 + 2.00E+00,1.384E+00 2.50E+00,1.617E+00 3.00E+00,1.818E+00 3.50E+00,1.996E+00 + 4.00E+00,2.156E+00 4.50E+00,2.301E+00 5.00E+00,2.433E+00 5.50E+00,2.556E+00 + 6.00E+00,2.671E+00 7.00E+00,2.879E+00 8.00E+00,3.067E+00 9.00E+00,3.237E+00 + 1.00E+01,3.394E+00 1.25E+01,3.742E+00 1.50E+01,4.041E+00 1.75E+01,4.304E+00 + 2.00E+01,4.539E+00 2.50E+01,4.943E+00 3.00E+01,5.282E+00 3.50E+01,5.574E+00 + 4.00E+01,5.829E+00 4.50E+01,6.057E+00 5.00E+01,6.261E+00 5.50E+01,6.447E+00 + 6.00E+01,6.617E+00 7.00E+01,6.920E+00 8.00E+01,7.183E+00 9.00E+01,7.416E+00 + 1.00E+02,7.624E+00 1.25E+02,8.067E+00 1.50E+02,8.429E+00 1.75E+02,8.736E+00 + 2.00E+02,9.002E+00 2.50E+02,9.447E+00 3.00E+02,9.811E+00 3.50E+02,1.012E+01 + 4.00E+02,1.039E+01 4.50E+02,1.062E+01 5.00E+02,1.083E+01 5.50E+02,1.102E+01 + 6.00E+02,1.119E+01 7.00E+02,1.150E+01 8.00E+02,1.177E+01 9.00E+02,1.201E+01 + 1.00E+03,1.222E+01 1.25E+03,1.266E+01 1.50E+03,1.303E+01 1.75E+03,1.333E+01 + 2.00E+03,1.360E+01 2.50E+03,1.405E+01 3.00E+03,1.441E+01 3.50E+03,1.472E+01 + 4.00E+03,1.499E+01 4.50E+03,1.522E+01 5.00E+03,1.543E+01 5.50E+03,1.562E+01 + 6.00E+03,1.580E+01 7.00E+03,1.611E+01 8.00E+03,1.637E+01 9.00E+03,1.661E+01 + 1.00E+04,1.682E+01 diff --git a/HEN_HOUSE/pegs4/density_corrections/elements/carbon_graphite_ICRU90_2.0g_cm3.density b/HEN_HOUSE/pegs4/density_corrections/elements/carbon_graphite_ICRU90_2.0g_cm3.density new file mode 100644 index 000000000..8825d36ad --- /dev/null +++ b/HEN_HOUSE/pegs4/density_corrections/elements/carbon_graphite_ICRU90_2.0g_cm3.density @@ -0,0 +1,32 @@ +c2.265 + 113 81.00000 2.000 1 + 6 1.000000 + 1.00E-03,2.261E-04 1.25E-03,2.835E-04 1.50E-03,3.413E-04 1.75E-03,3.994E-04 + 2.00E-03,4.579E-04 2.50E-03,5.759E-04 3.00E-03,6.953E-04 3.50E-03,8.161E-04 + 4.00E-03,9.383E-04 4.50E-03,1.062E-03 5.00E-03,1.187E-03 5.50E-03,1.313E-03 + 6.00E-03,1.441E-03 7.00E-03,1.700E-03 8.00E-03,1.965E-03 9.00E-03,2.235E-03 + 1.00E-02,2.511E-03 1.25E-02,3.223E-03 1.50E-02,3.967E-03 1.75E-02,4.742E-03 + 2.00E-02,5.549E-03 2.50E-02,7.251E-03 3.00E-02,9.070E-03 3.50E-02,1.100E-02 + 4.00E-02,1.304E-02 4.50E-02,1.518E-02 5.00E-02,1.741E-02 5.50E-02,1.974E-02 + 6.00E-02,2.217E-02 7.00E-02,2.727E-02 8.00E-02,3.269E-02 9.00E-02,3.841E-02 + 1.00E-01,4.440E-02 1.25E-01,6.044E-02 1.50E-01,7.778E-02 1.75E-01,9.618E-02 + 2.00E-01,1.154E-01 2.50E-01,1.559E-01 3.00E-01,1.982E-01 3.50E-01,2.415E-01 + 4.00E-01,2.854E-01 4.50E-01,3.293E-01 5.00E-01,3.731E-01 5.50E-01,4.166E-01 + 6.00E-01,4.595E-01 7.00E-01,5.434E-01 8.00E-01,6.245E-01 9.00E-01,7.025E-01 + 1.00E+00,7.774E-01 1.25E+00,9.519E-01 1.50E+00,1.110E+00 1.75E+00,1.253E+00 + 2.00E+00,1.384E+00 2.50E+00,1.617E+00 3.00E+00,1.818E+00 3.50E+00,1.996E+00 + 4.00E+00,2.156E+00 4.50E+00,2.301E+00 5.00E+00,2.433E+00 5.50E+00,2.556E+00 + 6.00E+00,2.671E+00 7.00E+00,2.879E+00 8.00E+00,3.067E+00 9.00E+00,3.237E+00 + 1.00E+01,3.394E+00 1.25E+01,3.742E+00 1.50E+01,4.041E+00 1.75E+01,4.304E+00 + 2.00E+01,4.539E+00 2.50E+01,4.943E+00 3.00E+01,5.282E+00 3.50E+01,5.574E+00 + 4.00E+01,5.829E+00 4.50E+01,6.057E+00 5.00E+01,6.261E+00 5.50E+01,6.447E+00 + 6.00E+01,6.617E+00 7.00E+01,6.920E+00 8.00E+01,7.183E+00 9.00E+01,7.416E+00 + 1.00E+02,7.624E+00 1.25E+02,8.067E+00 1.50E+02,8.429E+00 1.75E+02,8.736E+00 + 2.00E+02,9.002E+00 2.50E+02,9.447E+00 3.00E+02,9.811E+00 3.50E+02,1.012E+01 + 4.00E+02,1.039E+01 4.50E+02,1.062E+01 5.00E+02,1.083E+01 5.50E+02,1.102E+01 + 6.00E+02,1.119E+01 7.00E+02,1.150E+01 8.00E+02,1.177E+01 9.00E+02,1.201E+01 + 1.00E+03,1.222E+01 1.25E+03,1.266E+01 1.50E+03,1.303E+01 1.75E+03,1.333E+01 + 2.00E+03,1.360E+01 2.50E+03,1.405E+01 3.00E+03,1.441E+01 3.50E+03,1.472E+01 + 4.00E+03,1.499E+01 4.50E+03,1.522E+01 5.00E+03,1.543E+01 5.50E+03,1.562E+01 + 6.00E+03,1.580E+01 7.00E+03,1.611E+01 8.00E+03,1.637E+01 9.00E+03,1.661E+01 + 1.00E+04,1.682E+01 diff --git a/HEN_HOUSE/src/egs_utilities.mortran b/HEN_HOUSE/src/egs_utilities.mortran index cbd77fb02..cd3963791 100644 --- a/HEN_HOUSE/src/egs_utilities.mortran +++ b/HEN_HOUSE/src/egs_utilities.mortran @@ -896,12 +896,14 @@ $set-region-by-region-defaults; eii_flag = 0; "No EII by default. " eii_xfile = 'Off'; eii_L_factor = 1.0; "No L-shell EII xsection scaling by default" - -xsec_out = $XSEC-DEFAULT; "see egsnrc.macros for default" - -"DO i=1,len(photon_xsections) [ photon_xsections(i:i) = ' '; ]" -photon_xsections = $XDATA-DEFAULT;"default photon xsection, see egsnrc.macros" +"=========================================" +"See egsnrc.macros for defaults used below" +"=========================================" +xsec_out = $XSEC-DEFAULT; +photon_xsections = $XDATA-DEFAULT;"default photon xsection" comp_xsections = $COMP-XDATA-DEFAULT; +eadl_relax = $EADL-RELAX-DEFAULT; +mcdf_pe_xsections = $MCDF-PE-DEFAULT; "Ali:photonuc, 2 lines" photonuc_xsections = $PHOTONUC-XDATA-DEFAULT; "EMH:emf" @@ -913,9 +915,8 @@ BxIN=$BxDEF;ByIN=$ByDEF;BzIN=$BzDEF; EMLMTIN=$EMLMTDEF; Bx=BxIN; By=ByIN; Bz=BzIN; Bx_new=Bx; By_new=By; Bz_new=Bz; -"DO i=1,len(eii_xfile) [ eii_xfile(i:i) = ' '; ]" DO i=1,$MXMED [ - iraylm(i) = 0; "Rayleigh data available?" + iraylm(i) = 0; "Rayleigh data available?" DO j=1,len(iray_ff_file(i)) [ iray_ff_file(i)(j:j) = ' ';] DO j=1,len(iray_ff_media(i)) [ iray_ff_media(i)(j:j) = ' ';] " set all thresholds to zero " @@ -1629,6 +1630,39 @@ IF( which = 1 ) [ ] return; end; +/* Print binding energies: In the case of the EPDL library + only energies above 1 keV are output for elements not part + of the current simlation. This is due to the fact that energies + below 1 keV are taken from the relaxation database only for + elements requested in the input when EPDL is used. This is not + an issue when using XCOM. + */ +subroutine egs_print_binding_energies; +implicit none; +$declare_max_medium; +;COMIN/EDGE,MEDIA,EGS-IO/; +$INTEGER i,j; +integer*4 lnblnk1; +character*3 labels(16); +data labels/' K',' L1',' L2',' L3', + ' M1',' M2',' M3',' M4',' M5', + ' N1',' N2',' N3',' N4',' N5',' N6',' N7'/; + +$egs_info('(a,a,a)', + 'Binding energies from ',$cstring(photon_xsections), + ' photon cross section library'); +DO j = 1,$MXELEMENT [ + DO i = 1,$MXPESHELL [ + IF ( binding_energies(i,j) > 0 ) [ + $egs_info('(a,i3,a,a,a,1pe12.4,a)', + ' Eb(',j,',',labels(i),') = ',binding_energies(i,j),' MeV'); + ] + ] +] + +return;end; + + "=============================================================================" " scale elastic scattering strength by a given factor " "=============================================================================" @@ -2227,3 +2261,25 @@ return; end; jrec = jrec + 4; egs_read_real = r_4; return; end; +"**************************************************************** +"* * +"* Function ibsearch(a, nsh, b) * +"* * +"* binary search for an element l of array b such that * +"* b[l] =< a < b[l+1], array must be monotonically increasing * +"* * +"**************************************************************** +$INTEGER function ibsearch(a, nsh, b); + implicit none; + $REAL a, b(*); + $INTEGER min,max,help,nsh; + $REAL x; + min = 1; max = nsh; x = a; + WHILE ( min < max-1 )[ + help = (max+min)/2; + IF ( b(help).le.x)[min = help;] + ELSE[max = help;] + ] + ibsearch = min; + return;end; + diff --git a/HEN_HOUSE/src/egsnrc.macros b/HEN_HOUSE/src/egsnrc.macros index bf87f55c9..e857a9518 100644 --- a/HEN_HOUSE/src/egsnrc.macros +++ b/HEN_HOUSE/src/egsnrc.macros @@ -673,25 +673,23 @@ REPLACE {$MAXSHELL} WITH {3000}"max. number of shells" REPLACE {$MAXRELAX} WITH {10000}"max. number of relaxations channels" REPLACE {$MAXVAC} WITH {100} "max. number of vacancies" "============================================================" -" Set the macro below to .false. to recover the previous approach -" which allows photoelectric interactions with and shells. -" See below for more details on the shells considered by different -" interactions depending on the value of $EADL_RELAX: +" Set input key 'Atomic relaxations' to 'simple' to recover original +" implementation which allows photoelectric interactions with and +" shells. See below for details on the shells considered by different +" interactions depending on the value of eadl_relax: " " Interaction .false. .true. " ----------------------------------------- " Compton all available shells " EII K,L1..L3 K,L1..L3 " Photoeffect K,L1..L3,, K,L1..L3 +" Shellwise +" Photoeffect N/A All shells > $RELAX-CUTOFF " Relaxation " initial vacancy K,L1..L3, K,L1..L3 +" (for new photoeffect) K, L1..L3, M1..M5, N1..N4 " final vacancy L1..L3,, L1..L3,M1..M5,N1..N7... " -" The current limitation of the initial vacancies will be removed -" when a detailed implementation of photoelectric and EII interactions -" becomes available. -"============================================================" -REPLACE {$EADL_RELAX} WITH {.true.} "============================================================" ; "---------- BUFFER FLUSH SEMICOLON ----------" @@ -747,10 +745,57 @@ RELAX-DATA,SHELL-DATA/; REPLACE {$COMIN-RELAX-EADL;} WITH { ;COMIN/RELAX-DATA,RELAX-FOR-USER,SHELL-DATA, -STACK,THRESH,EPCONT,USEFUL,UPHIOT,RANDOM,BOUNDS,EGS-IO,MISC,MEDIA/; +STACK,THRESH,EPCONT,USEFUL,UPHIOT,RANDOM,BOUNDS,EGS-IO,MISC,MEDIA, +X-OPTIONS/; +}; +; "---------- BUFFER FLUSH SEMICOLON ----------" + +"***************************************************************************" +" " +" -------------- shell-wise photoelectric cross section data ------------ " +" +" Cross sections taken from Sabbatucci and Salvat, " +" Theory and calculation of the atomic photoeffect " +" Radiat. Phys. and Chem., Volume 121, April 2016, Pages 122-140 " +" " +"***************************************************************************" +" Shell-wise photoelectric cross sections from 50 eV up to 1 GeV for elements +" from Z=1 up to Z=99. All K, L, M, and N shells are included. One can define +" a threshold energy separating inner from outer shells. By default this +" energy is set to 1 keV, but for accurate calculation of quantities that " +" require knowledge of which particle deposited the energy, one might need to +" use the a lower threshold. +"***************************************************************************" +"============================================================" +REPLACE {$RELAX-CUTOFF} WITH {0.001"threshold energy for outer shells"} +REPLACE {$MXPESHELL} WITH {16} "K,L1..L3,M1..M5,N1..N7 + outer shell" +REPLACE {$MXNE} WITH {500} "number of energy points per shell " +"============================================================" +REPLACE {;COMIN/PE-SHELL-DATA/;} WITH {; + + common/pe_shell_data/ + pe_xsection($MXNE,$MXELEMENT,0:$MXPESHELL), "shellwise cross sections" + pe_elem_prob($MXNE,$MXELEMENT,$MXMED), "prob. of interaction with an" + "element of a medium" + pe_energy($MXNE,$MXELEMENT), "energy grid" + pe_zsorted($MXELEMENT,$MXMED),"sorted array of Z for each medium" + pe_be($MXELEMENT,$MXPESHELL), "binding energies" + pe_nshell($MXELEMENT), "number of shells for each element" + pe_zpos($MXELEMENT), "position of each Z element" + pe_nge($MXELEMENT), "number of energy points for each element" + pe_ne; "number of elements in the simulation" + $REAL pe_be, pe_energy, pe_xsection, pe_elem_prob; + $INTEGER pe_zsorted, pe_nshell, pe_zpos, pe_nge, pe_ne; }; ; "---------- BUFFER FLUSH SEMICOLON ----------" +REPLACE {$COMIN-SHELLWISE-PE-INIT;} WITH { +;COMIN/BREMPR,EDGE,EGS-IO,MEDIA,PHOTIN,THRESH,X-OPTIONS,USEFUL, + PE-SHELL-DATA/; +}; + + +; "---------- BUFFER FLUSH SEMICOLON ----------" " Some macros for C-style syntax in mortran " " Unfortunately, /= doesn't work because of /v1,v2,...,vn/=value;" @@ -1095,7 +1140,12 @@ REPLACE {;COMIN/USER/;} WITH { ;} "DEFAULT IS NULL" ; "---------- BUFFER FLUSH SEMICOLON ----------" -REPLACE {;COMIN/X-OPTIONS/;} WITH {;}; +REPLACE {;COMIN/X-OPTIONS/;} WITH { + ; + common/x_options/eadl_relax, "Use EADL relaxation" + mcdf_pe_xsections;"Use Sabbatucci and Salvat PE xsections" + $LOGICAL eadl_relax, mcdf_pe_xsections; +}; "------------------------------------------------------------------" "*** MACROS TO DEFINE THE COMMONS USED IN EACH SUBPROGRAM. " @@ -1119,7 +1169,7 @@ UPHIIN,UPHIOT,USEFUL,USER,RANDOM,ET-Control,CH-Steps,EGS-IO, EGS-VARIANCE-REDUCTION,EMF-INPUTS/;} REPLACE {$COMIN-HATCH;} WITH { ;COMIN/DEBUG,BOUNDS,BREMPR,ELECIN,MEDIA,MISC,PHOTIN,STACK,THRESH, -UPHIIN,UPHIOT,USEFUL,USER,RANDOM,ET-Control,EGS-IO, +UPHIIN,UPHIOT,USEFUL,USER,RANDOM,ET-Control,EGS-IO,X-OPTIONS, EGS-VARIANCE-REDUCTION/;} REPLACE {$COMIN-MOLLER;} WITH { ;COMIN/DEBUG,STACK,THRESH,UPHIOT,USEFUL,RANDOM,EGS-IO, @@ -1129,7 +1179,7 @@ REPLACE {$COMIN-PAIR;} WITH { EPCONT,TRIPLET-DATA,EGS-VARIANCE-REDUCTION,EGS-IO/;} REPLACE {$COMIN-PHOTO;} WITH { ;COMIN/BOUNDS,BREMPR,DEBUG,EDGE,EPCONT,MEDIA,PHOTIN,RANDOM, - STACK,UPHIOT,USEFUL,EGS-IO, + STACK,UPHIOT,USEFUL,EGS-IO,X-OPTIONS, EGS-VARIANCE-REDUCTION,RELAX-DATA/;} REPLACE {$COMIN-PHOTON;} WITH { ;COMIN/DEBUG,BOUNDS,MEDIA,MISC,EPCONT,PHOTIN,STACK,THRESH,UPHIOT, @@ -1143,21 +1193,21 @@ REPLACE {$COMIN-BLOCK;} WITH { EPCONT,CH-Steps,ET-Control,MEDIA,MISC,PHOTIN,RANDOM,STACK, THRESH, UPHIIN,UPHIOT,USEFUL,EGS-VARIANCE-REDUCTION/;} REPLACE {$COMIN-RELAX;} WITH { - ;COMIN/EPCONT,STACK,BOUNDS,USEFUL,RANDOM,EDGE,EGS-IO,RELAX-DATA/;} + ;COMIN/EPCONT,STACK,BOUNDS,USEFUL,RANDOM,EDGE,EGS-IO,RELAX-DATA,X-OPTIONS/;} REPLACE {$COMIN-SET-DEFAULTS;} WITH { ;COMIN/BOUNDS,BREMPR,COMPTON-DATA,EDGE,ELECIN,EPCONT,CH-Steps,ET-Control, MEDIA,MISC,PHOTIN,RANDOM,STACK,THRESH, UPHIIN,UPHIOT,USEFUL, EGS-VARIANCE-REDUCTION,EGS-IO,Spin-Data,EII-DATA,rayleigh_inputs, - EMF-INPUTS/;}; + EMF-INPUTS,X-OPTIONS/;}; REPLACE {$COMIN-INIT-COMPT;} WITH { - ;COMIN/COMPTON-DATA,BREMPR,EDGE,MEDIA,MISC,USEFUL,EGS-IO/;}; + ;COMIN/COMPTON-DATA,BREMPR,EDGE,MEDIA,MISC,USEFUL,EGS-IO,X-OPTIONS/;}; REPLACE {$COMIN-MSCATI;} WITH { ;COMIN/BOUNDS,ELECIN,MEDIA,MISC,RANDOM,ET-Control,USEFUL,EGS-IO/;}; REPLACE {$COMIN-INIT-TRIPLET;} WITH { ;COMIN/BREMPR,EGS-IO,MEDIA,TRIPLET-DATA,USEFUL/;}; REPLACE {$COMIN-GET-TRANSPORTP;} WITH { ;COMIN/GetInput,BOUNDS,ET-Control,EDGE,COMPTON-DATA,MEDIA,MISC, - BREMPR,EII-DATA,EGS-IO,rayleigh_inputs,EMF-INPUTS/;}; + BREMPR,EII-DATA,EGS-IO,rayleigh_inputs,EMF-INPUTS,X-OPTIONS/;}; "Ali:photonuc, 1 block" REPLACE {$COMIN-PHOTONUC;} WITH {;COMIN/STACK,EPCONT,USEFUL/;}; @@ -2236,7 +2286,7 @@ REPLACE {$EXACT-BCA-DEFAULT} WITH {.true.} ; REPLACE {$SPIN-EFFECTS-DEFAULT} WITH {.true.} ; -REPLACE {$IRAYLR-DEFAULT} WITH {0} +REPLACE {$IRAYLR-DEFAULT} WITH {1} ; REPLACE {$AP-DEFAULT} WITH {-1} ; @@ -2248,6 +2298,12 @@ REPLACE {$XDATA-DEFAULT} WITH {'xcom'} ; REPLACE {$COMP-XDATA-DEFAULT} WITH {'default'} ; +"EADL relaxation is now the default" +REPLACE {$EADL-RELAX-DEFAULT} WITH {.true.} +; +"Sabbatucci and Salvat PE xsections not the default yet" +REPLACE {$MCDF-PE-DEFAULT} WITH {.false.} +; "Ali:photonuc, 2 lines" REPLACE {$IPHOTONUCR-DEFAULT} WITH {0} ; diff --git a/HEN_HOUSE/src/egsnrc.mortran b/HEN_HOUSE/src/egsnrc.mortran index b7b6aa4c0..8847a185d 100644 --- a/HEN_HOUSE/src/egsnrc.mortran +++ b/HEN_HOUSE/src/egsnrc.mortran @@ -2239,7 +2239,7 @@ call mscati; "Initialize new MS, step-sizes, etc, IK Oct 97" "Calling order of the subroutines below is important when using" "detailed atomic relaxation in order to use the binding energies" "corresponding to the requested photon cross section library" -IF ($EADL_RELAX & photon_xsections = 'xcom' )[ +IF ( eadl_relax & photon_xsections = 'xcom' )[ call init_compton; "Initialize bound Compton scattering" call EDGSET(1,1); "Initialize relaxations and photo-absorption data" ] @@ -2251,6 +2251,10 @@ ELSE[ ] +IF( xsec_out = 1 ) [ + call egs_print_binding_energies; +] + call fix_brems; "Re-calculate dl1,... for the different technique" "employed in BREMS. Note that the old EGS sampling" "technique for BREMS had a bug that shows up only" @@ -2271,11 +2275,11 @@ call init_triplet; " SETUP IS NOW COMPLETE" IF (NMED.EQ.1)[ - $egs_info(*,' EGSnrc SUCCESSFULLY ''HATCHED'' FOR ONE MEDIUM.'); + $egs_info(*,'EGSnrc SUCCESSFULLY ''HATCHED'' FOR ONE MEDIUM.'); ] ELSE[ $egs_info('(a,i5,a)', - ' EGSnrc SUCCESSFULLY ''HATCHED'' FOR ',NMED,' MEDIA.'); + 'EGSnrc SUCCESSFULLY ''HATCHED'' FOR ',NMED,' MEDIA.'); ] RETURN; @@ -2438,16 +2442,17 @@ $RADC_HATCH; $need_bound_compton_data(getd); IF( ~getd ) [ - IF($EADL_RELAX)[ - $egs_fatal('(a,/a)', - 'You must turn ON Compton binding corrections when using a', - 'detailed atomic relaxation ($EADL_RELAX=true)!'); + IF( eadl_relax & photon_xsections = 'xcom' )[ + $egs_fatal('(a,/a,/a)', + 'You must turn ON Compton binding corrections when using', + 'a detailed atomic relaxation (eadl_relax=true) since ', + 'binding energies taken from incoh.data below 1 keV!'); ] $egs_info('(a/)',' Bound Compton scattering not requested! '); return; ] -$egs_info('(a$)',' Bound Compton scattering requested, reading data ......'); +$egs_info('(/a$)','Bound Compton scattering requested, reading data ......'); rewind($INCOHUNIT); DO j=1,18 [ read($INCOHUNIT,*); ] "skip 1st 18 lines of comments" iz = 0; @@ -2463,7 +2468,7 @@ DO j=1,$MXTOTSH [ "For detailed atomic relaxations set shell type "to actual shell number and update binding energies "with values from the photo-electric cross sections - IF ($EADL_RELAX)[ + IF (eadl_relax)[ IF (iz_array(j) ~= iz)[ shn_array(j) = 1; iz = iz_array(j); ] @@ -2538,19 +2543,30 @@ DO medium = 1,nmed [ ] ] -$egs_info('(a/)','...... Done.'); +$egs_info('(a/)',' ...... Done.'); $need_relaxation_data(getd); IF( getd ) return; +$egs_fatal('(/a,/a,/a,/a)', +' In subroutine init_compton: ', +' Scattering off bound electrons creates atomic vacancies,', +' potentially starting an atomic relaxation cascade. ', +' Please turn ON atomic relaxations.'); +/* +Turning ON relaxations to setup relaxations for bound Compton +and then turning it back OFF seems inconsistent. One should have +relaxations for all interactions with atomic electrons. + $egs_info('(a/,a/,a/,a//)', ' In subroutine init_compton: ', ' fluorescence not set but relaxation data are required for ', ' bound Compton scattering. ', ' calling EDGSET. '); iedgfl(1) = 1; "This was (2) originally DR" +eadl_relax = .true.; call edgset(1,1); iedgfl(1) = 0; "This was (2) originally DR" - +*/ return; end; @@ -3241,7 +3257,7 @@ $declare_write_buffer; integer*4 i,j,k; -$egs_info('(a,$)',' Reading screened Rutherford MS data ............... '); +$egs_info('(/a,$)','Reading screened Rutherford MS data ............... '); rewind($MSCAT-DATAFILE); DO i=0,$MAXL_MS [ DO j=0,$MAXQ_MS [ @@ -5025,6 +5041,8 @@ $DEFINE-LOCAL-VARIABLES-PHOTO; /* $REAL ftot,iprob; */ data n_warning/0/; +IF ( mcdf_pe_xsections )[call egs_shellwise_photo();return;] + NPold = NP; "Set the old stack counter" PEIG=E(NP); irl = ir(np); IF( peig < edge_energies(2,1) ) [ @@ -5093,15 +5111,9 @@ IF( iedgfl(irl) ~= 0 ) [ " User requested atomic relaxations " " left for now as before, to be changed!!! " IF( peig <= binding_energies($MXSHELL,iZ) ) [ "Outer shells, no atomic relaxation" - IF ($EADL_RELAX)[ "EADL relax: Below M2-shell -> just emit e- " iq(np) = -1; e(np) = peig + prm; - ] - ELSE["Default: Below N-shell -> local energy deposition " - edep = peig; - e(np) = pzero; "wt(np) = 0;" - ] ] ELSE ["Above N-shell -> sample the shell the photon is interacting with" $RANDOMSET br; /* ftot = 1; */ @@ -5117,7 +5129,7 @@ IF( iedgfl(irl) ~= 0 ) [ " User requested atomic relaxations " "EADL APPROACH 1: Do not allow interaction below L3. Deviates" "**************** from previous EGSnrc approach as it doesn't" " generate e- nor x-rays from and shells." - IF ($EADL_RELAX & k > 4)[ + IF (eadl_relax & k > 4)[ "No initial vacancy below L3 for now, just emit e-" iq(np) = -1; e(np) = peig + prm; @@ -5157,6 +5169,365 @@ $PLAY RUSSIAN ROULETTE WITH ELECTRONS FROM NPold; " TO NP;" return; end; +%E +"******************************************************************" +subroutine egs_shellwise_photo; +"******************************************************************" +" Derived from PHOTO by I. Kawrakow and A.F. Bielajew " +" Shellwise implementation and " +" sampling optimizations " +"******************************************************************" + +$IMPLICIT-NONE; + +$COMIN-PHOTO; "default replacement is: + "COMIN/BOUNDS,DEBUG,EDGE,EGS-VARIANCE-REDUCTION,EPCONT," + "MEDIA,PHOTIN,RANDOM,STACK,UPHIOT,USEFUL/" +;COMIN/PE-SHELL-DATA/; + +$DEFINE-VARIABLES-FOR-SELECT-PHOTOELECTRON-DIRECTION; +$DEFINE-LOCAL-VARIABLES-PHOTO; +$REAL slope, logE, int_prob; +$INTEGER zpos, ibsearch; +data n_warning/0/; + +NPold = NP; "Set the old stack counter" +PEIG=E(NP); irl = ir(np); +do_relax = .false.; +IF( peig < $RELAX-CUTOFF ) [ + IF( n_warning < 100 ) [ + n_warning = n_warning + 1; + $egs_info(*,' Subroutine egs_shellwise_photo called with E = ', + peig,' which is below the current min. energy of ', + $RELAX-CUTOFF,' keV! '); + $egs_info(*,' Converting now this photon to an electron, '); + $egs_info(*,' but you should check your code! '); + ] + iq(np) = -1; + e(np) = peig + prm; + return; +] + +edep = pzero; + +IF( iedgfl(irl) ~= 0 ) [" User requested atomic relaxations " + " sample element and atomic shell for" + j = -1; " the interaction." + IF( nne(medium) = 1 ) [ + iZ = int( zelem(medium,1) + 0.5 ); zpos = pe_zpos(iZ); + IF( pe_nshell(zpos) > 0) [ + logE = log(peig); + j = ibsearch(logE,pe_nge(zpos),pe_energy(1,zpos)); + ] + ] + ELSE [ + $RANDOMSET br; logE = log(peig); + "DO k=1,nne(medium) [" + DO k=nne(medium),1,-1 [ + iZ = int( zelem(medium,k) + 0.5 );zpos = pe_zpos(iZ); + IF( iZ < 1 | iZ > $MXELEMENT ) [ + $egs_info(*,' Error in egs_shellwise_photo: '); + $egs_fatal(*,' Atomic number of element ',k, + ' in medium ',medium,' is not between 1 and ',$MXELEMENT); + ] + j = ibsearch(logE,pe_nge(zpos),pe_energy(1,zpos)); + slope = pe_elem_prob(j+1,k,medium) - pe_elem_prob(j,k,medium); + slope = slope/(pe_energy(j+1,zpos)-pe_energy(j,zpos)); + int_prob = pe_elem_prob(j,k,medium)+slope*(logE-pe_energy(j,zpos)); + br -= exp(int_prob); + IF ( br <= 0 ) EXIT; + ] + ] + " Now we know the atomic number (iZ) and the energy interval the " + " photon energy is in (j). It is time to sample the shell the photon " + " is interacting with. " + IF( peig < pe_be(zpos,pe_nshell(zpos)) | pe_nshell(zpos) = 0 ) + [ "no atomic relaxation, create photo-electron" + iq(np) = -1; + e(np) = peig + prm; + ] + ELSE ["sample the shell the photon is interacting with" + $RANDOMSET br; sigtot = 0; + DO k=1,pe_nshell(zpos) [ + IF( peig > pe_be(zpos,k) ) [ + slope = pe_xsection(j+1,zpos,k) - pe_xsection(j,zpos,k); + slope = slope/(pe_energy(j+1,zpos)-pe_energy(j,zpos)); + int_prob=pe_xsection(j,zpos,k)+slope*(logE-pe_energy(j,zpos)); + br -= exp(int_prob); sigtot += exp(int_prob); + IF ( br <= 0 ) EXIT; + ] + ] + IF (k > pe_nshell(zpos))["outer shell, create photo-electron" + iq(np) = -1; + e(np) = peig + prm; + ] + ELSE[ + e_vac = pe_be(zpos,k); + e(np) = peig - e_vac + prm; do_relax = .true.; + iq(np) = -1; + ] + ] +] +ELSE ["No atomic relaxations, just create photo-electron" + e(np) = peig + prm; iq(np) = -1; +] + +IF( iq(np) = -1 ) [ + $SELECT-PHOTOELECTRON-DIRECTION; "Samples photo-electron direction" +] + +IF( do_relax ) [ + call egs_eadl_relax(iZ,k); +] + +IF( EDEP > 0 ) [ $AUSCALL($PHOTXAUS); ] "generates IARG = 4 call" + +; +$PLAY RUSSIAN ROULETTE WITH ELECTRONS FROM NPold; " TO NP;" + +return; +end; + +"*************************************************************************" +subroutine egs_read_shellwise_pe; +"*************************************************************************" + +$IMPLICIT-NONE; + +$declare_max_medium; +$COMIN-SHELLWISE-PE-INIT; + +$INTEGER lnblnk1,egs_get_unit,pe_sw_unit,ounit,egs_open_file; +$INTEGER sorted($MXELEMENT),i,j,k,l,m; +$REAL z_sorted($MXELEMENT),pz_sorted($MXELEMENT); +$REAL rest_xs($MXNE,$MXELEMENT); +$REAL tmp_e($MXNE,$MXPESHELL), tmp_xs($MXNE,$MXPESHELL); +$REAL new_e($MXNE),deltaEb,slope; +$INTEGER zread($MXELEMENT),ib($MXPESHELL),ibsearch; +character data_dir*128,pe_sw_file*144; + +$INTEGER medio,iZ,iZpos,egs_read_int,pos,curr_rec; +real*4 egs_read_real,e_r, e_old,sigma_r; +integer*2 nz, egs_read_short,ish, i_nshell,i_nge; +$LOGICAL is_open, is_there, shift_required; + +character*3 labels(16); +data labels/' K',' L1',' L2',' L3', + ' M1',' M2',' M3',' M4',' M5', + ' N1',' N2',' N3',' N4',' N5',' N6',' N7'/; + +/********************************************* + Open PE shellwise data file photo_shellwise.data + **********************************************/ +$egs_info('(/a$)', +' Reading renormalized photoelectric cross sections ......'); +data_dir = $cstring(hen_house) // 'data' // $file_sep; +pe_sw_file = $cstring(data_dir) // 'photo_shellwise.data'; + +/* Open a fortran unit for reading data */ +pe_sw_unit = egs_get_unit(0); +IF( pe_sw_unit < 1 ) [ + $egs_fatal(*,'egs_init_shellwise_pe: failed to get a free Fortran I/O unit'); +] +open(pe_sw_unit,file=pe_sw_file,status='old', + form='UNFORMATTED',ACCESS='direct',recl=1, + err=:no-pe-sw-file:); +GOTO :read-pe-sw:; +:no-pe-sw-file: +$egs_fatal('(2a)','egs_init_shellwise_pe: failed to open ', + pe_sw_file); +:read-pe-sw: +is_open = .true.; +/****************************** + Array initialization + ******************************/ +DO medio = 1,nmed [ + DO i=1,nne(medio) [ + pe_nshell(i*medio) = 0; + pe_nge(i*medio) = 0; + pe_zsorted(i,medio) = 0; + ] +] +DO l = 1,$MXELEMENT [ + pe_zpos(l) = -1; + DO k = 1,$MXNE [ + pe_energy(k,l) = 0.0; + DO m = 1,$MXPESHELL [ + pe_xsection(k,l,m) = 0.0; + ] + ] + DO k = 1,$MXPESHELL [ + pe_be(l,k) = -99; + ] +] + +/****************************** + Get shellwise PE xsections + ******************************/ +curr_rec = 1; iZpos = 0; +nz = egs_read_short(pe_sw_unit,curr_rec); +"$egs_info('(a,i2,a)','PE shellwise data available for ',nz,' elements....');" +DO medio = 1,nmed [ + DO i=1,nne(medio) [z_sorted(i) = zelem(medio,i);] + call egs_heap_sort(nne(medio),z_sorted,sorted); + DO i=1,nne(medio) [pe_zsorted(i,medio) = z_sorted(i);] + DO i=1,nne(medio) [ + iZ = z_sorted(i); + "Now check whether we have already loaded the data for" + "this atomic number" + is_there = .false.; + DO j = 1,medio-1 [ + DO k = 1, nne(j)[ + IF( iZ = pe_zsorted(k,j) ) [ + is_there = .true.; EXIT; + ] + ] + ] + IF (is_there) NEXT; + "Read data for this element" + iZpos += 1; zread(iZpos) = iZ; + pe_zpos(iZ) = iZpos; + pos = 3 + (iZ-1)*4; + curr_rec = egs_read_int(pe_sw_unit,pos) + 1; + i_nge = egs_read_short(pe_sw_unit,curr_rec); + i_nshell = egs_read_short(pe_sw_unit,curr_rec); + "$egs_info('(a,i2,a,i3,a,i2,a)','Element ',iZ,' has ',i_nge, + " ' energy points and ',i_nshell,' shells'); + pe_nge(iZpos) = i_nge; pe_nshell(iZpos) = i_nshell; + e_old = -1.0; ish = 0; + DO j = 1,i_nge[ + e_r = egs_read_real(pe_sw_unit,curr_rec); + sigma_r = egs_read_real(pe_sw_unit,curr_rec); + pe_energy(j,iZpos) = e_r; + pe_xsection(j,iZpos,0) = sigma_r; + rest_xs(j,iZpos) = sigma_r; + DO k = 1, i_nshell[ + sigma_r = egs_read_real(pe_sw_unit,curr_rec); + pe_xsection(j,iZpos,k) = sigma_r; + rest_xs(j,iZpos) = rest_xs(j,iZpos) - sigma_r; + + ] + "Extract binding energies from the data base" + IF (e_r - e_old < 1e-15)[ + pe_be(iZpos,i_nshell-ish) = e_r; + ish += 1; + ] + e_old = e_r; + + ] + "DO k = 1, i_nshell[ + " $egs_info('(a2,a3,a2,1pe12.4,a4)', + " 'E(',labels(k),')=',pe_be(iZpos,k),' MeV'); + "] + + ] +] +pe_ne = iZpos; +"$egs_info('(a,i2,a/)','Finished processing ',iZpos,' elements!'); + +/********************************************* + Adjust xsections to current binding energies + Shift energy scale for each subshell when + required. + *********************************************/ +"$egs_info(*,' Adjusting cross sections to new binding energies ...'); +DO i=1,pe_ne[ + iZ = zread(i); + IF (pe_nshell(i) = 0)[ + DO j=1,pe_nge(i)[ + pe_energy(j,i) = log(pe_energy(j,i)); + ] + NEXT; + ] + "Shift energy scale for different binding energy sets" + DO l=1,pe_nshell(i)[ + IF ( pe_be(i,l) ~= binding_energies(l,iZ))[ + shift_required = .true.; + deltaEb = binding_energies(l,iZ)-pe_be(i,l); + "$egs_info('(2(a,a,a,1pe12.4),a,1pe12.4,a,e12.4)', + " 'Eb_p_',labels(l),' = ',pe_be(i,l), + " ' Eb_e_',labels(l),' = ',binding_energies(l,iZ), + " ' diff = ',pe_be(i,l)-binding_energies(l,iZ), + " ' -> ', 100*(1.0 - binding_energies(l,iZ) / pe_be(i,l)) + " ); + ] + ELSE[shift_required =.false.;] + is_there = .false.; + DO j=1,pe_nge(i)[ + tmp_e(j,l) = pe_energy(j,i); + tmp_xs(j,l) = pe_xsection(j,i,l); + IF ( shift_required & + pe_energy(j,i) => pe_be(i,l) )[ + tmp_e(j,l) += deltaEb; + "$egs_info(*,'Shifting ',pe_energy(j,i),' to ',tmp_e(j,l)); + "Determine edge position in energy array" + IF (pe_energy(j,i) = pe_be(i,l) & ~is_there)[ + ib(l) = j; is_there = .true.; + ] + "Update new energy grid" + IF (l = 1)[ + new_e(j) = tmp_e(j,l); + ] + "ELSE IF(tmp_e(j,l) < binding_energies(l-1,iZ))[ + ELSE IF(j < ib(l-1))[ + new_e(j) = tmp_e(j,l); + ] + ] + ] + pe_be(i,l) = binding_energies(l,iZ); + ] + "Re-compute sub-shell xsections for new energy grid new_e" + "Not needed for K shell" + DO l=2,pe_nshell(i)[ + DO j=1,pe_nge(i)[ + IF ( new_e(j) >= pe_be(i,l-1) )[ + m = ibsearch(new_e(j),pe_nge(i),tmp_e(1,l)); + slope = log(tmp_xs(m+1,l)/tmp_xs(m,l)); + slope = slope/log(tmp_e(m+1,l)/tmp_e(m,l)); + pe_xsection(j,i,l) = log(tmp_xs(m,l)); + pe_xsection(j,i,l) += slope*log(new_e(j)/tmp_e(m,l)); + pe_xsection(j,i,l) = exp(pe_xsection(j,i,l)); + ] + ] + ] + "Re-compute total xsections for new energy grid new_e" + "$egs_info(*,'-> Z = ',iZ);" + DO j=1,pe_nge(i)[ + IF ( j < ib(pe_nshell(i)))[ + new_e(j) = pe_energy(j,i); + ] + m = ibsearch(new_e(j),pe_nge(i),pe_energy(1,i)); + slope = log(rest_xs(m+1,i)/rest_xs(m,i)); + slope = slope/log(pe_energy(m+1,i)/pe_energy(m,i)); + pe_xsection(j,i,0) = log(rest_xs(m,i)); + pe_xsection(j,i,0) += slope*log(new_e(j)/pe_energy(m,i)); + pe_xsection(j,i,0) = exp(pe_xsection(j,i,0)); + "$egs_info('(1pe12.4,1x,1pe12.4,1x,1pe12.4,1x,1pe12.4)', + " new_e(j),pe_xsection(j,i,0),pe_energy(j,i),rest_xs(j,i)); + DO l=1,pe_nshell(i)[ + pe_xsection(j,i,0) += pe_xsection(j,i,l); + ] + "$egs_info('(1pe12.4,1x,1pe12.4,1x,1pe12.4)', + " new_e(j),pe_xsection(j,i,0)); + ] + "Normalize shell cross sections to total for sampling" + "and update energy grid of ith element." + " Prepare for log/log interpolation." + DO j=1,pe_nge(i)[ + pe_energy(j,i) = log(new_e(j)); + DO l=1,pe_nshell(i)[ + pe_xsection(j,i,l) = log(pe_xsection(j,i,l)/pe_xsection(j,i,0)); + ] + ] +] + +$egs_info('(a/)',' done'); + +IF( is_open ) close(pe_sw_unit); +return; +end; + + "******************************************************************" SUBROUTINE RELAX(energy,n,iZ); "******************************************************************" @@ -5172,9 +5543,6 @@ SUBROUTINE RELAX(energy,n,iZ); " " " Version 1: I. Kawrakow, December 1998 " "******************************************************************" - -REPLACE {$RELAX-CUTOFF} WITH {0.001} - implicit none; " Input variables " @@ -5223,8 +5591,6 @@ $REAL e_array($MXVAC), "array with vacancy energies " $REAL xphi,yphi,xphi2,yphi2,rhophi2, cphi,sphi; "for azimuthal angle selection" -$LOGICAL eadl_relax; - " Global EGS4 variables " "=======================" $COMIN-RELAX; @@ -5251,7 +5617,6 @@ data final_state/ "See the final_state definition above" save first_transition,last_transition,final_state; "to avoid problems with " "non-static compiler options" -eadl_relax=$EADL_RELAX; IF (eadl_relax)[ call egs_eadl_relax(iZ,n); return; @@ -5476,13 +5841,14 @@ DO iZ=1,$MXELEMENT[ /********************************************* Determine minimum energy **********************************************/ -Pc = 1e30; Ec = 1e30; +/*Pc = 1e30; Ec = 1e30; DO medio = 1,nmed [ tmp = AE(medio) - rm; Ec = min(Ec,tmp); Pc = min(Pc,ap(medio)); ] min_be = min(Ec,Pc); "This is the minimum binding energy for which we need" "relaxation data" - +*/ +min_be = $RELAX-CUTOFF; $egs_debug('(a)',' ************ relax_init **************** '); $egs_debug('(a,f10.7)', ' Minimum binding energy requiring relaxation data: ',min_be); @@ -5490,7 +5856,7 @@ $egs_debug('(a,f10.7)', /********************************************* Open EADL relaxation data file relax.data **********************************************/ -$egs_info('(/a)','Reading EADL relaxation data ......'); +$egs_info('(/a)',' Reading EADL relaxation data ......'); data_dir = $cstring(hen_house) // 'data' // $file_sep; relax_file = $cstring(data_dir) // 'relax.data'; @@ -5536,7 +5902,7 @@ DO medio = 1,nmed [ ' Increase the parameter $MAXSHELL and retry '); ] $egs_info('(a,i3,a,i2,a)', - '-> Element Z = ',iZ,' has ',nshell,' shells'); + ' Z = ',iZ,' has ',nshell,' shells'); DO ish=shell_ntot+1,shell_ntot+nshell[ curr_rec = curr_rec+1; read(relax_unit, rec=curr_rec) shell_type(ish); @@ -5560,7 +5926,7 @@ DO medio = 1,nmed [ read(relax_unit,rec=curr_rec) itmp(k); curr_rec = curr_rec+1; read(relax_unit,rec=curr_rec) prob_r;wtmp(k)=prob_r; - /* Relaxation data uses C-style array indexing. + /*Relaxation data uses C-style array indexing. Hence final relaxation state increased by 1 for radiative and by 65 for non-radiative transitions*/ IF (itmp(k)<64) [itmp(k) +=1;] @@ -5602,7 +5968,7 @@ DO medio = 1,nmed [ ] $egs_debug('(a,i4,a)',' There are ', shell_ntot,' shells in the list of shells '); -$egs_info('(a/)','...... Done.'); +$egs_info('(a/)',' ...... Done.'); /***********************************************/ IF( is_open ) close(relax_unit); @@ -5708,13 +6074,13 @@ rfu_j0 = shell; "shell # of vacancy that initiated cascade in long list" rfu_n0 = shell_num(shell); "same but number is shell number in element iZ" rfu_t0 = shell_type(shell); "the shell type (encodes the type of relaxation)" rfu_E0 = Evac; "B.E. of vacancy that initiated cascade" - +" "Local energy deposition for vacancies below L3 shell." "Added here for consistency in algorithm for as long as" " and shells considered instead of proper shells." "AUSGAB call should be updated when charge of particle" "creating vacancy becomes available." -IF (shell_egs > 4) [ +IF (shell_egs > 4 & ~mcdf_pe_xsections) [ edep = Evac; "add energy of vacancy to edep" edep_local = Evac; "set value of edep_local to energy of vacancy" @@ -5960,16 +6326,16 @@ IF( got_data ) return; "statement. $egs_info('(a/,a)', - ' Output from subroutine EDGSET:', - ' =============================='); + 'Output from subroutine EDGSET:', + '=============================='); $need_relaxation_data(do_relax); IF( ~do_relax ) [ - IF($EADL_RELAX)[ + IF(eadl_relax)[ $egs_fatal('(a,/a)', 'You must turn ON atomic relaxations when requesting', - 'detailed atomic relaxation ($EADL_RELAX=true)!'); + 'detailed atomic relaxation (eadl_relax=true)!'); ] $egs_info('(a/)',' Atomic relaxations not requested! '); return; @@ -5977,14 +6343,14 @@ IF( ~do_relax ) [ $egs_info('(a/)',' Atomic relaxations requested! '); -$egs_info('(/a$)',' Reading photo-absorption data .....'); +$egs_info('(a$)',' Reading simplified photo-absorption data .....'); got_data = .true.; rewind($PHOTOUNIT); DO i=1,$MXELEMENT [ - IF ($EADL_RELAX)[ + IF (eadl_relax)[ "Skip, using binding_energies from *_photo.data file read($PHOTOUNIT,*); ] @@ -6003,7 +6369,7 @@ DO i=1,$MXELEMENT ] $egs_info('(a)',' Done'); -$egs_info('(/a$)',' Reading relaxation data .....'); +$egs_info('(/a$)',' Reading simplified relaxation data .....'); read($PHOTOUNIT,*); DO i=1,$MXELEMENT [ read($PHOTOUNIT,*) j,(relaxation_prob(k,i),k=1,19); "K-shell" @@ -6025,7 +6391,7 @@ DO i=1,$MXELEMENT [ read($PHOTOUNIT,*) j,relaxation_prob(38,i); "M-shell" ] $egs_info('(a)',' Done'); -$egs_info('(/a$)',' Reading photo cross section data .....'); +$egs_info('(/a$)',' Reading parametrized XCOM photo cross section data .....'); rewind($PHOCSUNIT); DO i=1,$MXELEMENT [ read($PHOCSUNIT,*) j,edge_number(i); @@ -6034,9 +6400,9 @@ DO i=1,$MXELEMENT [ edge_d(j,i),edge_energies(j,i); ] ] -$egs_info('(a/)',' Done'); +$egs_info('(a)',' Done'); -IF ($EADL_RELAX)[ +IF (eadl_relax)[ call egs_init_relax; ] @@ -6610,10 +6976,10 @@ call gauss_legendre(0d0,1d0,x_gauss,w_gauss,ngauss); $egs_info(*,' '); IF (ibr_nist = 1) [ -$egs_info(*,' Using NIST brems cross sections! '); +$egs_info(*,'Using NIST brems cross sections! '); ] ELSE IF (ibr_nist = 2) [ - $egs_info(*,' Using NRC brems cross sections! '); + $egs_info(*,'Using NRC brems cross sections! '); ] $egs_info(*,' '); DO medium=1,nmed [ @@ -7093,23 +7459,14 @@ DO j=1,$MXELEMENT [ eii_nshells(j) = 0; ] DO j=1,$MXMED [ eii_nsh(j) = 0; ] IF( eii_flag = 0 ) [ return; ] -/***************************************************************** - Get relaxation data if using PEGS4 data set and eii_flag>0, - so that the binding energies are available in eii_init. -*****************************************************************/ - if (toUpper($cstring(photon_xsections)) = 'PEGS4')[ - $need_relaxation_data(getd); - IF( ~getd )[ - $egs_info('(a/,a/,a/,a//)', - ' In subroutine eii_init: ', - ' Fluorescence not set, but binding energies required for ', - ' EII when using PEGS4 data. ', - ' Verifying if binding energies available by calling EDGSET. '); - itmp = iedgfl(1);iedgfl(1) = 1; - call edgset(1,1); iedgfl(1) = itmp; - ] - ] -/*****************************************************************/ +$need_relaxation_data(getd); +IF( ~getd )[ + $egs_fatal('(/a,/a,/a,/a)', + ' In subroutine eii_init: ', + ' Scattering off bound electrons creates atomic vacancies,', + ' potentially starting an atomic relaxation cascade. ', + ' Please turn ON atomic relaxations.'); +] /* find minimum of all threshold energies @@ -7160,15 +7517,6 @@ nsh_tot = nsh; tmp_array(1) = 0; DO j=2,$MXELEMENT [ tmp_array(j) = tmp_array(j-1) + eii_nshells(j-1); ] -/* - Get relaxation data if necessary so that the binding energies - are available - */ -itmp = iedgfl(1); -iedgfl(1) = 1; -call edgset(1,1); -iedgfl(1) = itmp; - /* set EII active shells per medium and for each element */ DO imed=1,nmed [ nsh = 0; @@ -7657,16 +8005,15 @@ character data_dir*128,photo_file*140,pair_file*140,rayleigh_file*144, "Ali:photonuc, 1 line" photonuc_file*144; -$egs_info('(/a)','(Re)-initializing photon cross section data'); -$egs_info('(a,a/)','Using data files from the series ', +$egs_info('(/a$)','(Re)-initializing photon cross sections'); +$egs_info('(a,a/)',' with files from the series: ', prefix(:lnblnk1(prefix))); -$egs_info('(a,a,a,a)','Photon cross sections: ',$cstring(prefix), - ' Compton cross sections: ',$cstring(comp_prefix)); +$egs_info('(a,a)',' Compton cross sections: ',$cstring(comp_prefix)); "Ali:photonuc, 1 block" IF(iphotonuc = 1) [ - $egs_info('(a,a)','Photonuclear cross sections: ', + $egs_info('(a,a)',' Photonuclear cross sections: ', $cstring(photonuc_prefix)); input_photonuc_data = .false.; IF(lnblnk1(photonuc_prefix) > 0 & photonuc_prefix(1:7) ~= 'default') [ @@ -7692,7 +8039,7 @@ ELSE [ ] "Ali: I moved this info line from inside the IF statement " because it's useful to print the cross section file either way -$egs_info('(a,a)','Using Compton cross sections from ', +$egs_info('(a,a)',' Using Compton cross sections from ', $cstring(compton_file)); "Ali:photonuc, 1 block" @@ -7704,7 +8051,7 @@ IF(iphotonuc = 1) [ ELSE [ photonuc_file = $cstring(data_dir) // 'iaea_photonuc.data'; ] - $egs_info('(a,a)','Using photonuclear cross sections from ', + $egs_info('(a,a)',' Using photonuclear cross sections from ', $cstring(photonuc_file)); ] @@ -7762,18 +8109,20 @@ DO iz=1,100 [ k,' binding energies read exceeding array size of', $MXSHXSEC,'Increase $MXSHXSEC in egsnrc.macros!'); ] - IF( ~$EADL_RELAX & k >= 4 ) EXIT; + IF( ~eadl_relax & k >= 4 ) EXIT; ] ] ] +IF (mcdf_pe_xsections)[call egs_read_shellwise_pe();] + DO medium = 1,nmed [ mge(medium) = $MXGE; nge = $MXGE; ge1(medium) = nge-1; ge1(medium) = ge1(medium)/log(up(medium)/ap(medium)); ge0(medium) = 1 - ge1(medium)*log(ap(medium)); - $egs_info('(a,i3,a,$)','Working on medium ',medium,' ... '); + $egs_info('(a,i3,a,$)',' Working on medium ',medium,' ... '); IF( out = 1 ) [ write(ounit,'(/,2x,a,i3,a,24a1/)') 'Medium ',medium,': ', (media(k,medium),k=1,24); @@ -7789,8 +8138,14 @@ DO medium = 1,nmed [ call egs_heap_sort(nne(medium),z_sorted,sorted); DO i=1,nne(medium) [ pz_sorted(i) = pz(medium,sorted(i)); ] - call egsi_get_data(0,photo_unit,nge,nne(medium),z_sorted,pz_sorted, + IF (mcdf_pe_xsections)[ + call egsi_get_shell_data(medium,nge,nne(medium),z_sorted,pz_sorted, + ge1(medium),ge0(medium),sig_photo); + ] + ELSE[ + call egsi_get_data(0,photo_unit,nge,nne(medium),z_sorted,pz_sorted, ge1(medium),ge0(medium),sig_photo); + ] call egsi_get_data(0,rayleigh_unit,nge,nne(medium),z_sorted,pz_sorted, ge1(medium),ge0(medium),sig_rayleigh); call egsi_get_data(1,pair_unit,nge,nne(medium),z_sorted,pz_sorted, @@ -8042,7 +8397,7 @@ ELSE[ "To avoid log(0) = -INF, x=1E-4 corresponds ~ 0.014 deg for 10 keV" "and much less for higher energies" IF( xgrid(1,medium) < 1e-6 ) xgrid(1,medium) = 1e-4; - write(*,*) '\n -> ', nff, ' atomic ff values computed!'; + $egs_info('(/a,i4,a)',' -> ', nff, ' atomic ff values computed!'); ] close(ff_unit); /* call routine for preparing sampling data */ @@ -8140,7 +8495,8 @@ nff = $MXRAYFF; "To avoid log(0) = -INF, x=1E-4 corresponds ~ 0.014 deg for 10 keV" "and much less for higher energies" IF( xgrid(1,medium) < 1e-6 ) xgrid(1,medium) = 1e-4; -write(*,*) '\n -> ', nff, ' atomic ff values computed!'; +"write(*,*) '\n -> ', nff, ' atomic ff values computed!'; +$egs_info('(/a,i4,a)',' -> ', nff, ' atomic ff values computed!'); close(ff_unit); /* call routine for preparing sampling data */ @@ -8483,68 +8839,109 @@ DO i=1,ne [ ] ] -/* Ali:photonuc. The whole routine below is commented out and -re-written above to accommodate reading photonuclear cross sections. +return; -rewind(iunit); -iz_old = 0; -DO k=1,n [ data(k) = 0; ] +:user-data-failure:; +$egs_fatal(*,'Error while reading user photon cross sections from unit ', + iunit); + +return; end; + +/***************************************************************** + Prepare photoelectric cross section data base for medium imed + on same energy grid as for the other photon interactions. + Normalize elemental PE cross sections to medium cross section + to use for selecting element the photon interacts with. + *****************************************************************/ +"=============================================================================" +subroutine egsi_get_shell_data(imed,n,ne,zsorted,pz_sorted,ge1,ge0,data); +"==========================================================================" +implicit none; +$declare_max_medium; +;COMIN/EGS-IO,BREMPR,USEFUL,MEDIA,PE-SHELL-DATA/; +$INTEGER n, "number of data points requested" + ne, "number of elements in medium" + ndat;"number of data points from original grid" +$REAL ge1,ge0,zsorted(*),pz_sorted(*),data(*); +$REAL sigma($MXNE),sigmaMedium; +"$INTEGER sorted(*); +real*4 etmp($MXINPUT),ftmp($MXINPUT); +real*4 gle,sig,p; +$INTEGER i,j,k,kk,iz,zpos,imed; + +DO k=1,n [ data(k) = 0;] +DO k=1,ne [ sigma(k) = 0;] DO i=1,ne [ - iiz = int(zsorted(i)+0.5); - DO iz=iz_old+1,iiz [ - read(iunit,*,err=:user-data-failure:) ndat; - IF( ndat > $MXINPUT ) [ - $egs_fatal(*,'Too many input data points. Max. is ',$MXINPUT); - ] - IF( flag = 0 ) [ - read(iunit,*,err=:user-data-failure:) (etmp(k),ftmp(k),k=1,ndat); - ] - ELSE [ - read(iunit,*,err=:user-data-failure:) (etmp(k+1),ftmp(k+1), - k=1,ndat); - IF( flag = 1 ) [ eth = 2*rm; ] ELSE [ eth = 4*rm; ] - ndat = ndat + 1; - DO k=2,ndat [ - ftmp(k) = ftmp(k) - 3*log(1-eth/exp(etmp(k))); - ] - ftmp(1) = ftmp(2); etmp(1) = log(eth); - ] - ] - iz_old = iiz; + iz = int(zsorted(i)+0.5); + zpos = pe_zpos(iz); ndat = pe_nge(zpos); + "Total cross sections for a given element" + "on initial energy grid" + DO k=1,ndat[ + pe_elem_prob(k,i,imed) = pz_sorted(i)*pe_xsection(k,zpos,0); + "etmp(k) = log(pe_energy(k,zpos)); Done in egs_read_shellwise_pe + etmp(k) = pe_energy(k,zpos); + ftmp(k) = log(pe_xsection(k,zpos,0)); + ] + "Total cross sections for a given element" + "on requested energy grid" DO k=1,n [ - gle = (k - ge0)/ge1; e = exp(gle); + gle = (k - ge0)/ge1; IF( gle < etmp(1) | gle >= etmp(ndat) ) [ - IF( flag = 0 ) [ - $egs_fatal(*,'Energy ',exp(gle), - ' is outside the available data range of ', - exp(etmp(1)),exp(etmp(ndat))); - ] - ELSE [ - IF( gle < etmp(1) ) [ sig = 0; ] - ELSE [ sig = exp(ftmp(ndat)); ] - ] + $egs_fatal(*,'egsi_get_shell_data: Energy ',exp(gle), + ' is outside the available data range of ', + exp(etmp(1)),exp(etmp(ndat))); ] ELSE [ + "Find energy interval gle falls in" DO kk=1,ndat-1 [ IF( gle >= etmp(kk) & gle < etmp(kk+1) ) EXIT; ] + "log/log interpolation" p = (gle - etmp(kk))/(etmp(kk+1) - etmp(kk)); sig = exp(p*ftmp(kk+1) + (1-p)*ftmp(kk)); ] - IF( flag ~= 0 & e > eth ) sig = sig*(1-eth/e)**3; data(k) = data(k) + pz_sorted(i)*sig; + "data(k) = data(k) + pz(imed,sorted(i))*sig; ] ] -Ali:photonuclear. End of the commented-out routine */ +"Normalize elemental cross section to medium cross section" +"Prepare for log/log interpolation" +DO i=1,ne [ + iz = int(zsorted(i)+0.5); + zpos = pe_zpos(iz); ndat = pe_nge(zpos); + DO k=1,ndat[ + sig = sigmaMedium(imed,pe_energy(k,zpos)); + pe_elem_prob(k,i,imed) = log(pe_elem_prob(k,i,imed)/sig); + ] +] -return; +return; end; -:user-data-failure:; -$egs_fatal(*,'Error while reading user photon cross sections from unit ', - iunit); +/*************************************************************** + Calculate photoelectric cross section for medium imed using + log/log linear interpolation and binary search to find energy + interval + ***************************************************************/ +$REAL function sigmaMedium(imed, logE); +implicit none; +$declare_max_medium; +;COMIN/BREMPR,PE-SHELL-DATA/; +$REAL logE, slope, sigma; +$INTEGER k,imed,Z,zpos,m,ibsearch; + +sigmaMedium = 0; +DO k=1,nne(imed) [ + Z = int( zelem(imed,k) + 0.5 );zpos = pe_zpos(Z); + m = ibsearch(logE,pe_nge(zpos),pe_energy(1,zpos)); + slope = log(pe_xsection(m+1,zpos,0)/pe_xsection(m,zpos,0)); + slope = slope/(pe_energy(m+1,zpos)-pe_energy(m,zpos)); + sigma = log(pe_xsection(m,zpos,0)); + sigma += slope*(logE - pe_energy(m,zpos)); + sigma = exp(sigma); + sigmaMedium += pz(imed,k)*sigma; +] return; end; - "=============================================================================" subroutine egs_heap_sort(n,rarray,jarray); "************************************************************************ diff --git a/HEN_HOUSE/src/get_inputs.mortran b/HEN_HOUSE/src/get_inputs.mortran index 138abf7ab..2f89b3909 100644 --- a/HEN_HOUSE/src/get_inputs.mortran +++ b/HEN_HOUSE/src/get_inputs.mortran @@ -1182,7 +1182,8 @@ $LOGICAL ecut_inregions,pcut_inregions,smax_inregions, pe_inregions,aux_inregions,photonuc_inregions;"Ali:photonuc" character*15 output_strings(14);"Ali:photonuc, increased by 1" -save output_strings,line; +character*16 p_xsections; +save output_strings,line, p_xsections; save ecut_inregions,pcut_inregions,smax_inregions, incoh_inregions,coh_inregions,relax_inregions, pe_inregions,aux_inregions,photonuc_inregions, @@ -1286,6 +1287,8 @@ allowed_inputs(ival,0) = 'Off'; allowed_inputs(ival,1) = 'On'; allowed_inputs(ival,2) = $ON_IN_REGIONS; allowed_inputs(ival,3) = $OFF_IN_REGIONS; +allowed_inputs(ival,4) = 'eadl'; +allowed_inputs(ival,5) = 'simple'; " Photoelectron angular sampling " ival = ival + 1; @@ -1574,10 +1577,24 @@ IF( error_flags(num_eii) = 0 ) [ $cstring(eii_xfile)); ] ] -/**********************************************************/ +/****************************************************** + Use Sabbatucci and Salvat shellwise photoelectric + cross sections together with XCOM or EPDL cross sections + for pair production. + ******************************************************/ IF( error_flags(num_pxsec) = 0 ) [ - photon_xsections = char_value(num_pxsec,1); + p_xsections = char_value(num_pxsec,1); + IF ( toUpper( $cstring(p_xsections) ) = 'MCDF-XCOM' )[ + mcdf_pe_xsections = .true.; photon_xsections = 'xcom'; + ] + ELSEIF ( toUpper( $cstring(p_xsections) ) = 'MCDF-EPDL' )[ + mcdf_pe_xsections = .true.; photon_xsections = 'epdl'; + ] + ELSE[ + mcdf_pe_xsections = .false.; + ] ] +/**********************************************************/ IF( error_flags(num_pxsec_out) = 0 ) [ xsec_out = value(num_pxsec_out,1); ] @@ -1650,8 +1667,14 @@ IF( error_flags(num_coh) = 0 ) [ write(*,'(/)'); ] /***************************************************************/ - - +/* + Inputs allowing setting by region are checked with the macro + $TURN-ON/OFF-IN-REGIONS. For this type of input, the second + and third entries correspond to turning ON or OFF the interaction + in specific regions. If no input is found, the default value is + used. + */ +/***************************************************************/ $TURN-ON/OFF-IN-REGIONS(num_incoh, 'Bound Compton start region', 'Bound Compton stop region', @@ -1690,6 +1713,28 @@ write(i_errors,*) ' ******************** end input transport parameter *********************** '; write(i_errors,*); +"Check if EADL relaxation requested. Note that original relaxation" +"algorithm using and is only turned ON for all regions." +"Moved past the $TURN-ON/OFF-IN-REGIONS statement to catch the" +"default case." +IF(value(num_relax,1) > 0 & value(num_relax,1) < 5)[ + eadl_relax = .true.; + "Default relaxation is EADL" + IF (value(num_relax,1) = 1)[value(num_relax,1)=4;] +] +ELSE [ + IF (mcdf_pe_xsections & value(num_relax,1) = 5)[ + eadl_relax = .true.; value(num_relax,1)=4; + $egs_warning('(a/,a/,a/)', + ' Simplified atomic relaxation not allowed', + ' with shellwise PE cross sections. Resetting', + ' to detailed EADL atomic relaxation!!!'); + ] + ELSE [ + eadl_relax = .false.; + ] +] + " we put the information stored in allowed_inputs into " " output_strings just in case the user over-rides it " " before printing out the settings " @@ -1728,7 +1773,8 @@ write(ounit,'(a,/)') write(ounit,'(a,/)') line; /* initialized in egs_set_defaults */ -write(ounit,'(a,38x,a)') ' Photon cross sections', $cstring(photon_xsections); +write(ounit,'(a,38x,a)') ' Photon cross sections', + $cstring(p_xsections); write(ounit,'(a,37x,a)') ' Compton cross sections', $cstring(comp_xsections); write(ounit,'(a,$)') ' Photon transport cutoff(MeV)'; @@ -1764,7 +1810,6 @@ write(ounit,'(a,30x,a4)') ' Bremsstrahlung cross sections',output_strings(9); write(ounit,'(a,29x,a3)') ' Bremsstrahlung angular sampling',output_strings(6); IF( spin_effects ) [ write(ounit,'(a,48x,a)') ' Spin effects','On'; ] ELSE [ write(ounit,'(a,48x,a)') ' Spin effects','Off'; ] -"write(ounit,'(a,34x,a)') ' Electron Impact Ionization',output_strings(13);" write(ounit,'(a,34x,a)') ' Electron Impact Ionization',$cstring(eii_xfile); IF (eii_L_factor ~= 1.0) [ write(ounit,'(a,25x,f6.4)') @@ -1832,6 +1877,7 @@ $declare_max_medium; $INTEGER imed,ival,lnblnk1,nchanged; character*24 medname; +ounit = i_log; ounit = i_log; delimeter = $THE_DELIMETER; call get_input_set_error_level(0); diff --git a/HEN_HOUSE/src/transportp.macros b/HEN_HOUSE/src/transportp.macros index 3287de390..735336880 100644 --- a/HEN_HOUSE/src/transportp.macros +++ b/HEN_HOUSE/src/transportp.macros @@ -132,8 +132,15 @@ IF( error_flags({P1}) = 0 ) [ DO j=1,$MXREG [ {P4}(j) = itmp; ] ] ] -ELSE [ - value({P1},1) = {P4}(1); +ELSE ["Use default value" + "Check if default is other than ON or Off" + "and shift as no input by regions requested" + IF({P4}(1) = 2 | {P4}(1) = 3)[ + value({P1},1) = {P4}(1)+2; + ] + ELSE[ + value({P1},1) = {P4}(1); + ] ] }; diff --git a/HEN_HOUSE/user_codes/dosxyznrc/dosxyznrc.mortran b/HEN_HOUSE/user_codes/dosxyznrc/dosxyznrc.mortran index 402542c88..4dd3c8f4b 100644 --- a/HEN_HOUSE/user_codes/dosxyznrc/dosxyznrc.mortran +++ b/HEN_HOUSE/user_codes/dosxyznrc/dosxyznrc.mortran @@ -1017,6 +1017,10 @@ REPLACE {$AUSCALL(#);} WITH REPLACE {$IBRDST-DEFAULT} WITH {0} ; +"Rayleigh scattering= Off" +REPLACE {$IRAYLR-DEFAULT} WITH {0} +; + "Bound Compton scattering= Off" REPLACE {$IBCMP-DEFAULT} WITH {0} ; @@ -3416,7 +3420,7 @@ $REAL XSI; "used to perform russian roulette on higher-order photons when" "variables for howfarless option" $INTEGER irx,iry,irz,irold_loc,irnew_loc,outside,min_plane,iri,i,ii,jj,kk, - ibsearch,hwfl_rr,dirx,diry,dirz; + ibsearch4,hwfl_rr,dirx,diry,dirz; $REAL dist_loc, dist_tot, edep_loc, tv,t(2),xi,yi,zi,ui(2),vi(2),wi(2), dist_tmp,x_tmp,y_tmp,z_tmp,edep_save,x_loc,y_loc,z_loc, xpln,ypln,zpln,dnear_loc,t_tot,t_rnd,dist_max,uii,vii,wii, @@ -3468,9 +3472,9 @@ IF(howfarless & (iarg=15 | iarg=17 | iarg=19 | iarg=23) & ir(np)