| 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 6
7 #include "MCGIDI_fromTOM.h" 7 #include "MCGIDI_fromTOM.h"
8 #include "MCGIDI.h" 8 #include "MCGIDI.h"
9 #include "MCGIDI_misc.h" 9 #include "MCGIDI_misc.h"
10 10
11 #if defined __cplusplus 11 #if defined __cplusplus
12 namespace GIDI { 12 namespace GIDI {
13 using namespace GIDI; 13 using namespace GIDI;
14 #endif 14 #endif
15 15
16 static int MCGIDI_energyAngular_linear_parseFr 16 static int MCGIDI_energyAngular_linear_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution );
17 /* 17 /*
18 ********************************************** 18 ************************************************************
19 */ 19 */
20 MCGIDI_energyAngular *MCGIDI_energyAngular_new 20 MCGIDI_energyAngular *MCGIDI_energyAngular_new( statusMessageReporting *smr ) {
21 21
22 MCGIDI_energyAngular *energyAngular; 22 MCGIDI_energyAngular *energyAngular;
23 23
24 if( ( energyAngular = (MCGIDI_energyAngula 24 if( ( energyAngular = (MCGIDI_energyAngular *) smr_malloc2( smr, sizeof( MCGIDI_energyAngular ), 0, "energyAngular" ) ) == NULL ) return( NULL );
25 if( MCGIDI_energyAngular_initialize( smr, 25 if( MCGIDI_energyAngular_initialize( smr, energyAngular ) ) energyAngular = MCGIDI_energyAngular_free( smr, energyAngular );
26 return( energyAngular ); 26 return( energyAngular );
27 } 27 }
28 /* 28 /*
29 ********************************************** 29 ************************************************************
30 */ 30 */
31 int MCGIDI_energyAngular_initialize( statusMes 31 int MCGIDI_energyAngular_initialize( statusMessageReporting * /*smr*/, MCGIDI_energyAngular *energyAngular ) {
32 32
33 memset( energyAngular, 0, sizeof( MCGIDI_e 33 memset( energyAngular, 0, sizeof( MCGIDI_energyAngular ) );
34 return( 0 ); 34 return( 0 );
35 } 35 }
36 /* 36 /*
37 ********************************************** 37 ************************************************************
38 */ 38 */
39 MCGIDI_energyAngular *MCGIDI_energyAngular_fre 39 MCGIDI_energyAngular *MCGIDI_energyAngular_free( statusMessageReporting *smr, MCGIDI_energyAngular *energyAngular ) {
40 40
41 MCGIDI_energyAngular_release( smr, energyA 41 MCGIDI_energyAngular_release( smr, energyAngular );
42 smr_freeMemory( (void **) &energyAngular ) 42 smr_freeMemory( (void **) &energyAngular );
43 return( NULL ); 43 return( NULL );
44 } 44 }
45 /* 45 /*
46 ********************************************** 46 ************************************************************
47 */ 47 */
48 int MCGIDI_energyAngular_release( statusMessag 48 int MCGIDI_energyAngular_release( statusMessageReporting *smr, MCGIDI_energyAngular *energyAngular ) {
49 49
50 int i; 50 int i;
51 51
52 for( i = 0; i < energyAngular->pdfOfEpGive 52 for( i = 0; i < energyAngular->pdfOfEpGivenE.numberOfWs; i++ ) {
53 MCGIDI_sampling_pdfsOfXGivenW_release( 53 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(energyAngular->pdfOfMuGivenEAndEp[i]) );
54 } 54 }
55 smr_freeMemory( (void **) &(energyAngular- 55 smr_freeMemory( (void **) &(energyAngular->pdfOfMuGivenEAndEp) );
56 MCGIDI_sampling_pdfsOfXGivenW_release( smr 56 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(energyAngular->pdfOfEpGivenE) );
57 MCGIDI_energyAngular_initialize( smr, ener 57 MCGIDI_energyAngular_initialize( smr, energyAngular );
58 58
59 return( 0 ); 59 return( 0 );
60 } 60 }
61 /* 61 /*
62 ********************************************** 62 ************************************************************
63 */ 63 */
64 int MCGIDI_energyAngular_parseFromTOM( statusM 64 int MCGIDI_energyAngular_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution ) {
65 65
66 xDataTOM_element *energyAngularElement; 66 xDataTOM_element *energyAngularElement;
67 char const *nativeData; 67 char const *nativeData;
68 68
69 if( ( energyAngularElement = xDataTOME_get 69 if( ( energyAngularElement = xDataTOME_getOneElementByName( smr, element, "energyAngular", 1 ) ) == NULL ) goto err;
70 70
71 if( ( nativeData = xDataTOM_getAttributesV 71 if( ( nativeData = xDataTOM_getAttributesValueInElement( energyAngularElement, "nativeData" ) ) == NULL ) goto err;
72 if( strcmp( nativeData, "KalbachMann" ) == 72 if( strcmp( nativeData, "KalbachMann" ) == 0 ) {
73 return( MCGIDI_KalbachMann_parseFromTO 73 return( MCGIDI_KalbachMann_parseFromTOM( smr, energyAngularElement, distribution ) ); }
74 else if( strcmp( nativeData, "linear" ) == 74 else if( strcmp( nativeData, "linear" ) == 0 ) {
75 return( MCGIDI_energyAngular_linear_pa 75 return( MCGIDI_energyAngular_linear_parseFromTOM( smr, energyAngularElement, distribution ) ); }
76 else { 76 else {
77 smr_setReportError2( smr, smr_unknownI 77 smr_setReportError2( smr, smr_unknownID, 1, "energyAngular nativeData = '%s' not supported", nativeData );
78 goto err; 78 goto err;
79 } 79 }
80 80
81 return( 0 ); 81 return( 0 );
82 82
83 err: 83 err:
84 return( 1 ); 84 return( 1 );
85 } 85 }
86 /* 86 /*
87 ********************************************** 87 ************************************************************
88 */ 88 */
89 static int MCGIDI_energyAngular_linear_parseFr 89 static int MCGIDI_energyAngular_linear_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution ) {
90 90
91 int iV, iW; 91 int iV, iW;
92 double y, norm, energyInFactor, energyOutF 92 double y, norm, energyInFactor, energyOutFactor;
93 char const *energyUnit, *multiplicityProba 93 char const *energyUnit, *multiplicityProbabilityUnits[2] = { "", "1/MeV" };
94 xDataTOM_element *linear; 94 xDataTOM_element *linear;
95 ptwXY_interpolation interpolationXY, inter 95 ptwXY_interpolation interpolationXY, interpolationWY, interpolationVY;
96 xDataTOM_XYs *XYs; 96 xDataTOM_XYs *XYs;
97 xDataTOM_W_XYs *W_XYs; 97 xDataTOM_W_XYs *W_XYs;
98 xDataTOM_V_W_XYs *V_W_XYs; 98 xDataTOM_V_W_XYs *V_W_XYs;
99 MCGIDI_pdfsOfXGivenW *pdfOfEpGivenE, *pdfO 99 MCGIDI_pdfsOfXGivenW *pdfOfEpGivenE, *pdfOfMuGivenEAndEp = NULL, *pdfOfMuGivenEAndEp2 = NULL;
100 MCGIDI_energyAngular *energyAngular = NULL 100 MCGIDI_energyAngular *energyAngular = NULL;
101 ptwXYPoints *pdfXY1 = NULL, *pdfXY2 = NULL 101 ptwXYPoints *pdfXY1 = NULL, *pdfXY2 = NULL;
102 nfu_status status; 102 nfu_status status;
103 103
104 if( ( linear = xDataTOME_getOneElementByNa 104 if( ( linear = xDataTOME_getOneElementByName( smr, element, "linear", 1 ) ) == NULL ) goto err;
105 105
106 if( MCGIDI_fromTOM_interpolation( smr, lin 106 if( MCGIDI_fromTOM_interpolation( smr, linear, 0, &interpolationVY ) ) goto err;
107 if( MCGIDI_fromTOM_interpolation( smr, lin 107 if( MCGIDI_fromTOM_interpolation( smr, linear, 1, &interpolationWY ) ) goto err;
108 if( MCGIDI_fromTOM_interpolation( smr, lin 108 if( MCGIDI_fromTOM_interpolation( smr, linear, 2, &interpolationXY ) ) goto err;
109 109
110 if( ( energyAngular = MCGIDI_energyAngular 110 if( ( energyAngular = MCGIDI_energyAngular_new( smr ) ) == NULL ) goto err;
111 if( ( energyAngular->frame = MCGIDI_misc_g 111 if( ( energyAngular->frame = MCGIDI_misc_getProductFrame( smr, linear ) ) == xDataTOM_frame_invalid ) goto err;
112 112
113 pdfOfEpGivenE = &(energyAngular->pdfOfEpGi 113 pdfOfEpGivenE = &(energyAngular->pdfOfEpGivenE);
114 pdfOfEpGivenE->interpolationWY = interpola 114 pdfOfEpGivenE->interpolationWY = interpolationVY;
115 pdfOfEpGivenE->interpolationXY = interpola 115 pdfOfEpGivenE->interpolationXY = interpolationWY;
116 116
117 if( ( V_W_XYs = (xDataTOM_V_W_XYs *) xData 117 if( ( V_W_XYs = (xDataTOM_V_W_XYs *) xDataTOME_getXDataIfID( smr, linear, "V_W_XYs" ) ) == NULL ) goto err;
118 if( ( pdfOfEpGivenE->Ws = (double *) smr_m 118 if( ( pdfOfEpGivenE->Ws = (double *) smr_malloc2( smr, V_W_XYs->length * sizeof( double ), 1, "pdfOfEpGivenE->Ws" ) ) == NULL ) goto err;
119 if( ( pdfOfEpGivenE->dist = (MCGIDI_pdfOfX 119 if( ( pdfOfEpGivenE->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, V_W_XYs->length * sizeof( MCGIDI_pdfOfX ), 0, "pdfOfEpGivenE->dist" ) ) == NULL ) goto err;
120 if( ( pdfOfMuGivenEAndEp = (MCGIDI_pdfsOfX 120 if( ( pdfOfMuGivenEAndEp = (MCGIDI_pdfsOfXGivenW *) smr_malloc2( smr, V_W_XYs->length * sizeof( MCGIDI_pdfsOfXGivenW ), 1, "pdfOfMuGivenEAndEp" ) ) == NULL ) goto err;
121 121
122 energyUnit = xDataTOM_subAxes_getUnit( smr 122 energyUnit = xDataTOM_subAxes_getUnit( smr, &(V_W_XYs->subAxes), 0 );
123 if( !smr_isOk( smr ) ) goto err; 123 if( !smr_isOk( smr ) ) goto err;
124 energyInFactor = MCGIDI_misc_getUnitConver 124 energyInFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyUnit, "MeV" );
125 if( !smr_isOk( smr ) ) goto err; 125 if( !smr_isOk( smr ) ) goto err;
126 126
127 energyUnit = xDataTOM_subAxes_getUnit( smr 127 energyUnit = xDataTOM_subAxes_getUnit( smr, &(V_W_XYs->subAxes), 1 );
128 if( !smr_isOk( smr ) ) goto err; 128 if( !smr_isOk( smr ) ) goto err;
129 energyOutFactor = MCGIDI_misc_getUnitConve 129 energyOutFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyUnit, "MeV" );
130 if( !smr_isOk( smr ) ) goto err; 130 if( !smr_isOk( smr ) ) goto err;
131 131
132 for( iV = 0; iV < V_W_XYs->length; iV++ ) 132 for( iV = 0; iV < V_W_XYs->length; iV++ ) {
133 W_XYs = &(V_W_XYs->W_XYs[iV]); 133 W_XYs = &(V_W_XYs->W_XYs[iV]);
134 pdfOfMuGivenEAndEp2 = &(pdfOfMuGivenEA 134 pdfOfMuGivenEAndEp2 = &(pdfOfMuGivenEAndEp[iV]);
135 pdfOfMuGivenEAndEp2->interpolationWY = 135 pdfOfMuGivenEAndEp2->interpolationWY = interpolationWY;
136 pdfOfMuGivenEAndEp2->interpolationXY = 136 pdfOfMuGivenEAndEp2->interpolationXY = interpolationXY;
137 if( ( pdfXY2 = ptwXY_new( interpolatio 137 if( ( pdfXY2 = ptwXY_new( interpolationWY, NULL, 2., 1e-6, W_XYs->length, 10, &status, 0 ) ) == NULL ) goto errA;
138 if( ( pdfOfMuGivenEAndEp2->Ws = (doubl 138 if( ( pdfOfMuGivenEAndEp2->Ws = (double *) smr_malloc2( smr, W_XYs->length * sizeof( double ), 1, "pdfOfMuGivenEAndEp2->Ws" ) ) == NULL ) goto err;
139 if( ( pdfOfMuGivenEAndEp2->dist = (MCG 139 if( ( pdfOfMuGivenEAndEp2->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, W_XYs->length * sizeof( MCGIDI_pdfOfX ), 0, "pdfOfMuGivenEAndEp2->dist" ) ) == NULL ) goto err;
140 for( iW = 0; iW < W_XYs->length; iW++ 140 for( iW = 0; iW < W_XYs->length; iW++ ) {
141 XYs = &(W_XYs->XYs[iW]); 141 XYs = &(W_XYs->XYs[iW]);
142 if( ( pdfXY1 = MCGIDI_misc_dataFr 142 if( ( pdfXY1 = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, interpolationXY, multiplicityProbabilityUnits ) ) == NULL ) goto err;
143 y = ptwXY_integrateDomain( pdfXY1, 143 y = ptwXY_integrateDomain( pdfXY1, &status );
144 if( ( status = ptwXY_setValueAtX( 144 if( ( status = ptwXY_setValueAtX( pdfXY2, energyOutFactor * XYs->value, y ) ) != nfu_Okay ) goto errA;
>> 145 if( status != nfu_Okay ) goto errA;
145 146
146 if( y == 0 ) { 147 if( y == 0 ) {
147 if( ( status = ptwXY_add_doubl 148 if( ( status = ptwXY_add_double( pdfXY1, 0.5 ) ) != nfu_Okay ) goto errA;
148 } 149 }
149 pdfOfMuGivenEAndEp2->Ws[iW] = ener 150 pdfOfMuGivenEAndEp2->Ws[iW] = energyOutFactor * XYs->value;
150 if( MCGIDI_fromTOM_pdfOfX( smr, pd 151 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY1, &(pdfOfMuGivenEAndEp2->dist[iW]), &norm ) ) goto err;
151 pdfOfMuGivenEAndEp2->numberOfWs++; 152 pdfOfMuGivenEAndEp2->numberOfWs++;
152 153
153 pdfXY1 = ptwXY_free( pdfXY1 ); 154 pdfXY1 = ptwXY_free( pdfXY1 );
154 } 155 }
155 pdfOfEpGivenE->Ws[iV] = energyInFactor 156 pdfOfEpGivenE->Ws[iV] = energyInFactor * W_XYs->value;
156 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY2 157 if( MCGIDI_fromTOM_pdfOfX( smr, pdfXY2, &(pdfOfEpGivenE->dist[iV]), &norm ) ) goto err;
157 pdfOfEpGivenE->numberOfWs++; 158 pdfOfEpGivenE->numberOfWs++;
158 159
159 pdfXY2 = ptwXY_free( pdfXY2 ); 160 pdfXY2 = ptwXY_free( pdfXY2 );
160 } 161 }
161 energyAngular->pdfOfMuGivenEAndEp = pdfOfM 162 energyAngular->pdfOfMuGivenEAndEp = pdfOfMuGivenEAndEp;
162 distribution->energyAngular = energyAngula 163 distribution->energyAngular = energyAngular;
163 distribution->type = MCGIDI_distributionTy 164 distribution->type = MCGIDI_distributionType_energyAngular_e;
164 165
165 return( 0 ); 166 return( 0 );
166 167
167 errA: 168 errA:
168 smr_setReportError2( smr, smr_unknownID, 1 169 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_integrateDomain err = %d: %s\n", status, nfu_statusMessage( status ) );
169 err: 170 err:
170 if( pdfXY1 != NULL ) ptwXY_free( pdfXY1 ); 171 if( pdfXY1 != NULL ) ptwXY_free( pdfXY1 );
171 if( pdfXY2 != NULL ) ptwXY_free( pdfXY2 ); 172 if( pdfXY2 != NULL ) ptwXY_free( pdfXY2 );
172 if( energyAngular != NULL ) MCGIDI_energyA 173 if( energyAngular != NULL ) MCGIDI_energyAngular_free( smr, energyAngular );
173 /* ????????? Need to free pdfOfMuGivenEAndEp, 174 /* ????????? Need to free pdfOfMuGivenEAndEp, now may be handled by MCGIDI_energyAngular_free. Need to check. */
174 return( 1 ); 175 return( 1 );
175 } 176 }
176 /* 177 /*
177 ********************************************** 178 ************************************************************
178 */ 179 */
179 int MCGIDI_energyAngular_sampleDistribution( s 180 int MCGIDI_energyAngular_sampleDistribution( statusMessageReporting *smr, MCGIDI_distribution *distribution, MCGIDI_quantitiesLookupModes &modes,
180 MCGIDI_decaySamplingInfo *decaySamplin 181 MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
181 182
182 double Ep; 183 double Ep;
183 MCGIDI_energyAngular *energyAngular = dist 184 MCGIDI_energyAngular *energyAngular = distribution->energyAngular;
184 185
185 MCGIDI_sampling_doubleDistribution( smr, & 186 MCGIDI_sampling_doubleDistribution( smr, &(energyAngular->pdfOfEpGivenE), energyAngular->pdfOfMuGivenEAndEp, modes, decaySamplingInfo );
186 Ep = decaySamplingInfo->mu; 187 Ep = decaySamplingInfo->mu;
187 decaySamplingInfo->mu = decaySamplingInfo- 188 decaySamplingInfo->mu = decaySamplingInfo->Ep;
188 decaySamplingInfo->Ep = Ep; 189 decaySamplingInfo->Ep = Ep;
189 decaySamplingInfo->frame = energyAngular-> 190 decaySamplingInfo->frame = energyAngular->frame;
190 191
191 return( 0 ); 192 return( 0 );
192 } 193 }
193 194
194 #if defined __cplusplus 195 #if defined __cplusplus
195 } 196 }
196 #endif 197 #endif
197 198
198 199