| Geant4 Cross Reference |
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 #include <string.h>
6 #include <cmath>
7
8 #include "MCGIDI_fromTOM.h"
9 #include "MCGIDI_misc.h"
10 #include "MCGIDI_private.h"
11
12 #if defined __cplusplus
13 namespace GIDI {
14 using namespace GIDI;
15 #endif
16
17 /*
18 ************************************************************
19 */
20 MCGIDI_angular *MCGIDI_angular_new( statusMessageReporting *smr ) {
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 }
28 /*
29 ************************************************************
30 */
31 int MCGIDI_angular_initialize( statusMessageReporting * /*smr*/, MCGIDI_angular *angular ) {
32
33 memset( angular, 0, sizeof( MCGIDI_angular ) );
34 return( 0 );
35 }
36 /*
37 ************************************************************
38 */
39 MCGIDI_angular *MCGIDI_angular_free( statusMessageReporting *smr, MCGIDI_angular *angular ) {
40
41 MCGIDI_angular_release( smr, angular );
42 smr_freeMemory( (void **) &angular );
43 return( NULL );
44 }
45 /*
46 ************************************************************
47 */
48 int MCGIDI_angular_release( statusMessageReporting *smr, MCGIDI_angular *angular ) {
49
50
51 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(angular->dists) );
52
53 MCGIDI_angular_initialize( smr, angular );
54 return( 0 );
55 }
56 /*
57 ************************************************************
58 */
59 int MCGIDI_angular_setTwoBodyMasses( statusMessageReporting * /*smr*/, MCGIDI_angular *angular, double projectileMass_MeV, double targetMass_MeV,
60 double productMass_MeV, double residualMass_MeV ) {
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 }
69 /*
70 ************************************************************
71 */
72 int MCGIDI_angular_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution, ptwXYPoints *norms ) {
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 }
90 angular->type = MCGIDI_angularType_isotropic; }
91 else if( strcmp( nativeData, "recoil" ) == 0 ) { /* BRB. Needs work to get referenced product data?????? */
92 angular->type = MCGIDI_angularType_recoil; }
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 }
161 angular->type = MCGIDI_angularType_linear;
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;
169 distribution->type = MCGIDI_distributionType_angular_e;
170
171 return( 0 );
172
173 err:
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 }
179 /*
180 ************************************************************
181 */
182 int MCGIDI_angular_sampleMu( statusMessageReporting *smr, MCGIDI_angular *angular, MCGIDI_quantitiesLookupModes &modes,
183 MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
184
185 double randomMu = decaySamplingInfo->rng( decaySamplingInfo->rngState );
186 MCGIDI_pdfsOfXGivenW_sampled sampled;
187
188 switch( angular->type ) {
189 case MCGIDI_angularType_isotropic :
190 decaySamplingInfo->frame = angular->frame;
191 decaySamplingInfo->mu = 1. - 2. * decaySamplingInfo->rng( decaySamplingInfo->rngState );
192 break;
193 case MCGIDI_angularType_linear :
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;
200 case MCGIDI_angularType_recoil :
201 default :
202 smr_setReportError2( smr, smr_unknownID, 1, "angular type = %d not supported", angular->type );
203 }
204 return( !smr_isOk( smr ) );
205 }
206
207 #if defined __cplusplus
208 }
209 #endif
210
211