11 #define M_PI 3.141592653589793238463
14 #include "MCGIDI_fromTOM.h"
15 #include "MCGIDI_misc.h"
16 #include "MCGIDI_private.h"
17 #include <nf_specialFunctions.h>
19 #if defined __cplusplus
37 MCGIDI_distribution *distribution );
41 static int MCGIDI_energy_sampleWatt( statusMessageReporting *smr,
double e_in_U,
double Watt_a,
double Watt_b, MCGIDI_decaySamplingInfo *decaySamplingInfo );
43 MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo );
52 if( ( energy = (MCGIDI_energy *) smr_malloc2( smr,
sizeof( MCGIDI_energy ), 0,
"energy" ) ) == NULL )
return( NULL );
61 memset( energy, 0,
sizeof( MCGIDI_energy ) );
81 if( energy->theta ) energy->theta =
ptwXY_free( energy->theta );
82 if( energy->Watt_a ) energy->Watt_a =
ptwXY_free( energy->Watt_a );
83 if( energy->Watt_b ) energy->Watt_b =
ptwXY_free( energy->Watt_b );
84 if( ( energy->type == MCGIDI_energyType_generalEvaporation ) || ( energy->type == MCGIDI_energyType_NBodyPhaseSpace ) ) {
86 else if( energy->type == MCGIDI_energyType_weightedFunctional ) {
87 for( i = 0; i < energy->weightedFunctionals.numberOfWeights; i++ ) {
88 ptwXY_free( energy->weightedFunctionals.weightedFunctional[i].weight );
99 int MCGIDI_energy_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution, ptwXYPoints *norms,
100 enum MCGIDI_energyType energyType,
double gammaEnergy_MeV ) {
102 MCGIDI_energy *
energy = NULL;
103 xDataTOM_element *energyElement, *linearElement, *functional, *frameElement;
104 char const *nativeData;
105 double projectileMass_MeV, targetMass_MeV;
111 energy->e_inCOMFactor = targetMass_MeV / ( projectileMass_MeV + targetMass_MeV );
113 if( ( energyType == MCGIDI_energyType_primaryGamma ) || ( energyType == MCGIDI_energyType_discreteGamma ) ) {
114 energy->type = energyType;
115 energy->gammaEnergy_MeV = gammaEnergy_MeV;
116 energy->frame = xDataTOM_frame_lab;
117 if( energyType == MCGIDI_energyType_primaryGamma ) energy->primaryGammaMassFactor = energy->e_inCOMFactor; }
123 if( linearElement == NULL ) {
139 smr_setReportError2( smr, smr_unknownID, 1,
"unsupported energy type: nativeData = '%s'", nativeData );
142 frameElement = functional; }
144 char const *toUnits[3] = {
"MeV",
"MeV",
"1/MeV" };
146 frameElement = linearElement;
148 energy->type = MCGIDI_energyType_linear;
152 distribution->energy =
energy;
166 xDataTOM_element *child;
169 if( strcmp( child->name,
"weighted" ) )
goto err;
171 energy->weightedFunctionals.numberOfWeights++;
173 energy->type = MCGIDI_energyType_weightedFunctional;
184 xDataTOM_element *child;
185 MCGIDI_energy *
energy = NULL;
186 ptwXYPoints *weight = NULL;
187 char const *toUnits[2] = {
"MeV",
"" };
191 if( strcmp( child->name,
"weight" ) == 0 ) {
193 else if( strcmp( child->name,
"evaporation" ) == 0 ) {
196 smr_setReportError2( smr, smr_unknownID, 1,
"unsupported energy type = '%s' in weighted functional", child->name );
200 weightedFunctional->weight = weight;
201 weightedFunctional->energy =
energy;
215 xDataTOM_element *thetaTOM, *gTOM;
216 ptwXYPoints *theta = NULL, *
g = NULL;
217 char const *toUnits[2] = {
"MeV",
"MeV" };
229 if( std::fabs( 1. - norm ) > 0.001 ) printf(
"bad norm = %e\n", norm );
231 energy->type = MCGIDI_energyType_generalEvaporation;
232 energy->theta = theta;
245 char const *U, *toUnits[2] = {
"MeV",
"MeV" };
246 xDataTOM_element *thetaTOM;
249 smr_setReportError2( smr, smr_unknownID, 1,
"functional form '%s' missing 'U' attribute", element->name );
255 energy->type = MCGIDI_energyType_simpleMaxwellianFission;
266 char const *U, *toUnits[2] = {
"MeV",
"MeV" };
267 xDataTOM_element *thetaTOM;
270 smr_setReportError2( smr, smr_unknownID, 1,
"functional form '%s' missing 'U' attribute", element->name );
276 energy->type = MCGIDI_energyType_evaporation;
287 char const *U, *toUnits[2] = {
"MeV",
"MeV" };
288 xDataTOM_element *aOrBTOM;
291 smr_setReportError2( smr, smr_unknownID, 1,
"functional form '%s' missing 'U' attribute", element->name );
299 toUnits[1] =
"1/MeV";
303 energy->type = MCGIDI_energyType_Watt;
314 int iE, length, nXs, i1,
n;
315 double E, T_M, EFL, EFH, argList[3], xs[] = { 1e-5, 1e-3, 1e-1, 1
e1, 1
e3, 1e5, 3e7 }, norm;
316 ptwXYPoints *ptwXY_TM = NULL, *pdfXY = NULL;
318 ptwXPoints *cdfX = NULL;
319 nfu_status status = nfu_Okay;
320 xDataTOM_element *TM_TOM;
322 MCGIDI_pdfsOfXGivenW *dists = &(energy->dists);
324 char const *EF, *TMUnits[2] = {
"MeV",
"MeV" };
326 nXs =
sizeof( xs ) /
sizeof( xs[0] );
329 smr_setReportError2( smr, smr_unknownID, 1,
"MadlandNix '%s' missing 'EFL' attribute", functional->name );
336 smr_setReportError2( smr, smr_unknownID, 1,
"MadlandNix '%s' missing 'EFH' attribute", functional->name );
347 dists->interpolationWY = ptwXY_interpolationLinLin;
348 dists->interpolationXY = ptwXY_interpolationLinLin;
349 if( ( dists->Ws = (
double *) smr_malloc2( smr, length *
sizeof(
double ), 1,
"dists->Ws" ) ) == NULL )
goto err;
350 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, length *
sizeof( MCGIDI_pdfOfX ), 0,
"dists->dist" ) ) == NULL )
goto err;
352 for( iE = 0; iE < length; iE++ ) {
355 dist = &(dists->dist[iE]);
359 (
void *) argList, 1e-3, 0, 12, &status ) ) == NULL )
goto err;
361 smr_setReportError2( smr, smr_unknownID, 1,
"ptwXY_normalize err = %d: %s\n", status,
nfu_statusMessage( status ) );
368 if( ( dist->Xs = (
double *) smr_malloc2( smr, 3 * n *
sizeof(
double ), 0,
"dist->Xs" ) ) == NULL )
goto err;
370 dist->pdf = &(dist->Xs[
n]);
371 dist->cdf = &(dist->pdf[
n]);
373 for( i1 = 0; i1 <
n; i1++ ) {
375 dist->Xs[i1] = point->x;
376 dist->pdf[i1] = point->y;
380 smr_setReportError2( smr, smr_unknownID, 1,
"ptwXY_runningIntegral err = %d: %s\n", status,
nfu_statusMessage( status ) );
386 for( i1 = 0; i1 <
n; i1++ ) dist->pdf[i1] /= norm;
391 energy->type = MCGIDI_energyType_MadlandNix;
397 if( ptwXY_TM != NULL )
ptwXY_free( ptwXY_TM );
399 if( cdfX != NULL ) cdfX =
ptwX_free( cdfX );
408 double *parameters = (
double *) argList, EFL, EFH, T_M;
409 nfu_status status = nfu_Okay;
424 double u1, u2, E1, E2, gamma1, gamma2, signG = 1;
426 u1 = std::sqrt( Ep ) - std::sqrt( E_F );
428 u2 = std::sqrt( Ep ) + std::sqrt( E_F );
433 if( *status != nfu_Okay )
return( 0. );
442 if( *status != nfu_Okay )
return( 0. );
443 return( ( u2 * std::sqrt( u2 ) * E2 - u1 * std::sqrt( u1 ) * E1 + signG * ( gamma2 - gamma1 ) ) / ( 3 * std::sqrt( E_F * T_M ) ) );
449 MCGIDI_distribution *distribution ) {
452 double xs[2] = { 0.0, 1.0 }, productMass_MeV, norm;
453 ptwXYPoints *pdf = NULL;
459 smr_setReportError2( smr, smr_unknownID, 1,
"functional form '%s' missing 'mass' attribute", functional->name );
463 argList[0] = energy->NBodyPhaseSpace.numberOfProducts;
465 smr_setReportError2( smr, smr_unknownID, 1,
"creating NBodyPhaseSpace pdf failed with ptwXY_createFromFunction error = %d (%s)",
472 energy->NBodyPhaseSpace.massFactor = ( 1. - productMass_MeV / ( MCGIDI_AMU2MeV * energy->NBodyPhaseSpace.mass ) );
477 energy->type = MCGIDI_energyType_NBodyPhaseSpace;
490 int numberOfProducts = *((
int *) argList);
491 double e = 0.5 * ( 3 * numberOfProducts - 8 );
500 MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
504 double theta, randomEp, Watt_a, Watt_b, e_in = modes.getProjectileEnergy( );
505 MCGIDI_pdfsOfXGivenW_sampled sampled;
507 decaySamplingInfo->frame = energy->frame;
508 switch( energy->type ) {
509 case MCGIDI_energyType_primaryGamma :
510 decaySamplingInfo->Ep = energy->gammaEnergy_MeV + e_in * energy->primaryGammaMassFactor;
512 case MCGIDI_energyType_discreteGamma :
513 decaySamplingInfo->Ep = energy->gammaEnergy_MeV;
515 case MCGIDI_energyType_linear :
516 randomEp = decaySamplingInfo->rng( decaySamplingInfo->rngState );
520 decaySamplingInfo->Ep = sampled.x;
522 case MCGIDI_energyType_generalEvaporation :
523 sampled.interpolationXY = energy->gInterpolation;
526 decaySamplingInfo->Ep = theta * sampled.x;
528 case MCGIDI_energyType_simpleMaxwellianFission :
531 decaySamplingInfo->Ep *= theta;
533 case MCGIDI_energyType_evaporation :
536 decaySamplingInfo->Ep *= theta;
538 case MCGIDI_energyType_Watt :
543 case MCGIDI_energyType_MadlandNix :
545 decaySamplingInfo->Ep = sampled.x;
547 case MCGIDI_energyType_NBodyPhaseSpace :
549 decaySamplingInfo->Ep = ( energy->e_inCOMFactor * e_in + energy->NBodyPhaseSpace.Q_MeV ) * energy->NBodyPhaseSpace.massFactor * sampled.x;
551 case MCGIDI_energyType_weightedFunctional :
555 smr_setReportError2( smr, smr_unknownID, 1,
"energy type = %d not supported", energy->type );
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.;
568 sqrt_x = std::sqrt( a );
569 norm_a = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x *
G4Exp( -a );
570 b = norm_a * decaySamplingInfo->rng( decaySamplingInfo->rngState );
571 for( i1 = 0; i1 < 16; i1++ ) {
572 x = 0.5 * ( xMin + xMax );
573 sqrt_x = std::sqrt( x );
574 c = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x *
G4Exp( -x );
582 decaySamplingInfo->Ep =
x;
592 double a = e_in_U_theta, b, c,
x, norm_a, xMin = 0., xMax =
a;
594 norm_a = 1 - ( 1 +
a ) *
G4Exp( -a );
595 b = 1. - norm_a * decaySamplingInfo->rng( decaySamplingInfo->rngState );
596 for( i1 = 0; i1 < 16; i1++ ) {
597 x = 0.5 * ( xMin + xMax );
598 c = ( 1 +
x ) *
G4Exp( -x );
606 decaySamplingInfo->Ep =
x;
613 static int MCGIDI_energy_sampleWatt( statusMessageReporting * ,
double e_in_U,
double Watt_a,
double Watt_b, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
617 double WattMin = 0., WattMax = e_in_U,
x, y,
z, energyOut, rand1, rand2;
619 x = 1. + ( Watt_b / ( 8. * Watt_a ) );
620 y = ( x + std::sqrt( x * x - 1. ) ) / Watt_a;
623 G4int icounter_max=1024;
627 if ( icounter > icounter_max ) {
628 G4cout <<
"Loop-counter exceeded the threshold value at " << __LINE__ <<
"th line of " << __FILE__ <<
"." <<
G4endl;
631 rand1 = -
G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
632 rand2 = -
G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
633 energyOut = y * rand1;
635 while( ( ( rand2 - z * ( rand1 + 1. ) ) * ( rand2 - z * ( rand1 + 1. ) ) > Watt_b * y * rand1 ) || ( energyOut < WattMin ) || ( energyOut > WattMax ) );
636 decaySamplingInfo->Ep = energyOut;
644 MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
649 double rW = decaySamplingInfo->rng( decaySamplingInfo->rngState ), cumulativeW = 0., weight;
650 MCGIDI_energyWeightedFunctional *weightedFunctional = NULL;
652 for( iW = 0; iW < energy->weightedFunctionals.numberOfWeights; iW++ ) {
653 weightedFunctional = &(energy->weightedFunctionals.weightedFunctional[iW]);
655 cumulativeW += weight;
656 if( cumulativeW >= rW )
break;
661 #if defined __cplusplus
xDataTOM_element * xDataTOME_getOneElementByName(statusMessageReporting *smr, xDataTOM_element *element, char const *name, int required)
static G4Pow * GetInstance()
int MCGIDI_energy_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution, ptwXYPoints *norms, enum MCGIDI_energyType energyType, double gammaEnergy_MeV)
G4double powA(G4double A, G4double y) const
ptwXY_interpolation ptwXY_getInterpolation(ptwXYPoints *ptwXY)
int64_t ptwXY_length(ptwXYPoints *ptwXY)
nfu_status ptwXY_normalize(ptwXYPoints *ptwXY)
MCGIDI_energy * MCGIDI_energy_free(statusMessageReporting *smr, MCGIDI_energy *energy)
xDataTOM_element * xDataTOME_getFirstElement(xDataTOM_element *element)
int MCGIDI_energy_initialize(statusMessageReporting *, MCGIDI_energy *energy)
static int MCGIDI_energy_parseNBodyPhaseSpaceFromTOM(statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy, MCGIDI_distribution *distribution)
double MCGIDI_outputChannel_getQ_MeV(statusMessageReporting *, MCGIDI_outputChannel *outputChannel, double)
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
static int MCGIDI_energy_parseWattFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
double MCGIDI_product_getTargetMass_MeV(statusMessageReporting *smr, MCGIDI_product *product)
int MCGIDI_sampling_pdfsOfX_release(statusMessageReporting *, MCGIDI_pdfOfX *dist)
int MCGIDI_sampling_pdfsOfXGivenW_release(statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *dists)
static double MCGIDI_energy_parseMadlandNixFromTOM_callback_g(double Ep, double EFL, double T_M, nfu_status *status)
ptwXYPoints * MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf(statusMessageReporting *smr, xDataTOM_XYs *XYs, ptwXY_interpolation interpolation, char const *toUnits[2])
static nfu_status MCGIDI_energy_NBodyPhaseSpacePDF_callback(double x, double *y, void *argList)
double nf_incompleteGammaFunctionComplementary(double a, double x, nfu_status *status)
double MCGIDI_product_getMass_MeV(statusMessageReporting *, MCGIDI_product *product)
static int MCGIDI_energy_parseMadlandNixFromTOM(statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy)
double ptwX_getPointAtIndex_Unsafely(ptwXPoints *ptwX, int64_t index)
int MCGIDI_energy_release(statusMessageReporting *smr, MCGIDI_energy *energy)
int xDataTOME_convertAttributeToInteger(statusMessageReporting *smr, xDataTOM_element *element, char const *name, int *n)
static nfu_status MCGIDI_energy_parseMadlandNixFromTOM_callback(double x, double *y, void *argList)
G4GLOB_DLL std::ostream G4cout
int MCGIDI_fromTOM_pdfOfX(statusMessageReporting *smr, ptwXYPoints *pdfXY, MCGIDI_pdfOfX *dist, double *norm)
int MCGIDI_fromTOM_pdfsOfXGivenW(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_pdfsOfXGivenW *dists, ptwXYPoints *norms, char const *toUnits[3])
int MCGIDI_misc_PQUStringToDouble(statusMessageReporting *smr, char const *str, char const *unit, double conversion, double *value)
static int MCGIDI_energy_parseGeneralEvaporationFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy)
int smr_isOk(statusMessageReporting *smr)
void * xDataTOME_getXDataIfID(statusMessageReporting *smr, xDataTOM_element *TE, char const *ID)
void * smr_freeMemory(void **p)
static int MCGIDI_energy_sampleWeightedFunctional(statusMessageReporting *smr, MCGIDI_energy *energy, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
double MCGIDI_product_getProjectileMass_MeV(statusMessageReporting *smr, MCGIDI_product *product)
double nf_incompleteGammaFunction(double a, double x, nfu_status *status)
int MCGIDI_misc_PQUStringToDoubleInUnitOf(statusMessageReporting *smr, char const *str, char const *toUnit, double *value)
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
nfu_status ptwXY_getXYPairAtIndex(ptwXYPoints *ptwXY, int64_t index, double *x, double *y)
char const * xDataTOM_getAttributesValueInElement(xDataTOM_element *element, char const *name)
static int MCGIDI_energy_parseSimpleMaxwellianFissionFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy)
static int MCGIDI_energy_parseEvaporationFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
int MCGIDI_energy_sampleEnergy(statusMessageReporting *smr, MCGIDI_energy *energy, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
static int MCGIDI_energy_sampleEvaporation(statusMessageReporting *smr, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo)
xDataTOM_element * xDataTOME_getNextElement(xDataTOM_element *element)
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
G4double energy(const ThreeVector &p, const G4double m)
static int MCGIDI_energy_sampleWatt(statusMessageReporting *smr, double e_in_U, double Watt_a, double Watt_b, MCGIDI_decaySamplingInfo *decaySamplingInfo)
const G4double x[NPOINTSGL]
const char * nfu_statusMessage(nfu_status status)
int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW(MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue)
double MCGIDI_sampling_ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x1)
ptwXPoints * ptwXY_runningIntegral(ptwXYPoints *ptwXY, nfu_status *status)
static int MCGIDI_energy_sampleSimpleMaxwellianFission(statusMessageReporting *smr, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo)
ptwXYPoints * MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf(statusMessageReporting *smr, xDataTOM_element *linear, char const *toUnits[2])
ptwXYPoints * ptwXY_createFromFunction(int n, double *xs, ptwXY_createFromFunction_callback func, void *argList, double accuracy, int checkForRoots, int biSectionMax, nfu_status *status)
static int MCGIDI_energy_parseWeightedFunctionalsFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy)
MCGIDI_energy * MCGIDI_energy_new(statusMessageReporting *smr)
int MCGIDI_sampling_sampleX_from_pdfOfX(MCGIDI_pdfOfX *dist, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue)
enum xDataTOM_frame MCGIDI_misc_getProductFrame(statusMessageReporting *smr, xDataTOM_element *frameElement)
static int MCGIDI_energy_parseWeightFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energyWeightedFunctional *weightedFunctional)
double nf_exponentialIntegral(int n, double x, nfu_status *status)