| Geant4 Cross Reference |
1 /* 1 /*
2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>>
3 # <<END-copyright>> 3 # <<END-copyright>>
4 */ 4 */
5 #include <string.h> 5 #include <string.h>
6 #include <cmath> 6 #include <cmath>
7 7
8 #include "MCGIDI_fromTOM.h" 8 #include "MCGIDI_fromTOM.h"
9 #include "MCGIDI_misc.h" 9 #include "MCGIDI_misc.h"
10 #include "MCGIDI_private.h" 10 #include "MCGIDI_private.h"
11 11
12 #if defined __cplusplus 12 #if defined __cplusplus
13 namespace GIDI { 13 namespace GIDI {
14 using namespace GIDI; 14 using namespace GIDI;
15 #endif 15 #endif
16 16
17 /* 17 /*
18 ********************************************** 18 ************************************************************
19 */ 19 */
20 MCGIDI_angular *MCGIDI_angular_new( statusMess 20 MCGIDI_angular *MCGIDI_angular_new( statusMessageReporting *smr ) {
21 21
22 MCGIDI_angular *angular; 22 MCGIDI_angular *angular;
23 23
24 if( ( angular = (MCGIDI_angular *) smr_mal 24 if( ( angular = (MCGIDI_angular *) smr_malloc2( smr, sizeof( MCGIDI_angular ), 0, "angular" ) ) == NULL ) return( NULL );
25 if( MCGIDI_angular_initialize( smr, angula 25 if( MCGIDI_angular_initialize( smr, angular ) ) angular = MCGIDI_angular_free( smr, angular );
26 return( angular ); 26 return( angular );
27 } 27 }
28 /* 28 /*
29 ********************************************** 29 ************************************************************
30 */ 30 */
31 int MCGIDI_angular_initialize( statusMessageRe 31 int MCGIDI_angular_initialize( statusMessageReporting * /*smr*/, MCGIDI_angular *angular ) {
32 32
33 memset( angular, 0, sizeof( MCGIDI_angular 33 memset( angular, 0, sizeof( MCGIDI_angular ) );
34 return( 0 ); 34 return( 0 );
35 } 35 }
36 /* 36 /*
37 ********************************************** 37 ************************************************************
38 */ 38 */
39 MCGIDI_angular *MCGIDI_angular_free( statusMes 39 MCGIDI_angular *MCGIDI_angular_free( statusMessageReporting *smr, MCGIDI_angular *angular ) {
40 40
41 MCGIDI_angular_release( smr, angular ); 41 MCGIDI_angular_release( smr, angular );
42 smr_freeMemory( (void **) &angular ); 42 smr_freeMemory( (void **) &angular );
43 return( NULL ); 43 return( NULL );
44 } 44 }
45 /* 45 /*
46 ********************************************** 46 ************************************************************
47 */ 47 */
48 int MCGIDI_angular_release( statusMessageRepor 48 int MCGIDI_angular_release( statusMessageReporting *smr, MCGIDI_angular *angular ) {
49 49
50 50
51 MCGIDI_sampling_pdfsOfXGivenW_release( smr 51 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(angular->dists) );
52 52
53 MCGIDI_angular_initialize( smr, angular ); 53 MCGIDI_angular_initialize( smr, angular );
54 return( 0 ); 54 return( 0 );
55 } 55 }
56 /* 56 /*
57 ********************************************** 57 ************************************************************
58 */ 58 */
59 int MCGIDI_angular_setTwoBodyMasses( statusMes 59 int MCGIDI_angular_setTwoBodyMasses( statusMessageReporting * /*smr*/, MCGIDI_angular *angular, double projectileMass_MeV, double targetMass_MeV,
60 double productMass_MeV, double residualMas 60 double productMass_MeV, double residualMass_MeV ) {
61 61
62 if( angular == NULL ) return( 0 ); 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 = projectileMa 63 angular->projectileMass_MeV = projectileMass_MeV;
64 angular->targetMass_MeV = targetMass_MeV; 64 angular->targetMass_MeV = targetMass_MeV;
65 angular->productMass_MeV = productMass_MeV 65 angular->productMass_MeV = productMass_MeV;
66 angular->residualMass_MeV = residualMass_M 66 angular->residualMass_MeV = residualMass_MeV;
67 return( 0 ); 67 return( 0 );
68 } 68 }
69 /* 69 /*
70 ********************************************** 70 ************************************************************
71 */ 71 */
72 int MCGIDI_angular_parseFromTOM( statusMessage 72 int MCGIDI_angular_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution, ptwXYPoints *norms ) {
73 73
74 MCGIDI_angular *angular = NULL; 74 MCGIDI_angular *angular = NULL;
75 xDataTOM_element *angularElement, *linearE 75 xDataTOM_element *angularElement, *linearElement, *frameElement = NULL;
76 char const *nativeData; 76 char const *nativeData;
77 ptwXYPoints *pdfXY = NULL; 77 ptwXYPoints *pdfXY = NULL;
78 ptwXPoints *cdfX = NULL; 78 ptwXPoints *cdfX = NULL;
79 ptwXY_interpolation interpolationXY, inter 79 ptwXY_interpolation interpolationXY, interpolationWY;
80 80
81 if( ( angularElement = xDataTOME_getOneEle 81 if( ( angularElement = xDataTOME_getOneElementByName( smr, element, "angular", 1 ) ) == NULL ) goto err;
82 if( ( angular = MCGIDI_angular_new( smr ) 82 if( ( angular = MCGIDI_angular_new( smr ) ) == NULL ) goto err;
83 83
84 if( ( nativeData = xDataTOM_getAttributesV 84 if( ( nativeData = xDataTOM_getAttributesValueInElement( angularElement, "nativeData" ) ) == NULL ) goto err;
85 if( strcmp( nativeData, "isotropic" ) == 0 85 if( strcmp( nativeData, "isotropic" ) == 0 ) {
86 if( ( frameElement = xDataTOME_getOneE 86 if( ( frameElement = xDataTOME_getOneElementByName( smr, angularElement, "isotropic", 1 ) ) == NULL ) {
87 smr_setReportError2( smr, smr_unkn 87 smr_setReportError2( smr, smr_unknownID, 1, "angular type missing for nativeData = '%s'", nativeData );
88 goto err; 88 goto err;
89 } 89 }
90 angular->type = MCGIDI_angularType_iso 90 angular->type = MCGIDI_angularType_isotropic; }
91 else if( strcmp( nativeData, "recoil" ) == 91 else if( strcmp( nativeData, "recoil" ) == 0 ) { /* BRB. Needs work to get referenced product data?????? */
92 angular->type = MCGIDI_angularType_rec 92 angular->type = MCGIDI_angularType_recoil; }
93 else { 93 else {
94 int i, j, n; 94 int i, j, n;
95 double norm, energyFactor; 95 double norm, energyFactor;
96 nfu_status status; 96 nfu_status status;
97 xDataTOM_XYs *XYs; 97 xDataTOM_XYs *XYs;
98 xDataTOM_W_XYs *W_XYs; 98 xDataTOM_W_XYs *W_XYs;
99 ptwXYPoint *point; 99 ptwXYPoint *point;
100 MCGIDI_pdfsOfXGivenW *dists = &(angula 100 MCGIDI_pdfsOfXGivenW *dists = &(angular->dists);
101 MCGIDI_pdfOfX *dist; 101 MCGIDI_pdfOfX *dist;
102 char const *energyUnit, *multiplicityP 102 char const *energyUnit, *multiplicityProbabilityUnits[2] = { "", "" };
103 103
104 if( ( linearElement = xDataTOME_getOne 104 if( ( linearElement = xDataTOME_getOneElementByName( NULL, angularElement, "linear", 0 ) ) == NULL ) {
105 if( ( linearElement = xDataTOME_ge 105 if( ( linearElement = xDataTOME_getOneElementByName( smr, angularElement, "pointwise", 1 ) ) == NULL ) {
106 smr_setReportError2( smr, smr_ 106 smr_setReportError2( smr, smr_unknownID, 1, "unsupported angular type: nativeData = '%s'", nativeData );
107 goto err; 107 goto err;
108 } 108 }
109 } 109 }
110 frameElement = linearElement; 110 frameElement = linearElement;
111 111
112 if( MCGIDI_fromTOM_interpolation( smr, 112 if( MCGIDI_fromTOM_interpolation( smr, linearElement, 0, &interpolationWY ) ) goto err;
113 if( MCGIDI_fromTOM_interpolation( smr, 113 if( MCGIDI_fromTOM_interpolation( smr, linearElement, 1, &interpolationXY ) ) goto err;
114 dists->interpolationWY = interpolation 114 dists->interpolationWY = interpolationWY;
115 dists->interpolationXY = interpolation 115 dists->interpolationXY = interpolationXY;
116 116
117 if( ( W_XYs = (xDataTOM_W_XYs *) xData 117 if( ( W_XYs = (xDataTOM_W_XYs *) xDataTOME_getXDataIfID( smr, linearElement, "W_XYs" ) ) == NULL ) goto err;
118 if( ( dists->Ws = (double *) smr_mallo 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 *) 119 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, W_XYs->length * sizeof( MCGIDI_pdfOfX ), 0, "dists->dist" ) ) == NULL ) goto err;
120 120
121 energyUnit = xDataTOM_subAxes_getUnit( 121 energyUnit = xDataTOM_subAxes_getUnit( smr, &(W_XYs->subAxes), 0 );
122 if( !smr_isOk( smr ) ) goto err; 122 if( !smr_isOk( smr ) ) goto err;
123 energyFactor = MCGIDI_misc_getUnitConv 123 energyFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyUnit, "MeV" );
124 if( !smr_isOk( smr ) ) goto err; 124 if( !smr_isOk( smr ) ) goto err;
125 125
126 for( i = 0; i < W_XYs->length; i++ ) { 126 for( i = 0; i < W_XYs->length; i++ ) {
127 XYs = &(W_XYs->XYs[i]); 127 XYs = &(W_XYs->XYs[i]);
128 dist = &(dists->dist[i]); 128 dist = &(dists->dist[i]);
129 dists->Ws[i] = XYs->value * energy 129 dists->Ws[i] = XYs->value * energyFactor;
130 if( ( pdfXY = MCGIDI_misc_dataFrom 130 if( ( pdfXY = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, interpolationXY, multiplicityProbabilityUnits ) ) == NULL ) goto err;
131 if( ptwXY_simpleCoalescePoints( pd 131 if( ptwXY_simpleCoalescePoints( pdfXY ) != nfu_Okay ) goto err;
132 dist->numberOfXs = n = (int) ptwXY 132 dist->numberOfXs = n = (int) ptwXY_length( pdfXY );
133 133
134 if( ( dist->Xs = (double *) smr_ma 134 if( ( dist->Xs = (double *) smr_malloc2( smr, 3 * n * sizeof( double ), 0, "dist->Xs" ) ) == NULL ) goto err;
135 dists->numberOfWs++; 135 dists->numberOfWs++;
136 dist->pdf = &(dist->Xs[n]); 136 dist->pdf = &(dist->Xs[n]);
137 dist->cdf = &(dist->pdf[n]); 137 dist->cdf = &(dist->pdf[n]);
138 138
139 for( j = 0; j < n; j++ ) { 139 for( j = 0; j < n; j++ ) {
140 point = ptwXY_getPointAtIndex_ 140 point = ptwXY_getPointAtIndex_Unsafely( pdfXY, j );
141 dist->Xs[j] = point->x; 141 dist->Xs[j] = point->x;
142 dist->pdf[j] = point->y; 142 dist->pdf[j] = point->y;
143 } 143 }
144 144
145 if( ( cdfX = ptwXY_runningIntegral 145 if( ( cdfX = ptwXY_runningIntegral( pdfXY, &status ) ) == NULL ) {
146 smr_setReportError2( smr, smr_ 146 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_runningIntegral err = %d: %s\n", status, nfu_statusMessage( status ) );
147 goto err; 147 goto err;
148 } 148 }
149 norm = ptwX_getPointAtIndex_Unsafe 149 norm = ptwX_getPointAtIndex_Unsafely( cdfX, n - 1 );
150 if( norms != NULL ) { 150 if( norms != NULL ) {
151 ptwXY_setValueAtX( norms, XYs- 151 ptwXY_setValueAtX( norms, XYs->value, norm ); }
152 else if( std::fabs( 1. - norm ) > 152 else if( std::fabs( 1. - norm ) > 0.99 ) {
153 smr_setReportError2( smr, smr_ 153 smr_setReportError2( smr, smr_unknownID, 1, "bad norm = %e for angular.linear data", norm );
154 goto err; 154 goto err;
155 } 155 }
156 for( j = 0; j < n; j++ ) dist->cdf 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 157 for( j = 0; j < n; j++ ) dist->pdf[j] /= norm;
158 pdfXY = ptwXY_free( pdfXY ); 158 pdfXY = ptwXY_free( pdfXY );
159 cdfX = ptwX_free( cdfX ); 159 cdfX = ptwX_free( cdfX );
160 } 160 }
161 angular->type = MCGIDI_angularType_lin 161 angular->type = MCGIDI_angularType_linear;
162 } 162 }
163 163
164 if( frameElement != NULL ) { 164 if( frameElement != NULL ) {
165 if( ( angular->frame = MCGIDI_misc_get 165 if( ( angular->frame = MCGIDI_misc_getProductFrame( smr, frameElement ) ) == xDataTOM_frame_invalid ) goto err;
166 } 166 }
167 167
168 distribution->angular = angular; 168 distribution->angular = angular;
169 distribution->type = MCGIDI_distributionTy 169 distribution->type = MCGIDI_distributionType_angular_e;
170 170
171 return( 0 ); 171 return( 0 );
172 172
173 err: 173 err:
174 if( pdfXY != NULL ) ptwXY_free( pdfXY ); 174 if( pdfXY != NULL ) ptwXY_free( pdfXY );
175 if( cdfX != NULL ) cdfX = ptwX_free( cdfX 175 if( cdfX != NULL ) cdfX = ptwX_free( cdfX );
176 if ( angular != NULL ) MCGIDI_angular_free 176 if ( angular != NULL ) MCGIDI_angular_free( smr, angular );
177 return( 1 ); 177 return( 1 );
178 } 178 }
179 /* 179 /*
180 ********************************************** 180 ************************************************************
181 */ 181 */
182 int MCGIDI_angular_sampleMu( statusMessageRepo 182 int MCGIDI_angular_sampleMu( statusMessageReporting *smr, MCGIDI_angular *angular, MCGIDI_quantitiesLookupModes &modes,
183 MCGIDI_decaySamplingInfo *decaySamplin 183 MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
184 184
185 double randomMu = decaySamplingInfo->rng( 185 double randomMu = decaySamplingInfo->rng( decaySamplingInfo->rngState );
186 MCGIDI_pdfsOfXGivenW_sampled sampled; 186 MCGIDI_pdfsOfXGivenW_sampled sampled;
187 187
188 switch( angular->type ) { 188 switch( angular->type ) {
189 case MCGIDI_angularType_isotropic : 189 case MCGIDI_angularType_isotropic :
190 decaySamplingInfo->frame = angular->fr 190 decaySamplingInfo->frame = angular->frame;
191 decaySamplingInfo->mu = 1. - 2. * deca 191 decaySamplingInfo->mu = 1. - 2. * decaySamplingInfo->rng( decaySamplingInfo->rngState );
192 break; 192 break;
193 case MCGIDI_angularType_linear : 193 case MCGIDI_angularType_linear :
194 decaySamplingInfo->frame = angular->fr 194 decaySamplingInfo->frame = angular->frame;
195 sampled.smr = smr; 195 sampled.smr = smr;
196 sampled.w = modes.getProjectileEnergy( 196 sampled.w = modes.getProjectileEnergy( );
197 MCGIDI_sampling_sampleX_from_pdfsOfXGi 197 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(angular->dists), &sampled, randomMu );
198 decaySamplingInfo->mu = sampled.x; 198 decaySamplingInfo->mu = sampled.x;
199 break; 199 break;
200 case MCGIDI_angularType_recoil : 200 case MCGIDI_angularType_recoil :
201 default : 201 default :
202 smr_setReportError2( smr, smr_unknownI 202 smr_setReportError2( smr, smr_unknownID, 1, "angular type = %d not supported", angular->type );
203 } 203 }
204 return( !smr_isOk( smr ) ); 204 return( !smr_isOk( smr ) );
205 } 205 }
206 206
207 #if defined __cplusplus 207 #if defined __cplusplus
208 } 208 }
209 #endif 209 #endif
210 210
211 211