| Geant4 Cross Reference |
1 /* 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( statusMess
21
22 MCGIDI_angular *angular;
23
24 if( ( angular = (MCGIDI_angular *) smr_mal
25 if( MCGIDI_angular_initialize( smr, angula
26 return( angular );
27 }
28 /*
29 **********************************************
30 */
31 int MCGIDI_angular_initialize( statusMessageRe
32
33 memset( angular, 0, sizeof( MCGIDI_angular
34 return( 0 );
35 }
36 /*
37 **********************************************
38 */
39 MCGIDI_angular *MCGIDI_angular_free( statusMes
40
41 MCGIDI_angular_release( smr, angular );
42 smr_freeMemory( (void **) &angular );
43 return( NULL );
44 }
45 /*
46 **********************************************
47 */
48 int MCGIDI_angular_release( statusMessageRepor
49
50
51 MCGIDI_sampling_pdfsOfXGivenW_release( smr
52
53 MCGIDI_angular_initialize( smr, angular );
54 return( 0 );
55 }
56 /*
57 **********************************************
58 */
59 int MCGIDI_angular_setTwoBodyMasses( statusMes
60 double productMass_MeV, double residualMas
61
62 if( angular == NULL ) return( 0 );
63 angular->projectileMass_MeV = projectileMa
64 angular->targetMass_MeV = targetMass_MeV;
65 angular->productMass_MeV = productMass_MeV
66 angular->residualMass_MeV = residualMass_M
67 return( 0 );
68 }
69 /*
70 **********************************************
71 */
72 int MCGIDI_angular_parseFromTOM( statusMessage
73
74 MCGIDI_angular *angular = NULL;
75 xDataTOM_element *angularElement, *linearE
76 char const *nativeData;
77 ptwXYPoints *pdfXY = NULL;
78 ptwXPoints *cdfX = NULL;
79 ptwXY_interpolation interpolationXY, inter
80
81 if( ( angularElement = xDataTOME_getOneEle
82 if( ( angular = MCGIDI_angular_new( smr )
83
84 if( ( nativeData = xDataTOM_getAttributesV
85 if( strcmp( nativeData, "isotropic" ) == 0
86 if( ( frameElement = xDataTOME_getOneE
87 smr_setReportError2( smr, smr_unkn
88 goto err;
89 }
90 angular->type = MCGIDI_angularType_iso
91 else if( strcmp( nativeData, "recoil" ) ==
92 angular->type = MCGIDI_angularType_rec
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 = &(angula
101 MCGIDI_pdfOfX *dist;
102 char const *energyUnit, *multiplicityP
103
104 if( ( linearElement = xDataTOME_getOne
105 if( ( linearElement = xDataTOME_ge
106 smr_setReportError2( smr, smr_
107 goto err;
108 }
109 }
110 frameElement = linearElement;
111
112 if( MCGIDI_fromTOM_interpolation( smr,
113 if( MCGIDI_fromTOM_interpolation( smr,
114 dists->interpolationWY = interpolation
115 dists->interpolationXY = interpolation
116
117 if( ( W_XYs = (xDataTOM_W_XYs *) xData
118 if( ( dists->Ws = (double *) smr_mallo
119 if( ( dists->dist = (MCGIDI_pdfOfX *)
120
121 energyUnit = xDataTOM_subAxes_getUnit(
122 if( !smr_isOk( smr ) ) goto err;
123 energyFactor = MCGIDI_misc_getUnitConv
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 * energy
130 if( ( pdfXY = MCGIDI_misc_dataFrom
131 if( ptwXY_simpleCoalescePoints( pd
132 dist->numberOfXs = n = (int) ptwXY
133
134 if( ( dist->Xs = (double *) smr_ma
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_
141 dist->Xs[j] = point->x;
142 dist->pdf[j] = point->y;
143 }
144
145 if( ( cdfX = ptwXY_runningIntegral
146 smr_setReportError2( smr, smr_
147 goto err;
148 }
149 norm = ptwX_getPointAtIndex_Unsafe
150 if( norms != NULL ) {
151 ptwXY_setValueAtX( norms, XYs-
152 else if( std::fabs( 1. - norm ) >
153 smr_setReportError2( smr, smr_
154 goto err;
155 }
156 for( j = 0; j < n; j++ ) dist->cdf
157 for( j = 0; j < n; j++ ) dist->pdf
158 pdfXY = ptwXY_free( pdfXY );
159 cdfX = ptwX_free( cdfX );
160 }
161 angular->type = MCGIDI_angularType_lin
162 }
163
164 if( frameElement != NULL ) {
165 if( ( angular->frame = MCGIDI_misc_get
166 }
167
168 distribution->angular = angular;
169 distribution->type = MCGIDI_distributionTy
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
177 return( 1 );
178 }
179 /*
180 **********************************************
181 */
182 int MCGIDI_angular_sampleMu( statusMessageRepo
183 MCGIDI_decaySamplingInfo *decaySamplin
184
185 double randomMu = decaySamplingInfo->rng(
186 MCGIDI_pdfsOfXGivenW_sampled sampled;
187
188 switch( angular->type ) {
189 case MCGIDI_angularType_isotropic :
190 decaySamplingInfo->frame = angular->fr
191 decaySamplingInfo->mu = 1. - 2. * deca
192 break;
193 case MCGIDI_angularType_linear :
194 decaySamplingInfo->frame = angular->fr
195 sampled.smr = smr;
196 sampled.w = modes.getProjectileEnergy(
197 MCGIDI_sampling_sampleX_from_pdfsOfXGi
198 decaySamplingInfo->mu = sampled.x;
199 break;
200 case MCGIDI_angularType_recoil :
201 default :
202 smr_setReportError2( smr, smr_unknownI
203 }
204 return( !smr_isOk( smr ) );
205 }
206
207 #if defined __cplusplus
208 }
209 #endif
210
211