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