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