8 #include "MCGIDI_fromTOM.h"
9 #include "MCGIDI_misc.h"
10 #include "MCGIDI_private.h"
12 #if defined __cplusplus
22 MCGIDI_angular *angular;
24 if( ( angular = (MCGIDI_angular *) smr_malloc2( smr,
sizeof( MCGIDI_angular ), 0,
"angular" ) ) == NULL )
return( NULL );
33 memset( angular, 0,
sizeof( MCGIDI_angular ) );
60 double productMass_MeV,
double residualMass_MeV ) {
62 if( angular == NULL )
return( 0 );
63 angular->projectileMass_MeV = projectileMass_MeV;
64 angular->targetMass_MeV = targetMass_MeV;
65 angular->productMass_MeV = productMass_MeV;
66 angular->residualMass_MeV = residualMass_MeV;
72 int MCGIDI_angular_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution, ptwXYPoints *norms ) {
74 MCGIDI_angular *angular = NULL;
75 xDataTOM_element *angularElement, *linearElement, *frameElement = NULL;
76 char const *nativeData;
77 ptwXYPoints *pdfXY = NULL;
78 ptwXPoints *cdfX = NULL;
79 ptwXY_interpolation interpolationXY, interpolationWY;
85 if( strcmp( nativeData,
"isotropic" ) == 0 ) {
87 smr_setReportError2( smr, smr_unknownID, 1,
"angular type missing for nativeData = '%s'", nativeData );
90 angular->type = MCGIDI_angularType_isotropic; }
91 else if( strcmp( nativeData,
"recoil" ) == 0 ) {
92 angular->type = MCGIDI_angularType_recoil; }
95 double norm, energyFactor;
98 xDataTOM_W_XYs *W_XYs;
100 MCGIDI_pdfsOfXGivenW *dists = &(angular->dists);
102 char const *energyUnit, *multiplicityProbabilityUnits[2] = {
"",
"" };
106 smr_setReportError2( smr, smr_unknownID, 1,
"unsupported angular type: nativeData = '%s'", nativeData );
110 frameElement = linearElement;
114 dists->interpolationWY = interpolationWY;
115 dists->interpolationXY = interpolationXY;
117 if( ( W_XYs = (xDataTOM_W_XYs *)
xDataTOME_getXDataIfID( smr, linearElement,
"W_XYs" ) ) == NULL )
goto err;
118 if( ( dists->Ws = (
double *) smr_malloc2( smr, W_XYs->length *
sizeof(
double ), 1,
"dists->Ws" ) ) == NULL )
goto err;
119 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, W_XYs->length *
sizeof( MCGIDI_pdfOfX ), 0,
"dists->dist" ) ) == NULL )
goto err;
126 for( i = 0; i < W_XYs->length; i++ ) {
127 XYs = &(W_XYs->XYs[i]);
128 dist = &(dists->dist[i]);
129 dists->Ws[i] = XYs->value * energyFactor;
134 if( ( dist->Xs = (
double *) smr_malloc2( smr, 3 * n *
sizeof(
double ), 0,
"dist->Xs" ) ) == NULL )
goto err;
136 dist->pdf = &(dist->Xs[
n]);
137 dist->cdf = &(dist->pdf[
n]);
139 for( j = 0; j <
n; j++ ) {
141 dist->Xs[j] = point->x;
142 dist->pdf[j] = point->y;
146 smr_setReportError2( smr, smr_unknownID, 1,
"ptwXY_runningIntegral err = %d: %s\n", status,
nfu_statusMessage( status ) );
150 if( norms != NULL ) {
152 else if( std::fabs( 1. - norm ) > 0.99 ) {
153 smr_setReportError2( smr, smr_unknownID, 1,
"bad norm = %e for angular.linear data", norm );
157 for( j = 0; j <
n; j++ ) dist->pdf[j] /= norm;
161 angular->type = MCGIDI_angularType_linear;
164 if( frameElement != NULL ) {
168 distribution->angular = angular;
169 distribution->type = MCGIDI_distributionType_angular_e;
175 if( cdfX != NULL ) cdfX =
ptwX_free( cdfX );
183 MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
185 double randomMu = decaySamplingInfo->rng( decaySamplingInfo->rngState );
186 MCGIDI_pdfsOfXGivenW_sampled sampled;
188 switch( angular->type ) {
189 case MCGIDI_angularType_isotropic :
190 decaySamplingInfo->frame = angular->frame;
191 decaySamplingInfo->mu = 1. - 2. * decaySamplingInfo->rng( decaySamplingInfo->rngState );
193 case MCGIDI_angularType_linear :
194 decaySamplingInfo->frame = angular->frame;
196 sampled.w = modes.getProjectileEnergy( );
198 decaySamplingInfo->mu = sampled.x;
200 case MCGIDI_angularType_recoil :
202 smr_setReportError2( smr, smr_unknownID, 1,
"angular type = %d not supported", angular->type );
207 #if defined __cplusplus
xDataTOM_element * xDataTOME_getOneElementByName(statusMessageReporting *smr, xDataTOM_element *element, char const *name, int required)
int MCGIDI_angular_sampleMu(statusMessageReporting *smr, MCGIDI_angular *angular, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
int64_t ptwXY_length(ptwXYPoints *ptwXY)
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
int MCGIDI_sampling_pdfsOfXGivenW_release(statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *dists)
ptwXYPoints * MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf(statusMessageReporting *smr, xDataTOM_XYs *XYs, ptwXY_interpolation interpolation, char const *toUnits[2])
double ptwX_getPointAtIndex_Unsafely(ptwXPoints *ptwX, int64_t index)
char const * xDataTOM_subAxes_getUnit(statusMessageReporting *smr, xDataTOM_subAxes *subAxes, int index)
int smr_isOk(statusMessageReporting *smr)
void * xDataTOME_getXDataIfID(statusMessageReporting *smr, xDataTOM_element *TE, char const *ID)
void * smr_freeMemory(void **p)
double MCGIDI_misc_getUnitConversionFactor(statusMessageReporting *smr, char const *fromUnit, char const *toUnit)
int MCGIDI_angular_setTwoBodyMasses(statusMessageReporting *, MCGIDI_angular *angular, double projectileMass_MeV, double targetMass_MeV, double productMass_MeV, double residualMass_MeV)
int MCGIDI_angular_release(statusMessageReporting *smr, MCGIDI_angular *angular)
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
char const * xDataTOM_getAttributesValueInElement(xDataTOM_element *element, char const *name)
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
MCGIDI_angular * MCGIDI_angular_free(statusMessageReporting *smr, MCGIDI_angular *angular)
const char * nfu_statusMessage(nfu_status status)
int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW(MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue)
ptwXPoints * ptwXY_runningIntegral(ptwXYPoints *ptwXY, nfu_status *status)
int MCGIDI_angular_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution, ptwXYPoints *norms)
enum xDataTOM_frame MCGIDI_misc_getProductFrame(statusMessageReporting *smr, xDataTOM_element *frameElement)
MCGIDI_angular * MCGIDI_angular_new(statusMessageReporting *smr)
int MCGIDI_angular_initialize(statusMessageReporting *, MCGIDI_angular *angular)
int MCGIDI_fromTOM_interpolation(statusMessageReporting *smr, xDataTOM_element *element, int index, ptwXY_interpolation *interpolation)