Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
MCGIDI_angular.cc File Reference
#include <string.h>
#include <cmath>
#include "MCGIDI_fromTOM.h"
#include "MCGIDI_misc.h"
#include "MCGIDI_private.h"

Go to the source code of this file.

Functions

MCGIDI_angularMCGIDI_angular_new (statusMessageReporting *smr)
 
int MCGIDI_angular_initialize (statusMessageReporting *, MCGIDI_angular *angular)
 
MCGIDI_angularMCGIDI_angular_free (statusMessageReporting *smr, MCGIDI_angular *angular)
 
int MCGIDI_angular_release (statusMessageReporting *smr, MCGIDI_angular *angular)
 
int MCGIDI_angular_setTwoBodyMasses (statusMessageReporting *, MCGIDI_angular *angular, double projectileMass_MeV, double targetMass_MeV, double productMass_MeV, double residualMass_MeV)
 
int MCGIDI_angular_parseFromTOM (statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution, ptwXYPoints *norms)
 
int MCGIDI_angular_sampleMu (statusMessageReporting *smr, MCGIDI_angular *angular, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
 

Function Documentation

◆ MCGIDI_angular_free()

MCGIDI_angular * MCGIDI_angular_free ( statusMessageReporting smr,
MCGIDI_angular angular 
)

Definition at line 39 of file MCGIDI_angular.cc.

39 {
40
41 MCGIDI_angular_release( smr, angular );
42 smr_freeMemory( (void **) &angular );
43 return( NULL );
44}
int MCGIDI_angular_release(statusMessageReporting *smr, MCGIDI_angular *angular)
void * smr_freeMemory(void **p)

Referenced by MCGIDI_angular_new(), MCGIDI_angular_parseFromTOM(), MCGIDI_distribution_release(), and MCGIDI_LLNLAngular_angularEnergy_parseFromTOM().

◆ MCGIDI_angular_initialize()

int MCGIDI_angular_initialize ( statusMessageReporting smr,
MCGIDI_angular angular 
)

Definition at line 31 of file MCGIDI_angular.cc.

31 {
32
33 memset( angular, 0, sizeof( MCGIDI_angular ) );
34 return( 0 );
35}

Referenced by MCGIDI_angular_new(), and MCGIDI_angular_release().

◆ MCGIDI_angular_new()

MCGIDI_angular * MCGIDI_angular_new ( statusMessageReporting smr)

Definition at line 20 of file MCGIDI_angular.cc.

20 {
21
22 MCGIDI_angular *angular;
23
24 if( ( angular = (MCGIDI_angular *) smr_malloc2( smr, sizeof( MCGIDI_angular ), 0, "angular" ) ) == NULL ) return( NULL );
25 if( MCGIDI_angular_initialize( smr, angular ) ) angular = MCGIDI_angular_free( smr, angular );
26 return( angular );
27}
int MCGIDI_angular_initialize(statusMessageReporting *, MCGIDI_angular *angular)
MCGIDI_angular * MCGIDI_angular_free(statusMessageReporting *smr, MCGIDI_angular *angular)
#define smr_malloc2(smr, size, zero, forItem)

Referenced by MCGIDI_angular_parseFromTOM().

◆ MCGIDI_angular_parseFromTOM()

int MCGIDI_angular_parseFromTOM ( statusMessageReporting smr,
xDataTOM_element element,
MCGIDI_distribution distribution,
ptwXYPoints norms 
)

Definition at line 72 of file MCGIDI_angular.cc.

72 {
73
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;
80
81 if( ( angularElement = xDataTOME_getOneElementByName( smr, element, "angular", 1 ) ) == NULL ) goto err;
82 if( ( angular = MCGIDI_angular_new( smr ) ) == NULL ) goto err;
83
84 if( ( nativeData = xDataTOM_getAttributesValueInElement( angularElement, "nativeData" ) ) == NULL ) goto err;
85 if( strcmp( nativeData, "isotropic" ) == 0 ) {
86 if( ( frameElement = xDataTOME_getOneElementByName( smr, angularElement, "isotropic", 1 ) ) == NULL ) {
87 smr_setReportError2( smr, smr_unknownID, 1, "angular type missing for nativeData = '%s'", nativeData );
88 goto err;
89 }
91 else if( strcmp( nativeData, "recoil" ) == 0 ) { /* BRB. Needs work to get referenced product data?????? */
93 else {
94 int i, j, n;
95 double norm, energyFactor;
96 nfu_status status;
97 xDataTOM_XYs *XYs;
98 xDataTOM_W_XYs *W_XYs;
99 ptwXYPoint *point;
100 MCGIDI_pdfsOfXGivenW *dists = &(angular->dists);
101 MCGIDI_pdfOfX *dist;
102 char const *energyUnit, *multiplicityProbabilityUnits[2] = { "", "" };
103
104 if( ( linearElement = xDataTOME_getOneElementByName( NULL, angularElement, "linear", 0 ) ) == NULL ) {
105 if( ( linearElement = xDataTOME_getOneElementByName( smr, angularElement, "pointwise", 1 ) ) == NULL ) {
106 smr_setReportError2( smr, smr_unknownID, 1, "unsupported angular type: nativeData = '%s'", nativeData );
107 goto err;
108 }
109 }
110 frameElement = linearElement;
111
112 if( MCGIDI_fromTOM_interpolation( smr, linearElement, 0, &interpolationWY ) ) goto err;
113 if( MCGIDI_fromTOM_interpolation( smr, linearElement, 1, &interpolationXY ) ) goto err;
114 dists->interpolationWY = interpolationWY;
115 dists->interpolationXY = interpolationXY;
116
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;
120
121 energyUnit = xDataTOM_subAxes_getUnit( smr, &(W_XYs->subAxes), 0 );
122 if( !smr_isOk( smr ) ) goto err;
123 energyFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyUnit, "MeV" );
124 if( !smr_isOk( smr ) ) goto err;
125
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;
130 if( ( pdfXY = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, interpolationXY, multiplicityProbabilityUnits ) ) == NULL ) goto err;
131 if( ptwXY_simpleCoalescePoints( pdfXY ) != nfu_Okay ) goto err;
132 dist->numberOfXs = n = (int) ptwXY_length( pdfXY );
133
134 if( ( dist->Xs = (double *) smr_malloc2( smr, 3 * n * sizeof( double ), 0, "dist->Xs" ) ) == NULL ) goto err;
135 dists->numberOfWs++;
136 dist->pdf = &(dist->Xs[n]);
137 dist->cdf = &(dist->pdf[n]);
138
139 for( j = 0; j < n; j++ ) {
140 point = ptwXY_getPointAtIndex_Unsafely( pdfXY, j );
141 dist->Xs[j] = point->x;
142 dist->pdf[j] = point->y;
143 }
144
145 if( ( cdfX = ptwXY_runningIntegral( pdfXY, &status ) ) == NULL ) {
146 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_runningIntegral err = %d: %s\n", status, nfu_statusMessage( status ) );
147 goto err;
148 }
149 norm = ptwX_getPointAtIndex_Unsafely( cdfX, n - 1 );
150 if( norms != NULL ) {
151 ptwXY_setValueAtX( norms, XYs->value, norm ); }
152 else if( std::fabs( 1. - norm ) > 0.99 ) {
153 smr_setReportError2( smr, smr_unknownID, 1, "bad norm = %e for angular.linear data", norm );
154 goto err;
155 }
156 for( j = 0; j < n; j++ ) dist->cdf[j] = ptwX_getPointAtIndex_Unsafely( cdfX, j ) / norm;
157 for( j = 0; j < n; j++ ) dist->pdf[j] /= norm;
158 pdfXY = ptwXY_free( pdfXY );
159 cdfX = ptwX_free( cdfX );
160 }
162 }
163
164 if( frameElement != NULL ) {
165 if( ( angular->frame = MCGIDI_misc_getProductFrame( smr, frameElement ) ) == xDataTOM_frame_invalid ) goto err;
166 }
167
168 distribution->angular = angular;
170
171 return( 0 );
172
173err:
174 if( pdfXY != NULL ) ptwXY_free( pdfXY );
175 if( cdfX != NULL ) cdfX = ptwX_free( cdfX );
176 if ( angular != NULL ) MCGIDI_angular_free( smr, angular );
177 return( 1 );
178}
@ MCGIDI_distributionType_angular_e
Definition: MCGIDI.h:208
@ MCGIDI_angularType_isotropic
Definition: MCGIDI.h:212
@ MCGIDI_angularType_recoil
Definition: MCGIDI.h:212
@ MCGIDI_angularType_linear
Definition: MCGIDI.h:212
MCGIDI_angular * MCGIDI_angular_new(statusMessageReporting *smr)
int MCGIDI_fromTOM_interpolation(statusMessageReporting *smr, xDataTOM_element *element, int index, enum ptwXY_interpolation_e *interpolation)
ptwXYPoints * MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf(statusMessageReporting *smr, xDataTOM_XYs *XYs, ptwXY_interpolation interpolation, char const *units[2])
Definition: MCGIDI_misc.cc:405
double MCGIDI_misc_getUnitConversionFactor(statusMessageReporting *smr, char const *fromUnit, char const *toUnit)
Definition: MCGIDI_misc.cc:381
enum xDataTOM_frame MCGIDI_misc_getProductFrame(statusMessageReporting *smr, xDataTOM_element *frameElement)
Definition: MCGIDI_misc.cc:315
@ nfu_Okay
Definition: nf_utilities.h:25
enum nfu_status_e nfu_status
const char * nfu_statusMessage(nfu_status status)
Definition: nf_utilities.cc:76
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
Definition: ptwXY_core.cc:876
nfu_status ptwXY_simpleCoalescePoints(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:529
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
Definition: ptwXY_core.cc:684
enum ptwXY_interpolation_e ptwXY_interpolation
int64_t ptwXY_length(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:583
ptwXPoints * ptwXY_runningIntegral(ptwXYPoints *ptwXY, nfu_status *status)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:574
double ptwX_getPointAtIndex_Unsafely(ptwXPoints *ptwX, int64_t index)
Definition: ptwX_core.cc:215
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
Definition: ptwX_core.cc:158
#define smr_setReportError2(smr, libraryID, code, fmt,...)
int smr_isOk(statusMessageReporting *smr)
#define smr_unknownID
MCGIDI_pdfsOfXGivenW dists
Definition: MCGIDI.h:322
enum xDataTOM_frame frame
Definition: MCGIDI.h:319
enum MCGIDI_angularType type
Definition: MCGIDI.h:320
MCGIDI_angular * angular
Definition: MCGIDI.h:383
enum MCGIDI_distributionType type
Definition: MCGIDI.h:382
double * Xs
Definition: MCGIDI.h:299
double * pdf
Definition: MCGIDI.h:300
int numberOfXs
Definition: MCGIDI.h:298
double * cdf
Definition: MCGIDI.h:301
ptwXY_interpolation interpolationXY
Definition: MCGIDI.h:306
MCGIDI_pdfOfX * dist
Definition: MCGIDI.h:308
ptwXY_interpolation interpolationWY
Definition: MCGIDI.h:306
double y
Definition: ptwXY.h:62
double x
Definition: ptwXY.h:62
xDataTOM_subAxes subAxes
Definition: xDataTOM.h:96
xDataTOM_XYs * XYs
Definition: xDataTOM.h:97
double value
Definition: xDataTOM.h:82
xDataTOM_element * xDataTOME_getOneElementByName(statusMessageReporting *smr, xDataTOM_element *element, char const *name, int required)
Definition: xDataTOM.cc:246
char const * xDataTOM_getAttributesValueInElement(xDataTOM_element *element, char const *name)
Definition: xDataTOM.cc:286
char const * xDataTOM_subAxes_getUnit(statusMessageReporting *smr, xDataTOM_subAxes *subAxes, int index)
@ xDataTOM_frame_invalid
Definition: xDataTOM.h:23
void * xDataTOME_getXDataIfID(statusMessageReporting *smr, xDataTOM_element *TE, char const *ID)
Definition: xDataTOM.cc:512

Referenced by MCGIDI_distribution_parseFromTOM(), MCGIDI_LLNLAngular_angularEnergy_parseFromTOM(), and MCGIDI_uncorrelated_parseFromTOM().

◆ MCGIDI_angular_release()

int MCGIDI_angular_release ( statusMessageReporting smr,
MCGIDI_angular angular 
)

Definition at line 48 of file MCGIDI_angular.cc.

48 {
49
50
52
53 MCGIDI_angular_initialize( smr, angular );
54 return( 0 );
55}
int MCGIDI_sampling_pdfsOfXGivenW_release(statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *dists)

Referenced by MCGIDI_angular_free().

◆ MCGIDI_angular_sampleMu()

int MCGIDI_angular_sampleMu ( statusMessageReporting smr,
MCGIDI_angular angular,
MCGIDI_quantitiesLookupModes modes,
MCGIDI_decaySamplingInfo decaySamplingInfo 
)

Definition at line 182 of file MCGIDI_angular.cc.

183 {
184
185 double randomMu = decaySamplingInfo->rng( decaySamplingInfo->rngState );
187
188 switch( angular->type ) {
190 decaySamplingInfo->frame = angular->frame;
191 decaySamplingInfo->mu = 1. - 2. * decaySamplingInfo->rng( decaySamplingInfo->rngState );
192 break;
194 decaySamplingInfo->frame = angular->frame;
195 sampled.smr = smr;
196 sampled.w = modes.getProjectileEnergy( );
197 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(angular->dists), &sampled, randomMu );
198 decaySamplingInfo->mu = sampled.x;
199 break;
201 default :
202 smr_setReportError2( smr, smr_unknownID, 1, "angular type = %d not supported", angular->type );
203 }
204 return( !smr_isOk( smr ) );
205}
int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW(MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double r)
double getProjectileEnergy(void) const
Definition: MCGIDI.h:97
double(* rng)(void *)
Definition: MCGIDI.h:258
enum xDataTOM_frame frame
Definition: MCGIDI.h:256
statusMessageReporting * smr
Definition: MCGIDI.h:312

Referenced by MCGIDI_outputChannel_sampleProductsAtE(), MCGIDI_product_sampleMu(), and MCGIDI_uncorrelated_sampleDistribution().

◆ MCGIDI_angular_setTwoBodyMasses()

int MCGIDI_angular_setTwoBodyMasses ( statusMessageReporting smr,
MCGIDI_angular angular,
double  projectileMass_MeV,
double  targetMass_MeV,
double  productMass_MeV,
double  residualMass_MeV 
)

Definition at line 59 of file MCGIDI_angular.cc.

60 {
61
62 if( angular == NULL ) return( 0 ); /* ???????? This needs work. Happens when first product of a two-body reaction as no distribution. */
63 angular->projectileMass_MeV = projectileMass_MeV;
64 angular->targetMass_MeV = targetMass_MeV;
65 angular->productMass_MeV = productMass_MeV;
66 angular->residualMass_MeV = residualMass_MeV;
67 return( 0 );
68}
double residualMass_MeV
Definition: MCGIDI.h:323
double productMass_MeV
Definition: MCGIDI.h:323
double projectileMass_MeV
Definition: MCGIDI.h:323
double targetMass_MeV
Definition: MCGIDI.h:323

Referenced by MCGIDI_product_setTwoBodyMasses().