| Geant4 Cross Reference |
1 /* 1 /*
2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>>
3 # <<END-copyright>> 3 # <<END-copyright>>
4 */ 4 */
5 #include <string.h> 5 #include <string.h>
6 #define _USE_MATH_DEFINES <<
7 #include <cmath> 6 #include <cmath>
8 7
>> 8 #ifdef WIN32
>> 9 //#include "erf.h"
>> 10 #include "amp_math.h"
>> 11 #define M_PI 3.141592653589793238463
>> 12 #endif
9 13
10 #include "MCGIDI_fromTOM.h" 14 #include "MCGIDI_fromTOM.h"
11 #include "MCGIDI_misc.h" 15 #include "MCGIDI_misc.h"
12 #include "MCGIDI_private.h" 16 #include "MCGIDI_private.h"
13 #include <nf_specialFunctions.h> 17 #include <nf_specialFunctions.h>
14 18
15 #if defined __cplusplus 19 #if defined __cplusplus
16 #include "G4Exp.hh" 20 #include "G4Exp.hh"
17 #include "G4Log.hh" 21 #include "G4Log.hh"
18 #include "G4Pow.hh" 22 #include "G4Pow.hh"
19 namespace GIDI { 23 namespace GIDI {
20 using namespace GIDI; 24 using namespace GIDI;
21 #endif 25 #endif
22 26
23 static int MCGIDI_energy_parseWeightFromTOM( s 27 static int MCGIDI_energy_parseWeightFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energyWeightedFunctional *weightedFunctional );
24 static int MCGIDI_energy_parseWeightedFunction 28 static int MCGIDI_energy_parseWeightedFunctionalsFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
25 static int MCGIDI_energy_parseGeneralEvaporati 29 static int MCGIDI_energy_parseGeneralEvaporationFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
26 static int MCGIDI_energy_parseEvaporationFromT 30 static int MCGIDI_energy_parseEvaporationFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
27 static int MCGIDI_energy_parseWattFromTOM( sta 31 static int MCGIDI_energy_parseWattFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
28 static int MCGIDI_energy_parseSimpleMaxwellian 32 static int MCGIDI_energy_parseSimpleMaxwellianFissionFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
29 static int MCGIDI_energy_parseMadlandNixFromTO 33 static int MCGIDI_energy_parseMadlandNixFromTOM( statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy );
30 static nfu_status MCGIDI_energy_parseMadlandNi 34 static nfu_status MCGIDI_energy_parseMadlandNixFromTOM_callback( double x, double *y, void *argList );
31 static double MCGIDI_energy_parseMadlandNixFro 35 static double MCGIDI_energy_parseMadlandNixFromTOM_callback_g( double Ep, double EFL, double T_M, nfu_status *status );
32 static int MCGIDI_energy_parseNBodyPhaseSpaceF 36 static int MCGIDI_energy_parseNBodyPhaseSpaceFromTOM( statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy,
33 MCGIDI_distribution *distribution ); 37 MCGIDI_distribution *distribution );
34 38
35 static int MCGIDI_energy_sampleSimpleMaxwellia 39 static int MCGIDI_energy_sampleSimpleMaxwellianFission( statusMessageReporting *smr, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo );
36 static int MCGIDI_energy_sampleEvaporation( st 40 static int MCGIDI_energy_sampleEvaporation( statusMessageReporting *smr, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo );
37 static int MCGIDI_energy_sampleWatt( statusMes 41 static int MCGIDI_energy_sampleWatt( statusMessageReporting *smr, double e_in_U, double Watt_a, double Watt_b, MCGIDI_decaySamplingInfo *decaySamplingInfo );
38 static int MCGIDI_energy_sampleWeightedFunctio 42 static int MCGIDI_energy_sampleWeightedFunctional( statusMessageReporting *smr, MCGIDI_energy *energy,
39 MCGIDI_quantitiesLookupModes &modes, MCGID 43 MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo );
40 static nfu_status MCGIDI_energy_NBodyPhaseSpac 44 static nfu_status MCGIDI_energy_NBodyPhaseSpacePDF_callback( double x, double *y, void *argList );
41 /* 45 /*
42 ********************************************** 46 ************************************************************
43 */ 47 */
44 MCGIDI_energy *MCGIDI_energy_new( statusMessag 48 MCGIDI_energy *MCGIDI_energy_new( statusMessageReporting *smr ) {
45 49
46 MCGIDI_energy *energy; 50 MCGIDI_energy *energy;
47 51
48 if( ( energy = (MCGIDI_energy *) smr_mallo 52 if( ( energy = (MCGIDI_energy *) smr_malloc2( smr, sizeof( MCGIDI_energy ), 0, "energy" ) ) == NULL ) return( NULL );
49 if( MCGIDI_energy_initialize( smr, energy 53 if( MCGIDI_energy_initialize( smr, energy ) ) energy = MCGIDI_energy_free( smr, energy );
50 return( energy ); 54 return( energy );
51 } 55 }
52 /* 56 /*
53 ********************************************** 57 ************************************************************
54 */ 58 */
55 int MCGIDI_energy_initialize( statusMessageRep 59 int MCGIDI_energy_initialize( statusMessageReporting * /*smr*/, MCGIDI_energy *energy ) {
56 60
57 memset( energy, 0, sizeof( MCGIDI_energy ) 61 memset( energy, 0, sizeof( MCGIDI_energy ) );
58 return( 0 ); 62 return( 0 );
59 } 63 }
60 /* 64 /*
61 ********************************************** 65 ************************************************************
62 */ 66 */
63 MCGIDI_energy *MCGIDI_energy_free( statusMessa 67 MCGIDI_energy *MCGIDI_energy_free( statusMessageReporting *smr, MCGIDI_energy *energy ) {
64 68
65 MCGIDI_energy_release( smr, energy ); 69 MCGIDI_energy_release( smr, energy );
66 smr_freeMemory( (void **) &energy ); 70 smr_freeMemory( (void **) &energy );
67 return( NULL ); 71 return( NULL );
68 } 72 }
69 /* 73 /*
70 ********************************************** 74 ************************************************************
71 */ 75 */
72 int MCGIDI_energy_release( statusMessageReport 76 int MCGIDI_energy_release( statusMessageReporting *smr, MCGIDI_energy *energy ) {
73 77
74 int i; 78 int i;
75 79
76 MCGIDI_sampling_pdfsOfXGivenW_release( smr 80 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(energy->dists) );
77 if( energy->theta ) energy->theta = ptwXY_ 81 if( energy->theta ) energy->theta = ptwXY_free( energy->theta );
78 if( energy->Watt_a ) energy->Watt_a = ptwX 82 if( energy->Watt_a ) energy->Watt_a = ptwXY_free( energy->Watt_a );
79 if( energy->Watt_b ) energy->Watt_b = ptwX 83 if( energy->Watt_b ) energy->Watt_b = ptwXY_free( energy->Watt_b );
80 if( ( energy->type == MCGIDI_energyType_ge 84 if( ( energy->type == MCGIDI_energyType_generalEvaporation ) || ( energy->type == MCGIDI_energyType_NBodyPhaseSpace ) ) {
81 MCGIDI_sampling_pdfsOfX_release( smr, 85 MCGIDI_sampling_pdfsOfX_release( smr, &(energy->g) ); }
82 else if( energy->type == MCGIDI_energyType 86 else if( energy->type == MCGIDI_energyType_weightedFunctional ) {
83 for( i = 0; i < energy->weightedFuncti 87 for( i = 0; i < energy->weightedFunctionals.numberOfWeights; i++ ) {
84 ptwXY_free( energy->weightedFuncti 88 ptwXY_free( energy->weightedFunctionals.weightedFunctional[i].weight );
85 MCGIDI_energy_free( smr, energy->w 89 MCGIDI_energy_free( smr, energy->weightedFunctionals.weightedFunctional[i].energy );
86 } 90 }
87 } 91 }
88 92
89 MCGIDI_energy_initialize( smr, energy ); 93 MCGIDI_energy_initialize( smr, energy );
90 return( 0 ); 94 return( 0 );
91 } 95 }
92 /* 96 /*
93 ********************************************** 97 ************************************************************
94 */ 98 */
95 int MCGIDI_energy_parseFromTOM( statusMessageR 99 int MCGIDI_energy_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution, ptwXYPoints *norms,
96 enum MCGIDI_energyType energyType, dou 100 enum MCGIDI_energyType energyType, double gammaEnergy_MeV ) {
97 101
98 MCGIDI_energy *energy = NULL; 102 MCGIDI_energy *energy = NULL;
99 xDataTOM_element *energyElement, *linearEl 103 xDataTOM_element *energyElement, *linearElement, *functional, *frameElement;
100 char const *nativeData; 104 char const *nativeData;
101 double projectileMass_MeV, targetMass_MeV; 105 double projectileMass_MeV, targetMass_MeV;
102 106
103 if( ( energy = MCGIDI_energy_new( smr ) ) 107 if( ( energy = MCGIDI_energy_new( smr ) ) == NULL ) goto err;
104 108
105 projectileMass_MeV = MCGIDI_product_getPro 109 projectileMass_MeV = MCGIDI_product_getProjectileMass_MeV( smr, distribution->product );
106 targetMass_MeV = MCGIDI_product_getTargetM 110 targetMass_MeV = MCGIDI_product_getTargetMass_MeV( smr, distribution->product );
107 energy->e_inCOMFactor = targetMass_MeV / ( 111 energy->e_inCOMFactor = targetMass_MeV / ( projectileMass_MeV + targetMass_MeV );
108 112
109 if( ( energyType == MCGIDI_energyType_prim 113 if( ( energyType == MCGIDI_energyType_primaryGamma ) || ( energyType == MCGIDI_energyType_discreteGamma ) ) {
110 energy->type = energyType; 114 energy->type = energyType;
111 energy->gammaEnergy_MeV = gammaEnergy_ 115 energy->gammaEnergy_MeV = gammaEnergy_MeV;
112 energy->frame = xDataTOM_frame_lab; 116 energy->frame = xDataTOM_frame_lab; /* BRB. This should not be hardwired?????? Probably needs to be changed in GND also. */
113 if( energyType == MCGIDI_energyType_pr 117 if( energyType == MCGIDI_energyType_primaryGamma ) energy->primaryGammaMassFactor = energy->e_inCOMFactor; }
114 else { 118 else {
115 if( ( energyElement = xDataTOME_getOne 119 if( ( energyElement = xDataTOME_getOneElementByName( smr, element, "energy", 1 ) ) == NULL ) goto err;
116 if( ( nativeData = xDataTOM_getAttribu 120 if( ( nativeData = xDataTOM_getAttributesValueInElement( energyElement, "nativeData" ) ) == NULL ) goto err;
117 if( ( linearElement = xDataTOME_getOne 121 if( ( linearElement = xDataTOME_getOneElementByName( NULL, energyElement, "linear", 0 ) ) == NULL )
118 linearElement = xDataTOME_getOneEl 122 linearElement = xDataTOME_getOneElementByName( NULL, energyElement, "pointwise", 0 );
119 if( linearElement == NULL ) { 123 if( linearElement == NULL ) {
120 if( ( functional = xDataTOME_getOn 124 if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "generalEvaporation", 0 ) ) != NULL ) {
121 if( MCGIDI_energy_parseGeneral 125 if( MCGIDI_energy_parseGeneralEvaporationFromTOM( smr, functional, energy ) ) goto err; }
122 else if( ( functional = xDataTOME_ 126 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "simpleMaxwellianFission", 0 ) ) != NULL ) {
123 if( MCGIDI_energy_parseSimpleM 127 if( MCGIDI_energy_parseSimpleMaxwellianFissionFromTOM( smr, functional, energy ) ) goto err; }
124 else if( ( functional = xDataTOME_ 128 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "evaporation", 0 ) ) != NULL ) {
125 if( MCGIDI_energy_parseEvapora 129 if( MCGIDI_energy_parseEvaporationFromTOM( smr, functional, energy ) ) goto err; }
126 else if( ( functional = xDataTOME_ 130 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "Watt", 0 ) ) != NULL ) {
127 if( MCGIDI_energy_parseWattFro 131 if( MCGIDI_energy_parseWattFromTOM( smr, functional, energy ) ) goto err; }
128 else if( ( functional = xDataTOME_ 132 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "MadlandNix", 0 ) ) != NULL ) {
129 if( MCGIDI_energy_parseMadland 133 if( MCGIDI_energy_parseMadlandNixFromTOM( smr, functional, energy ) ) goto err; }
130 else if( ( functional = xDataTOME_ 134 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "NBodyPhaseSpace", 0 ) ) != NULL ) {
131 if( MCGIDI_energy_parseNBodyPh 135 if( MCGIDI_energy_parseNBodyPhaseSpaceFromTOM( smr, functional, energy, distribution ) ) goto err; }
132 else if( ( functional = xDataTOME_ 136 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "weightedFunctionals", 0 ) ) != NULL ) {
133 if( MCGIDI_energy_parseWeighte 137 if( MCGIDI_energy_parseWeightedFunctionalsFromTOM( smr, functional, energy ) ) goto err; }
134 else { 138 else {
135 smr_setReportError2( smr, smr_ 139 smr_setReportError2( smr, smr_unknownID, 1, "unsupported energy type: nativeData = '%s'", nativeData );
136 goto err; 140 goto err;
137 } 141 }
138 frameElement = functional; } 142 frameElement = functional; }
139 else { 143 else {
140 char const *toUnits[3] = { "MeV", 144 char const *toUnits[3] = { "MeV", "MeV", "1/MeV" };
141 145
142 frameElement = linearElement; 146 frameElement = linearElement;
143 if( MCGIDI_fromTOM_pdfsOfXGivenW( 147 if( MCGIDI_fromTOM_pdfsOfXGivenW( smr, linearElement, &(energy->dists), norms, toUnits ) ) goto err;
144 energy->type = MCGIDI_energyType_l 148 energy->type = MCGIDI_energyType_linear;
145 } 149 }
146 if( ( energy->frame = MCGIDI_misc_getP 150 if( ( energy->frame = MCGIDI_misc_getProductFrame( smr, frameElement ) ) == xDataTOM_frame_invalid ) goto err;
147 } 151 }
148 distribution->energy = energy; 152 distribution->energy = energy;
149 153
150 return( 0 ); 154 return( 0 );
151 155
152 err: 156 err:
153 if( energy != NULL ) MCGIDI_energy_free( s 157 if( energy != NULL ) MCGIDI_energy_free( smr, energy );
154 return( 1 ); 158 return( 1 );
155 } 159 }
156 /* 160 /*
157 ********************************************** 161 ************************************************************
158 */ 162 */
159 static int MCGIDI_energy_parseWeightedFunction 163 static int MCGIDI_energy_parseWeightedFunctionalsFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
160 164
161 int i; 165 int i;
162 xDataTOM_element *child; 166 xDataTOM_element *child;
163 167
164 for( i = 0, child = xDataTOME_getFirstElem 168 for( i = 0, child = xDataTOME_getFirstElement( element ); child != NULL; i++, child = xDataTOME_getNextElement( child ) ) {
165 if( strcmp( child->name, "weighted" ) 169 if( strcmp( child->name, "weighted" ) ) goto err;
166 if( MCGIDI_energy_parseWeightFromTOM( 170 if( MCGIDI_energy_parseWeightFromTOM( smr, child, &(energy->weightedFunctionals.weightedFunctional[i]) ) ) goto err;
167 energy->weightedFunctionals.numberOfWe 171 energy->weightedFunctionals.numberOfWeights++;
168 } 172 }
169 energy->type = MCGIDI_energyType_weightedF 173 energy->type = MCGIDI_energyType_weightedFunctional;
170 return( 0 ); 174 return( 0 );
171 175
172 err: 176 err:
173 return( 1 ); 177 return( 1 );
174 } 178 }
175 /* 179 /*
176 ********************************************** 180 ************************************************************
177 */ 181 */
178 static int MCGIDI_energy_parseWeightFromTOM( s 182 static int MCGIDI_energy_parseWeightFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energyWeightedFunctional *weightedFunctional ) {
179 183
180 xDataTOM_element *child; 184 xDataTOM_element *child;
181 MCGIDI_energy *energy = NULL; 185 MCGIDI_energy *energy = NULL;
182 ptwXYPoints *weight = NULL; 186 ptwXYPoints *weight = NULL;
183 char const *toUnits[2] = { "MeV", "" }; 187 char const *toUnits[2] = { "MeV", "" };
184 188
185 if( ( energy = MCGIDI_energy_new( smr ) ) 189 if( ( energy = MCGIDI_energy_new( smr ) ) == NULL ) goto err;
186 for( child = xDataTOME_getFirstElement( el 190 for( child = xDataTOME_getFirstElement( element ); child != NULL; child = xDataTOME_getNextElement( child ) ) {
187 if( strcmp( child->name, "weight" ) == 191 if( strcmp( child->name, "weight" ) == 0 ) {
188 if( ( weight = MCGIDI_misc_dataFro 192 if( ( weight = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, child, toUnits ) ) == NULL ) goto err; }
189 else if( strcmp( child->name, "evapora 193 else if( strcmp( child->name, "evaporation" ) == 0 ) {
190 if( MCGIDI_energy_parseEvaporation 194 if( MCGIDI_energy_parseEvaporationFromTOM( smr, child, energy ) ) goto err; }
191 else { 195 else {
192 smr_setReportError2( smr, smr_unkn 196 smr_setReportError2( smr, smr_unknownID, 1, "unsupported energy type = '%s' in weighted functional", child->name );
193 goto err; 197 goto err;
194 } 198 }
195 } 199 }
196 weightedFunctional->weight = weight; 200 weightedFunctional->weight = weight;
197 weightedFunctional->energy = energy; 201 weightedFunctional->energy = energy;
198 return( 0 ); 202 return( 0 );
199 203
200 err: 204 err:
201 if( weight != NULL ) ptwXY_free( weight ); 205 if( weight != NULL ) ptwXY_free( weight );
202 if( energy != NULL ) MCGIDI_energy_free( s 206 if( energy != NULL ) MCGIDI_energy_free( smr, energy );
203 return( 1 ); 207 return( 1 );
204 } 208 }
205 /* 209 /*
206 ********************************************** 210 ************************************************************
207 */ 211 */
208 static int MCGIDI_energy_parseGeneralEvaporati 212 static int MCGIDI_energy_parseGeneralEvaporationFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
209 213
210 double norm; 214 double norm;
211 xDataTOM_element *thetaTOM, *gTOM; 215 xDataTOM_element *thetaTOM, *gTOM;
212 ptwXYPoints *theta = NULL, *g = NULL; 216 ptwXYPoints *theta = NULL, *g = NULL;
213 char const *toUnits[2] = { "MeV", "MeV" }; 217 char const *toUnits[2] = { "MeV", "MeV" };
214 218
215 if( ( thetaTOM = xDataTOME_getOneElementBy 219 if( ( thetaTOM = xDataTOME_getOneElementByName( smr, element, "theta", 1 ) ) == NULL ) goto err;
216 if( ( theta = MCGIDI_misc_dataFromElement2 220 if( ( theta = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, thetaTOM, toUnits ) ) == NULL ) goto err;
217 221
218 if( ( gTOM = xDataTOME_getOneElementByName 222 if( ( gTOM = xDataTOME_getOneElementByName( smr, element, "g", 1 ) ) == NULL ) goto err;
219 toUnits[0] = ""; 223 toUnits[0] = "";
220 toUnits[1] = ""; 224 toUnits[1] = "";
221 if( ( g = MCGIDI_misc_dataFromElement2ptwX 225 if( ( g = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, gTOM, toUnits ) ) == NULL ) goto err;
222 if( MCGIDI_fromTOM_pdfOfX( smr, g, &(energ 226 if( MCGIDI_fromTOM_pdfOfX( smr, g, &(energy->g), &norm ) ) goto err;
223 energy->gInterpolation = ptwXY_getInterpol 227 energy->gInterpolation = ptwXY_getInterpolation( g );
224 g = ptwXY_free( g ); 228 g = ptwXY_free( g );
225 if( std::fabs( 1. - norm ) > 0.001 ) print 229 if( std::fabs( 1. - norm ) > 0.001 ) printf( "bad norm = %e\n", norm );
226 230
227 energy->type = MCGIDI_energyType_generalEv 231 energy->type = MCGIDI_energyType_generalEvaporation;
228 energy->theta = theta; 232 energy->theta = theta;
229 return( 0 ); 233 return( 0 );
230 234
231 err: 235 err:
232 if( theta != NULL ) ptwXY_free( theta ); 236 if( theta != NULL ) ptwXY_free( theta );
233 if( g != NULL ) ptwXY_free( g ); 237 if( g != NULL ) ptwXY_free( g );
234 return( 1 ); 238 return( 1 );
235 } 239 }
236 /* 240 /*
237 ********************************************** 241 ************************************************************
238 */ 242 */
239 static int MCGIDI_energy_parseSimpleMaxwellian 243 static int MCGIDI_energy_parseSimpleMaxwellianFissionFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
240 244
241 char const *U, *toUnits[2] = { "MeV", "MeV 245 char const *U, *toUnits[2] = { "MeV", "MeV" };
242 xDataTOM_element *thetaTOM; 246 xDataTOM_element *thetaTOM;
243 247
244 if( ( U = xDataTOM_getAttributesValueInEle 248 if( ( U = xDataTOM_getAttributesValueInElement( element, "U" ) ) == NULL ) {
245 smr_setReportError2( smr, smr_unknownI 249 smr_setReportError2( smr, smr_unknownID, 1, "functional form '%s' missing 'U' attribute", element->name );
246 goto err; 250 goto err;
247 } 251 }
248 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 252 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, U, "MeV", &(energy->U) ) != 0 ) goto err;
249 if( ( thetaTOM = xDataTOME_getOneElementBy 253 if( ( thetaTOM = xDataTOME_getOneElementByName( smr, element, "theta", 1 ) ) == NULL ) goto err;
250 if( ( energy->theta = MCGIDI_misc_dataFrom 254 if( ( energy->theta = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, thetaTOM, toUnits ) ) == NULL ) goto err;
251 energy->type = MCGIDI_energyType_simpleMax 255 energy->type = MCGIDI_energyType_simpleMaxwellianFission;
252 return( 0 ); 256 return( 0 );
253 257
254 err: 258 err:
255 return( 1 ); 259 return( 1 );
256 } 260 }
257 /* 261 /*
258 ********************************************** 262 ************************************************************
259 */ 263 */
260 static int MCGIDI_energy_parseEvaporationFromT 264 static int MCGIDI_energy_parseEvaporationFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
261 265
262 char const *U, *toUnits[2] = { "MeV", "MeV 266 char const *U, *toUnits[2] = { "MeV", "MeV" };
263 xDataTOM_element *thetaTOM; 267 xDataTOM_element *thetaTOM;
264 268
265 if( ( U = xDataTOM_getAttributesValueInEle 269 if( ( U = xDataTOM_getAttributesValueInElement( element, "U" ) ) == NULL ) {
266 smr_setReportError2( smr, smr_unknownI 270 smr_setReportError2( smr, smr_unknownID, 1, "functional form '%s' missing 'U' attribute", element->name );
267 goto err; 271 goto err;
268 } 272 }
269 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 273 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, U, "MeV", &(energy->U) ) != 0 ) goto err;
270 if( ( thetaTOM = xDataTOME_getOneElementBy 274 if( ( thetaTOM = xDataTOME_getOneElementByName( smr, element, "theta", 1 ) ) == NULL ) goto err;
271 if( ( energy->theta = MCGIDI_misc_dataFrom 275 if( ( energy->theta = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, thetaTOM, toUnits ) ) == NULL ) goto err;
272 energy->type = MCGIDI_energyType_evaporati 276 energy->type = MCGIDI_energyType_evaporation;
273 return( 0 ); 277 return( 0 );
274 278
275 err: 279 err:
276 return( 1 ); 280 return( 1 );
277 } 281 }
278 /* 282 /*
279 ********************************************** 283 ************************************************************
280 */ 284 */
281 static int MCGIDI_energy_parseWattFromTOM( sta 285 static int MCGIDI_energy_parseWattFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
282 286
283 char const *U, *toUnits[2] = { "MeV", "MeV 287 char const *U, *toUnits[2] = { "MeV", "MeV" };
284 xDataTOM_element *aOrBTOM; 288 xDataTOM_element *aOrBTOM;
285 289
286 if( ( U = xDataTOM_getAttributesValueInEle 290 if( ( U = xDataTOM_getAttributesValueInElement( element, "U" ) ) == NULL ) {
287 smr_setReportError2( smr, smr_unknownI 291 smr_setReportError2( smr, smr_unknownID, 1, "functional form '%s' missing 'U' attribute", element->name );
288 goto err; 292 goto err;
289 } 293 }
290 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 294 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, U, "MeV", &(energy->U) ) != 0 ) goto err;
291 295
292 if( ( aOrBTOM = xDataTOME_getOneElementByN 296 if( ( aOrBTOM = xDataTOME_getOneElementByName( smr, element, "a", 1 ) ) == NULL ) goto err;
293 if( ( energy->Watt_a = MCGIDI_misc_dataFro 297 if( ( energy->Watt_a = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, aOrBTOM, toUnits ) ) == NULL ) goto err;
294 298
295 toUnits[1] = "1/MeV"; 299 toUnits[1] = "1/MeV";
296 if( ( aOrBTOM = xDataTOME_getOneElementByN 300 if( ( aOrBTOM = xDataTOME_getOneElementByName( smr, element, "b", 1 ) ) == NULL ) goto err;
297 if( ( energy->Watt_b = MCGIDI_misc_dataFro 301 if( ( energy->Watt_b = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, aOrBTOM, toUnits ) ) == NULL ) goto err;
298 302
299 energy->type = MCGIDI_energyType_Watt; 303 energy->type = MCGIDI_energyType_Watt;
300 return( 0 ); 304 return( 0 );
301 305
302 err: 306 err:
303 return( 1 ); 307 return( 1 );
304 } 308 }
305 /* 309 /*
306 ********************************************** 310 ************************************************************
307 */ 311 */
308 static int MCGIDI_energy_parseMadlandNixFromTO 312 static int MCGIDI_energy_parseMadlandNixFromTOM( statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy ) {
309 313
310 int iE, length, nXs, i1, n; 314 int iE, length, nXs, i1, n;
311 double E=0., T_M=0., EFL=0., EFH=0., argLi << 315 double E, T_M, EFL, EFH, argList[3], xs[] = { 1e-5, 1e-3, 1e-1, 1e1, 1e3, 1e5, 3e7 }, norm;
312 xs[] = { 1e-5, 1e-3, 1e-1, 1e1, 1e3 <<
313 ptwXYPoints *ptwXY_TM = NULL, *pdfXY = NUL 316 ptwXYPoints *ptwXY_TM = NULL, *pdfXY = NULL;
314 ptwXYPoint *point; 317 ptwXYPoint *point;
315 ptwXPoints *cdfX = NULL; 318 ptwXPoints *cdfX = NULL;
316 nfu_status status = nfu_Okay; 319 nfu_status status = nfu_Okay;
317 xDataTOM_element *TM_TOM; 320 xDataTOM_element *TM_TOM;
318 xDataTOM_XYs *XYs; 321 xDataTOM_XYs *XYs;
319 MCGIDI_pdfsOfXGivenW *dists = &(energy->di 322 MCGIDI_pdfsOfXGivenW *dists = &(energy->dists);
320 MCGIDI_pdfOfX *dist; 323 MCGIDI_pdfOfX *dist;
321 char const *EF, *TMUnits[2] = { "MeV", "Me 324 char const *EF, *TMUnits[2] = { "MeV", "MeV" };
322 325
323 nXs = sizeof( xs ) / sizeof( xs[0] ); 326 nXs = sizeof( xs ) / sizeof( xs[0] );
324 327
325 if( ( EF = xDataTOM_getAttributesValueInEl 328 if( ( EF = xDataTOM_getAttributesValueInElement( functional, "EFL" ) ) == NULL ) {
326 smr_setReportError2( smr, smr_unknownI 329 smr_setReportError2( smr, smr_unknownID, 1, "MadlandNix '%s' missing 'EFL' attribute", functional->name );
327 goto err; 330 goto err;
328 } 331 }
329 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 332 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, EF, TMUnits[0], &EFL ) != 0 ) goto err;
330 argList[0] = EFL; 333 argList[0] = EFL;
331 334
332 if( ( EF = xDataTOM_getAttributesValueInEl 335 if( ( EF = xDataTOM_getAttributesValueInElement( functional, "EFH" ) ) == NULL ) {
333 smr_setReportError2( smr, smr_unknownI 336 smr_setReportError2( smr, smr_unknownID, 1, "MadlandNix '%s' missing 'EFH' attribute", functional->name );
334 goto err; 337 goto err;
335 } 338 }
336 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 339 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, EF, TMUnits[0], &EFH ) != 0 ) goto err;
337 argList[1] = EFH; 340 argList[1] = EFH;
338 341
339 if( ( TM_TOM = xDataTOME_getOneElementByNa 342 if( ( TM_TOM = xDataTOME_getOneElementByName( smr, functional, "T_M", 1 ) ) == NULL ) goto err;
340 if( ( XYs = (xDataTOM_XYs *) xDataTOME_get 343 if( ( XYs = (xDataTOM_XYs *) xDataTOME_getXDataIfID( smr, TM_TOM, "XYs" ) ) == NULL ) goto err;
341 if( ( ptwXY_TM = MCGIDI_misc_dataFromXYs2p 344 if( ( ptwXY_TM = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, ptwXY_interpolationLinLin, TMUnits ) ) == NULL ) goto err;
342 345
343 length = (int) ptwXY_length( ptwXY_TM ); 346 length = (int) ptwXY_length( ptwXY_TM );
344 dists->interpolationWY = ptwXY_interpolati 347 dists->interpolationWY = ptwXY_interpolationLinLin;
345 dists->interpolationXY = ptwXY_interpolati 348 dists->interpolationXY = ptwXY_interpolationLinLin; /* Ignoring what the data says as it is probably wrong. */
346 if( ( dists->Ws = (double *) smr_malloc2( 349 if( ( dists->Ws = (double *) smr_malloc2( smr, length * sizeof( double ), 1, "dists->Ws" ) ) == NULL ) goto err;
347 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_ 350 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, length * sizeof( MCGIDI_pdfOfX ), 0, "dists->dist" ) ) == NULL ) goto err;
348 351
349 for( iE = 0; iE < length; iE++ ) { 352 for( iE = 0; iE < length; iE++ ) {
350 ptwXY_getXYPairAtIndex( ptwXY_TM, iE, 353 ptwXY_getXYPairAtIndex( ptwXY_TM, iE, &E, &T_M );
351 argList[2] = T_M; 354 argList[2] = T_M;
352 dist = &(dists->dist[iE]); 355 dist = &(dists->dist[iE]);
353 dists->Ws[iE] = E; 356 dists->Ws[iE] = E;
354 357
355 if( ( pdfXY = ptwXY_createFromFunction 358 if( ( pdfXY = ptwXY_createFromFunction( nXs, xs, (ptwXY_createFromFunction_callback) MCGIDI_energy_parseMadlandNixFromTOM_callback,
356 (void *) argList, 1e-3, 0, 12, &st 359 (void *) argList, 1e-3, 0, 12, &status ) ) == NULL ) goto err;
357 if( ( status = ptwXY_normalize( pdfXY 360 if( ( status = ptwXY_normalize( pdfXY ) ) != nfu_Okay ) {
358 smr_setReportError2( smr, smr_unkn 361 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_normalize err = %d: %s\n", status, nfu_statusMessage( status ) );
359 goto err; 362 goto err;
360 } 363 }
361 364
362 if( ptwXY_simpleCoalescePoints( pdfXY 365 if( ptwXY_simpleCoalescePoints( pdfXY ) != nfu_Okay ) goto err;
363 dist->numberOfXs = n = (int) ptwXY_len 366 dist->numberOfXs = n = (int) ptwXY_length( pdfXY );
364 367
365 if( ( dist->Xs = (double *) smr_malloc 368 if( ( dist->Xs = (double *) smr_malloc2( smr, 3 * n * sizeof( double ), 0, "dist->Xs" ) ) == NULL ) goto err;
366 dists->numberOfWs++; 369 dists->numberOfWs++;
367 dist->pdf = &(dist->Xs[n]); 370 dist->pdf = &(dist->Xs[n]);
368 dist->cdf = &(dist->pdf[n]); 371 dist->cdf = &(dist->pdf[n]);
369 372
370 for( i1 = 0; i1 < n; i1++ ) { 373 for( i1 = 0; i1 < n; i1++ ) {
371 point = ptwXY_getPointAtIndex_Unsa 374 point = ptwXY_getPointAtIndex_Unsafely( pdfXY, i1 );
372 dist->Xs[i1] = point->x; 375 dist->Xs[i1] = point->x;
373 dist->pdf[i1] = point->y; 376 dist->pdf[i1] = point->y;
374 } 377 }
375 378
376 if( ( cdfX = ptwXY_runningIntegral( pd 379 if( ( cdfX = ptwXY_runningIntegral( pdfXY, &status ) ) == NULL ) {
377 smr_setReportError2( smr, smr_unkn 380 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_runningIntegral err = %d: %s\n", status, nfu_statusMessage( status ) );
378 goto err; 381 goto err;
379 } 382 }
380 383
381 norm = ptwX_getPointAtIndex_Unsafely( 384 norm = ptwX_getPointAtIndex_Unsafely( cdfX, n - 1 );
382 for( i1 = 0; i1 < n; i1++ ) dist->cdf[ 385 for( i1 = 0; i1 < n; i1++ ) dist->cdf[i1] = ptwX_getPointAtIndex_Unsafely( cdfX, i1 ) / norm;
383 for( i1 = 0; i1 < n; i1++ ) dist->pdf[ 386 for( i1 = 0; i1 < n; i1++ ) dist->pdf[i1] /= norm;
384 pdfXY = ptwXY_free( pdfXY ); 387 pdfXY = ptwXY_free( pdfXY );
385 cdfX = ptwX_free( cdfX ); 388 cdfX = ptwX_free( cdfX );
386 } 389 }
387 390
388 energy->type = MCGIDI_energyType_MadlandNi 391 energy->type = MCGIDI_energyType_MadlandNix;
389 392
390 ptwXY_free( ptwXY_TM ); 393 ptwXY_free( ptwXY_TM );
391 return( 0 ); 394 return( 0 );
392 395
393 err: 396 err:
394 if( ptwXY_TM != NULL ) ptwXY_free( ptwXY_T 397 if( ptwXY_TM != NULL ) ptwXY_free( ptwXY_TM );
395 if( pdfXY != NULL ) ptwXY_free( pdfXY ); 398 if( pdfXY != NULL ) ptwXY_free( pdfXY );
396 if( cdfX != NULL ) cdfX = ptwX_free( cdfX 399 if( cdfX != NULL ) cdfX = ptwX_free( cdfX );
397 400
398 return( 1 ); 401 return( 1 );
399 } 402 }
400 /* 403 /*
401 ********************************************** 404 ************************************************************
402 */ 405 */
403 static nfu_status MCGIDI_energy_parseMadlandNi 406 static nfu_status MCGIDI_energy_parseMadlandNixFromTOM_callback( double Ep, double *y, void *argList ) {
404 407
405 double *parameters = (double *) argList, E 408 double *parameters = (double *) argList, EFL, EFH, T_M;
406 nfu_status status = nfu_Okay; 409 nfu_status status = nfu_Okay;
407 410
408 EFL = parameters[0]; 411 EFL = parameters[0];
409 EFH = parameters[1]; 412 EFH = parameters[1];
410 T_M = parameters[2]; 413 T_M = parameters[2];
411 *y = MCGIDI_energy_parseMadlandNixFromTOM_ 414 *y = MCGIDI_energy_parseMadlandNixFromTOM_callback_g( Ep, EFL, T_M, &status );
412 if( status == nfu_Okay ) *y += MCGIDI_ener 415 if( status == nfu_Okay ) *y += MCGIDI_energy_parseMadlandNixFromTOM_callback_g( Ep, EFH, T_M, &status );
413 *y *= 0.5; 416 *y *= 0.5;
414 return( status ); 417 return( status );
415 } 418 }
416 /* 419 /*
417 ********************************************** 420 ************************************************************
418 */ 421 */
419 static double MCGIDI_energy_parseMadlandNixFro 422 static double MCGIDI_energy_parseMadlandNixFromTOM_callback_g( double Ep, double E_F, double T_M, nfu_status *status ) {
420 423
421 double u1, u2, E1, E2 = 0., gamma1 = 0., g << 424 double u1, u2, E1, E2, gamma1, gamma2, signG = 1;
422 425
423 u1 = std::sqrt( Ep ) - std::sqrt( E_F ); 426 u1 = std::sqrt( Ep ) - std::sqrt( E_F );
424 u1 *= u1 / T_M; 427 u1 *= u1 / T_M;
425 u2 = std::sqrt( Ep ) + std::sqrt( E_F ); 428 u2 = std::sqrt( Ep ) + std::sqrt( E_F );
426 u2 *= u2 / T_M; 429 u2 *= u2 / T_M;
427 E1 = 0; /* u1^3/2 * E 430 E1 = 0; /* u1^3/2 * E1 is zero for u1 = 0. but E1 is infinity, whence, the next test. */
428 if( u1 != 0 ) E1 = nf_exponentialIntegral( 431 if( u1 != 0 ) E1 = nf_exponentialIntegral( 1, u1, status );
429 if( *status == nfu_Okay ) E2 = nf_exponent 432 if( *status == nfu_Okay ) E2 = nf_exponentialIntegral( 1, u2, status );
430 if( *status != nfu_Okay ) return( 0. ); 433 if( *status != nfu_Okay ) return( 0. );
431 if( u1 > 2. ) { 434 if( u1 > 2. ) {
432 signG = -1; 435 signG = -1;
433 gamma1 = nf_incompleteGammaFunctionCom 436 gamma1 = nf_incompleteGammaFunctionComplementary( 1.5, u1, status );
434 if( *status == nfu_Okay ) gamma2 = nf_ 437 if( *status == nfu_Okay ) gamma2 = nf_incompleteGammaFunctionComplementary( 1.5, u2, status ); }
435 else { 438 else {
436 gamma1 = nf_incompleteGammaFunction( 1 439 gamma1 = nf_incompleteGammaFunction( 1.5, u1, status );
437 if( *status == nfu_Okay ) gamma2 = nf_ 440 if( *status == nfu_Okay ) gamma2 = nf_incompleteGammaFunction( 1.5, u2, status );
438 } 441 }
439 if( *status != nfu_Okay ) return( 0. ); 442 if( *status != nfu_Okay ) return( 0. );
440 return( ( u2 * std::sqrt( u2 ) * E2 - u1 * 443 return( ( u2 * std::sqrt( u2 ) * E2 - u1 * std::sqrt( u1 ) * E1 + signG * ( gamma2 - gamma1 ) ) / ( 3 * std::sqrt( E_F * T_M ) ) );
441 } 444 }
442 /* 445 /*
443 ********************************************** 446 ************************************************************
444 */ 447 */
445 static int MCGIDI_energy_parseNBodyPhaseSpaceF 448 static int MCGIDI_energy_parseNBodyPhaseSpaceFromTOM( statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy,
446 MCGIDI_distribution *distribution ) { 449 MCGIDI_distribution *distribution ) {
447 450
448 int argList[1]; 451 int argList[1];
449 double xs[2] = { 0.0, 1.0 }, productMass_M 452 double xs[2] = { 0.0, 1.0 }, productMass_MeV, norm;
450 ptwXYPoints *pdf = NULL; 453 ptwXYPoints *pdf = NULL;
451 nfu_status status; 454 nfu_status status;
452 char const *mass; 455 char const *mass;
453 456
454 if( xDataTOME_convertAttributeToInteger( N 457 if( xDataTOME_convertAttributeToInteger( NULL, functional, "numberOfProducts", &(energy->NBodyPhaseSpace.numberOfProducts) ) != 0 ) goto err;
455 if( ( mass = xDataTOM_getAttributesValueIn 458 if( ( mass = xDataTOM_getAttributesValueInElement( functional, "mass" ) ) == NULL ) {
456 smr_setReportError2( smr, smr_unknownI 459 smr_setReportError2( smr, smr_unknownID, 1, "functional form '%s' missing 'mass' attribute", functional->name );
457 goto err; 460 goto err;
458 } 461 }
459 if( MCGIDI_misc_PQUStringToDouble( smr, ma 462 if( MCGIDI_misc_PQUStringToDouble( smr, mass, "amu", MCGIDI_AMU2MeV, &(energy->NBodyPhaseSpace.mass) ) ) goto err;
460 argList[0] = energy->NBodyPhaseSpace.numbe 463 argList[0] = energy->NBodyPhaseSpace.numberOfProducts;
461 if( ( pdf = ptwXY_createFromFunction( 2, x 464 if( ( pdf = ptwXY_createFromFunction( 2, xs, MCGIDI_energy_NBodyPhaseSpacePDF_callback, (void *) argList, 1e-3, 0, 16, &status ) ) == NULL ) {
462 smr_setReportError2( smr, smr_unknownI 465 smr_setReportError2( smr, smr_unknownID, 1, "creating NBodyPhaseSpace pdf failed with ptwXY_createFromFunction error = %d (%s)",
463 status, nfu_statusMessage( status 466 status, nfu_statusMessage( status ) );
464 goto err; 467 goto err;
465 } 468 }
466 if( MCGIDI_fromTOM_pdfOfX( smr, pdf, &(ene 469 if( MCGIDI_fromTOM_pdfOfX( smr, pdf, &(energy->g), &norm ) ) goto err;
467 productMass_MeV = MCGIDI_product_getMass_M 470 productMass_MeV = MCGIDI_product_getMass_MeV( smr, distribution->product );
468 if( !smr_isOk( smr ) ) goto err; 471 if( !smr_isOk( smr ) ) goto err;
469 energy->NBodyPhaseSpace.massFactor = ( 1. 472 energy->NBodyPhaseSpace.massFactor = ( 1. - productMass_MeV / ( MCGIDI_AMU2MeV * energy->NBodyPhaseSpace.mass ) ); /* ??????? Hardwired MCGIDI_AMU2MeV */
470 energy->NBodyPhaseSpace.Q_MeV = MCGIDI_out 473 energy->NBodyPhaseSpace.Q_MeV = MCGIDI_outputChannel_getQ_MeV( smr, distribution->product->outputChannel, 0. );
471 if( !smr_isOk( smr ) ) goto err; 474 if( !smr_isOk( smr ) ) goto err;
472 475
473 ptwXY_free( pdf ); 476 ptwXY_free( pdf );
474 energy->type = MCGIDI_energyType_NBodyPhas 477 energy->type = MCGIDI_energyType_NBodyPhaseSpace;
475 478
476 return( 0 ); 479 return( 0 );
477 480
478 err: 481 err:
479 if( pdf != NULL ) ptwXY_free( pdf ); 482 if( pdf != NULL ) ptwXY_free( pdf );
480 return( 1 ); 483 return( 1 );
481 } 484 }
482 /* 485 /*
483 ********************************************** 486 ************************************************************
484 */ 487 */
485 static nfu_status MCGIDI_energy_NBodyPhaseSpac 488 static nfu_status MCGIDI_energy_NBodyPhaseSpacePDF_callback( double x, double *y, void *argList ) {
486 489
487 int numberOfProducts = *((int *) argList); 490 int numberOfProducts = *((int *) argList);
488 double e = 0.5 * ( 3 * numberOfProducts - 491 double e = 0.5 * ( 3 * numberOfProducts - 8 );
489 492
490 *y = std::sqrt( x ) * G4Pow::GetInstance() 493 *y = std::sqrt( x ) * G4Pow::GetInstance()->powA( 1.0 - x, e );
491 return( nfu_Okay ); 494 return( nfu_Okay );
492 } 495 }
493 /* 496 /*
494 ********************************************** 497 ************************************************************
495 */ 498 */
496 int MCGIDI_energy_sampleEnergy( statusMessageR 499 int MCGIDI_energy_sampleEnergy( statusMessageReporting *smr, MCGIDI_energy *energy, MCGIDI_quantitiesLookupModes &modes,
497 MCGIDI_decaySamplingInfo *decaySamplin 500 MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
498 /* 501 /*
499 * This function must be called before angula 502 * This function must be called before angular sampling as it sets the frame but does not test it.
500 */ 503 */
501 double theta, randomEp, Watt_a, Watt_b, e_ 504 double theta, randomEp, Watt_a, Watt_b, e_in = modes.getProjectileEnergy( );
502 MCGIDI_pdfsOfXGivenW_sampled sampled; 505 MCGIDI_pdfsOfXGivenW_sampled sampled;
503 506
504 decaySamplingInfo->frame = energy->frame; 507 decaySamplingInfo->frame = energy->frame;
505 switch( energy->type ) { 508 switch( energy->type ) {
506 case MCGIDI_energyType_primaryGamma : 509 case MCGIDI_energyType_primaryGamma :
507 decaySamplingInfo->Ep = energy->gammaE 510 decaySamplingInfo->Ep = energy->gammaEnergy_MeV + e_in * energy->primaryGammaMassFactor;
508 break; 511 break;
509 case MCGIDI_energyType_discreteGamma : 512 case MCGIDI_energyType_discreteGamma :
510 decaySamplingInfo->Ep = energy->gammaE 513 decaySamplingInfo->Ep = energy->gammaEnergy_MeV;
511 break; 514 break;
512 case MCGIDI_energyType_linear : 515 case MCGIDI_energyType_linear :
513 randomEp = decaySamplingInfo->rng( dec 516 randomEp = decaySamplingInfo->rng( decaySamplingInfo->rngState );
514 sampled.smr = smr; 517 sampled.smr = smr;
515 sampled.w = e_in; 518 sampled.w = e_in;
516 MCGIDI_sampling_sampleX_from_pdfsOfXGi 519 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(energy->dists), &sampled, randomEp );
517 decaySamplingInfo->Ep = sampled.x; 520 decaySamplingInfo->Ep = sampled.x;
518 break; 521 break;
519 case MCGIDI_energyType_generalEvaporation 522 case MCGIDI_energyType_generalEvaporation :
520 sampled.interpolationXY = energy->gInt 523 sampled.interpolationXY = energy->gInterpolation;
521 MCGIDI_sampling_sampleX_from_pdfOfX( & 524 MCGIDI_sampling_sampleX_from_pdfOfX( &(energy->g), &sampled, decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
522 theta = MCGIDI_sampling_ptwXY_getValue 525 theta = MCGIDI_sampling_ptwXY_getValueAtX( energy->theta, e_in );
523 decaySamplingInfo->Ep = theta * sample 526 decaySamplingInfo->Ep = theta * sampled.x;
524 break; 527 break;
525 case MCGIDI_energyType_simpleMaxwellianFis 528 case MCGIDI_energyType_simpleMaxwellianFission :
526 theta = MCGIDI_sampling_ptwXY_getValue 529 theta = MCGIDI_sampling_ptwXY_getValueAtX( energy->theta, e_in );
527 MCGIDI_energy_sampleSimpleMaxwellianFi 530 MCGIDI_energy_sampleSimpleMaxwellianFission( smr, ( e_in - energy->U ) / theta, decaySamplingInfo );
528 decaySamplingInfo->Ep *= theta; 531 decaySamplingInfo->Ep *= theta;
529 break; 532 break;
530 case MCGIDI_energyType_evaporation : 533 case MCGIDI_energyType_evaporation :
531 theta = MCGIDI_sampling_ptwXY_getValue 534 theta = MCGIDI_sampling_ptwXY_getValueAtX( energy->theta, e_in );
532 MCGIDI_energy_sampleEvaporation( smr, 535 MCGIDI_energy_sampleEvaporation( smr, ( e_in - energy->U ) / theta, decaySamplingInfo );
533 decaySamplingInfo->Ep *= theta; 536 decaySamplingInfo->Ep *= theta;
534 break; 537 break;
535 case MCGIDI_energyType_Watt : 538 case MCGIDI_energyType_Watt :
536 Watt_a = MCGIDI_sampling_ptwXY_getValu 539 Watt_a = MCGIDI_sampling_ptwXY_getValueAtX( energy->Watt_a, e_in );
537 Watt_b = MCGIDI_sampling_ptwXY_getValu 540 Watt_b = MCGIDI_sampling_ptwXY_getValueAtX( energy->Watt_b, e_in );
538 MCGIDI_energy_sampleWatt( smr, e_in - 541 MCGIDI_energy_sampleWatt( smr, e_in - energy->U, Watt_a, Watt_b, decaySamplingInfo );
539 break; 542 break;
540 case MCGIDI_energyType_MadlandNix : 543 case MCGIDI_energyType_MadlandNix :
541 MCGIDI_sampling_sampleX_from_pdfsOfXGi 544 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(energy->dists), &sampled, decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
542 decaySamplingInfo->Ep = sampled.x; 545 decaySamplingInfo->Ep = sampled.x;
543 break; 546 break;
544 case MCGIDI_energyType_NBodyPhaseSpace : 547 case MCGIDI_energyType_NBodyPhaseSpace :
545 MCGIDI_sampling_sampleX_from_pdfOfX( & 548 MCGIDI_sampling_sampleX_from_pdfOfX( &(energy->g), &sampled, decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
546 decaySamplingInfo->Ep = ( energy->e_in 549 decaySamplingInfo->Ep = ( energy->e_inCOMFactor * e_in + energy->NBodyPhaseSpace.Q_MeV ) * energy->NBodyPhaseSpace.massFactor * sampled.x;
547 break; 550 break;
548 case MCGIDI_energyType_weightedFunctional 551 case MCGIDI_energyType_weightedFunctional :
549 MCGIDI_energy_sampleWeightedFunctional 552 MCGIDI_energy_sampleWeightedFunctional( smr, energy, modes, decaySamplingInfo );
550 break; 553 break;
551 default : 554 default :
552 smr_setReportError2( smr, smr_unknownI 555 smr_setReportError2( smr, smr_unknownID, 1, "energy type = %d not supported", energy->type );
553 } 556 }
554 557
555 return( !smr_isOk( smr ) ); 558 return( !smr_isOk( smr ) );
556 } 559 }
557 /* 560 /*
558 ********************************************** 561 ************************************************************
559 */ 562 */
560 static int MCGIDI_energy_sampleSimpleMaxwellia 563 static int MCGIDI_energy_sampleSimpleMaxwellianFission( statusMessageReporting * /*smr*/, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
561 564
562 int i1; 565 int i1;
563 double a = e_in_U_theta, b, c, x, norm_a, 566 double a = e_in_U_theta, b, c, x, norm_a, xMin = 0., xMax = a, sqrt_x, sqrt_pi_2 = std::sqrt( M_PI ) / 2.;
564 567
565 sqrt_x = std::sqrt( a ); 568 sqrt_x = std::sqrt( a );
566 norm_a = sqrt_pi_2 * erf( sqrt_x ) - sqrt_ 569 norm_a = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x * G4Exp( -a );
567 b = norm_a * decaySamplingInfo->rng( decay 570 b = norm_a * decaySamplingInfo->rng( decaySamplingInfo->rngState );
568 for( i1 = 0; i1 < 16; i1++ ) { 571 for( i1 = 0; i1 < 16; i1++ ) {
569 x = 0.5 * ( xMin + xMax ); 572 x = 0.5 * ( xMin + xMax );
570 sqrt_x = std::sqrt( x ); 573 sqrt_x = std::sqrt( x );
571 c = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x 574 c = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x * G4Exp( -x );
572 if( b < c ) { 575 if( b < c ) {
573 xMax = x; } 576 xMax = x; }
574 else { 577 else {
575 xMin = x; 578 xMin = x;
576 } 579 }
577 } 580 }
578 /* To order e, the correct x is x + e 581 /* To order e, the correct x is x + e where e = 1 + ( 1 - b * exp( x ) ) / x. */
579 decaySamplingInfo->Ep = x; 582 decaySamplingInfo->Ep = x;
580 583
581 return( 0 ); 584 return( 0 );
582 } 585 }
583 /* 586 /*
584 ********************************************** 587 ************************************************************
585 */ 588 */
586 static int MCGIDI_energy_sampleEvaporation( st 589 static int MCGIDI_energy_sampleEvaporation( statusMessageReporting * /*smr*/, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
587 590
588 int i1; 591 int i1;
589 double a = e_in_U_theta, b, c, x, norm_a, 592 double a = e_in_U_theta, b, c, x, norm_a, xMin = 0., xMax = a;
590 593
591 norm_a = 1 - ( 1 + a ) * G4Exp( -a ); 594 norm_a = 1 - ( 1 + a ) * G4Exp( -a );
592 b = 1. - norm_a * decaySamplingInfo->rng( 595 b = 1. - norm_a * decaySamplingInfo->rng( decaySamplingInfo->rngState );
593 for( i1 = 0; i1 < 16; i1++ ) { 596 for( i1 = 0; i1 < 16; i1++ ) {
594 x = 0.5 * ( xMin + xMax ); 597 x = 0.5 * ( xMin + xMax );
595 c = ( 1 + x ) * G4Exp( -x ); 598 c = ( 1 + x ) * G4Exp( -x );
596 if( b > c ) { 599 if( b > c ) {
597 xMax = x; } 600 xMax = x; }
598 else { 601 else {
599 xMin = x; 602 xMin = x;
600 } 603 }
601 } 604 }
602 /* To order e, the correct x is x + e 605 /* To order e, the correct x is x + e where e = 1 + ( 1 - b * std::exp( x ) ) / x. */
603 decaySamplingInfo->Ep = x; 606 decaySamplingInfo->Ep = x;
604 607
605 return( 0 ); 608 return( 0 );
606 } 609 }
607 /* 610 /*
608 ********************************************** 611 ************************************************************
609 */ 612 */
610 static int MCGIDI_energy_sampleWatt( statusMes 613 static int MCGIDI_energy_sampleWatt( statusMessageReporting * /*smr*/, double e_in_U, double Watt_a, double Watt_b, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
611 /* 614 /*
612 * From MCAPM via Sample Watt Spectrum as in 615 * From MCAPM via Sample Watt Spectrum as in TART ( Kalos algorithm ).
613 */ 616 */
614 double WattMin = 0., WattMax = e_in_U, x, << 617 double WattMin = 0., WattMax = e_in_U, x, y, z, energyOut, rand1, rand2;
615 618
616 x = 1. + ( Watt_b / ( 8. * Watt_a ) ); 619 x = 1. + ( Watt_b / ( 8. * Watt_a ) );
617 y = ( x + std::sqrt( x * x - 1. ) ) / Watt 620 y = ( x + std::sqrt( x * x - 1. ) ) / Watt_a;
618 z = Watt_a * y - 1.; 621 z = Watt_a * y - 1.;
619 G4int icounter=0; 622 G4int icounter=0;
620 G4int icounter_max=1024; 623 G4int icounter_max=1024;
621 do 624 do
622 { 625 {
623 icounter++; 626 icounter++;
624 if ( icounter > icounter_max ) { 627 if ( icounter > icounter_max ) {
625 G4cout << "Loop-counter exceeded the 628 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
626 break; 629 break;
627 } 630 }
628 rand1 = -G4Log( decaySamplingInfo->rng 631 rand1 = -G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
629 rand2 = -G4Log( decaySamplingInfo->rng 632 rand2 = -G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
630 energyOut = y * rand1; 633 energyOut = y * rand1;
631 } 634 }
632 while( ( ( rand2 - z * ( rand1 + 1. ) ) * 635 while( ( ( rand2 - z * ( rand1 + 1. ) ) * ( rand2 - z * ( rand1 + 1. ) ) > Watt_b * y * rand1 ) || ( energyOut < WattMin ) || ( energyOut > WattMax ) ); // Loop checking, 11.06.2015, T. Koi
633 decaySamplingInfo->Ep = energyOut; 636 decaySamplingInfo->Ep = energyOut;
634 637
635 return( 0 ); 638 return( 0 );
636 } 639 }
637 /* 640 /*
638 ********************************************** 641 ************************************************************
639 */ 642 */
640 static int MCGIDI_energy_sampleWeightedFunctio 643 static int MCGIDI_energy_sampleWeightedFunctional( statusMessageReporting *smr, MCGIDI_energy *energy,
641 MCGIDI_quantitiesLookupModes &modes, M 644 MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
642 /* 645 /*
643 c This routine assumes that the weights sum 646 c This routine assumes that the weights sum to 1.
644 */ 647 */
645 int iW; 648 int iW;
646 double rW = decaySamplingInfo->rng( decayS 649 double rW = decaySamplingInfo->rng( decaySamplingInfo->rngState ), cumulativeW = 0., weight;
647 MCGIDI_energyWeightedFunctional *weightedF 650 MCGIDI_energyWeightedFunctional *weightedFunctional = NULL;
648 651
649 for( iW = 0; iW < energy->weightedFunction 652 for( iW = 0; iW < energy->weightedFunctionals.numberOfWeights; iW++ ) {
650 weightedFunctional = &(energy->weighte 653 weightedFunctional = &(energy->weightedFunctionals.weightedFunctional[iW]);
651 weight = MCGIDI_sampling_ptwXY_getValu 654 weight = MCGIDI_sampling_ptwXY_getValueAtX( weightedFunctional->weight, modes.getProjectileEnergy( ) );
652 cumulativeW += weight; 655 cumulativeW += weight;
653 if( cumulativeW >= rW ) break; 656 if( cumulativeW >= rW ) break;
654 } 657 }
655 return( MCGIDI_energy_sampleEnergy( smr, w 658 return( MCGIDI_energy_sampleEnergy( smr, weightedFunctional->energy, modes, decaySamplingInfo ) );
656 } 659 }
657 660
658 #if defined __cplusplus 661 #if defined __cplusplus
659 } 662 }
660 #endif 663 #endif
661 664
662 665