| 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 #define _USE_MATH_DEFINES 6 #define _USE_MATH_DEFINES
7 #include <cmath> 7 #include <cmath>
8 8
9 9
10 #include "MCGIDI_fromTOM.h" 10 #include "MCGIDI_fromTOM.h"
11 #include "MCGIDI_misc.h" 11 #include "MCGIDI_misc.h"
12 #include "MCGIDI_private.h" 12 #include "MCGIDI_private.h"
13 #include <nf_specialFunctions.h> 13 #include <nf_specialFunctions.h>
14 14
15 #if defined __cplusplus 15 #if defined __cplusplus
16 #include "G4Exp.hh" 16 #include "G4Exp.hh"
17 #include "G4Log.hh" 17 #include "G4Log.hh"
18 #include "G4Pow.hh" 18 #include "G4Pow.hh"
19 namespace GIDI { 19 namespace GIDI {
20 using namespace GIDI; 20 using namespace GIDI;
21 #endif 21 #endif
22 22
23 static int MCGIDI_energy_parseWeightFromTOM( s 23 static int MCGIDI_energy_parseWeightFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energyWeightedFunctional *weightedFunctional );
24 static int MCGIDI_energy_parseWeightedFunction 24 static int MCGIDI_energy_parseWeightedFunctionalsFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
25 static int MCGIDI_energy_parseGeneralEvaporati 25 static int MCGIDI_energy_parseGeneralEvaporationFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
26 static int MCGIDI_energy_parseEvaporationFromT 26 static int MCGIDI_energy_parseEvaporationFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
27 static int MCGIDI_energy_parseWattFromTOM( sta 27 static int MCGIDI_energy_parseWattFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
28 static int MCGIDI_energy_parseSimpleMaxwellian 28 static int MCGIDI_energy_parseSimpleMaxwellianFissionFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy );
29 static int MCGIDI_energy_parseMadlandNixFromTO 29 static int MCGIDI_energy_parseMadlandNixFromTOM( statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy );
30 static nfu_status MCGIDI_energy_parseMadlandNi 30 static nfu_status MCGIDI_energy_parseMadlandNixFromTOM_callback( double x, double *y, void *argList );
31 static double MCGIDI_energy_parseMadlandNixFro 31 static double MCGIDI_energy_parseMadlandNixFromTOM_callback_g( double Ep, double EFL, double T_M, nfu_status *status );
32 static int MCGIDI_energy_parseNBodyPhaseSpaceF 32 static int MCGIDI_energy_parseNBodyPhaseSpaceFromTOM( statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy,
33 MCGIDI_distribution *distribution ); 33 MCGIDI_distribution *distribution );
34 34
35 static int MCGIDI_energy_sampleSimpleMaxwellia 35 static int MCGIDI_energy_sampleSimpleMaxwellianFission( statusMessageReporting *smr, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo );
36 static int MCGIDI_energy_sampleEvaporation( st 36 static int MCGIDI_energy_sampleEvaporation( statusMessageReporting *smr, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo );
37 static int MCGIDI_energy_sampleWatt( statusMes 37 static int MCGIDI_energy_sampleWatt( statusMessageReporting *smr, double e_in_U, double Watt_a, double Watt_b, MCGIDI_decaySamplingInfo *decaySamplingInfo );
38 static int MCGIDI_energy_sampleWeightedFunctio 38 static int MCGIDI_energy_sampleWeightedFunctional( statusMessageReporting *smr, MCGIDI_energy *energy,
39 MCGIDI_quantitiesLookupModes &modes, MCGID 39 MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo );
40 static nfu_status MCGIDI_energy_NBodyPhaseSpac 40 static nfu_status MCGIDI_energy_NBodyPhaseSpacePDF_callback( double x, double *y, void *argList );
41 /* 41 /*
42 ********************************************** 42 ************************************************************
43 */ 43 */
44 MCGIDI_energy *MCGIDI_energy_new( statusMessag 44 MCGIDI_energy *MCGIDI_energy_new( statusMessageReporting *smr ) {
45 45
46 MCGIDI_energy *energy; 46 MCGIDI_energy *energy;
47 47
48 if( ( energy = (MCGIDI_energy *) smr_mallo 48 if( ( energy = (MCGIDI_energy *) smr_malloc2( smr, sizeof( MCGIDI_energy ), 0, "energy" ) ) == NULL ) return( NULL );
49 if( MCGIDI_energy_initialize( smr, energy 49 if( MCGIDI_energy_initialize( smr, energy ) ) energy = MCGIDI_energy_free( smr, energy );
50 return( energy ); 50 return( energy );
51 } 51 }
52 /* 52 /*
53 ********************************************** 53 ************************************************************
54 */ 54 */
55 int MCGIDI_energy_initialize( statusMessageRep 55 int MCGIDI_energy_initialize( statusMessageReporting * /*smr*/, MCGIDI_energy *energy ) {
56 56
57 memset( energy, 0, sizeof( MCGIDI_energy ) 57 memset( energy, 0, sizeof( MCGIDI_energy ) );
58 return( 0 ); 58 return( 0 );
59 } 59 }
60 /* 60 /*
61 ********************************************** 61 ************************************************************
62 */ 62 */
63 MCGIDI_energy *MCGIDI_energy_free( statusMessa 63 MCGIDI_energy *MCGIDI_energy_free( statusMessageReporting *smr, MCGIDI_energy *energy ) {
64 64
65 MCGIDI_energy_release( smr, energy ); 65 MCGIDI_energy_release( smr, energy );
66 smr_freeMemory( (void **) &energy ); 66 smr_freeMemory( (void **) &energy );
67 return( NULL ); 67 return( NULL );
68 } 68 }
69 /* 69 /*
70 ********************************************** 70 ************************************************************
71 */ 71 */
72 int MCGIDI_energy_release( statusMessageReport 72 int MCGIDI_energy_release( statusMessageReporting *smr, MCGIDI_energy *energy ) {
73 73
74 int i; 74 int i;
75 75
76 MCGIDI_sampling_pdfsOfXGivenW_release( smr 76 MCGIDI_sampling_pdfsOfXGivenW_release( smr, &(energy->dists) );
77 if( energy->theta ) energy->theta = ptwXY_ 77 if( energy->theta ) energy->theta = ptwXY_free( energy->theta );
78 if( energy->Watt_a ) energy->Watt_a = ptwX 78 if( energy->Watt_a ) energy->Watt_a = ptwXY_free( energy->Watt_a );
79 if( energy->Watt_b ) energy->Watt_b = ptwX 79 if( energy->Watt_b ) energy->Watt_b = ptwXY_free( energy->Watt_b );
80 if( ( energy->type == MCGIDI_energyType_ge 80 if( ( energy->type == MCGIDI_energyType_generalEvaporation ) || ( energy->type == MCGIDI_energyType_NBodyPhaseSpace ) ) {
81 MCGIDI_sampling_pdfsOfX_release( smr, 81 MCGIDI_sampling_pdfsOfX_release( smr, &(energy->g) ); }
82 else if( energy->type == MCGIDI_energyType 82 else if( energy->type == MCGIDI_energyType_weightedFunctional ) {
83 for( i = 0; i < energy->weightedFuncti 83 for( i = 0; i < energy->weightedFunctionals.numberOfWeights; i++ ) {
84 ptwXY_free( energy->weightedFuncti 84 ptwXY_free( energy->weightedFunctionals.weightedFunctional[i].weight );
85 MCGIDI_energy_free( smr, energy->w 85 MCGIDI_energy_free( smr, energy->weightedFunctionals.weightedFunctional[i].energy );
86 } 86 }
87 } 87 }
88 88
89 MCGIDI_energy_initialize( smr, energy ); 89 MCGIDI_energy_initialize( smr, energy );
90 return( 0 ); 90 return( 0 );
91 } 91 }
92 /* 92 /*
93 ********************************************** 93 ************************************************************
94 */ 94 */
95 int MCGIDI_energy_parseFromTOM( statusMessageR 95 int MCGIDI_energy_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution, ptwXYPoints *norms,
96 enum MCGIDI_energyType energyType, dou 96 enum MCGIDI_energyType energyType, double gammaEnergy_MeV ) {
97 97
98 MCGIDI_energy *energy = NULL; 98 MCGIDI_energy *energy = NULL;
99 xDataTOM_element *energyElement, *linearEl 99 xDataTOM_element *energyElement, *linearElement, *functional, *frameElement;
100 char const *nativeData; 100 char const *nativeData;
101 double projectileMass_MeV, targetMass_MeV; 101 double projectileMass_MeV, targetMass_MeV;
102 102
103 if( ( energy = MCGIDI_energy_new( smr ) ) 103 if( ( energy = MCGIDI_energy_new( smr ) ) == NULL ) goto err;
104 104
105 projectileMass_MeV = MCGIDI_product_getPro 105 projectileMass_MeV = MCGIDI_product_getProjectileMass_MeV( smr, distribution->product );
106 targetMass_MeV = MCGIDI_product_getTargetM 106 targetMass_MeV = MCGIDI_product_getTargetMass_MeV( smr, distribution->product );
107 energy->e_inCOMFactor = targetMass_MeV / ( 107 energy->e_inCOMFactor = targetMass_MeV / ( projectileMass_MeV + targetMass_MeV );
108 108
109 if( ( energyType == MCGIDI_energyType_prim 109 if( ( energyType == MCGIDI_energyType_primaryGamma ) || ( energyType == MCGIDI_energyType_discreteGamma ) ) {
110 energy->type = energyType; 110 energy->type = energyType;
111 energy->gammaEnergy_MeV = gammaEnergy_ 111 energy->gammaEnergy_MeV = gammaEnergy_MeV;
112 energy->frame = xDataTOM_frame_lab; 112 energy->frame = xDataTOM_frame_lab; /* BRB. This should not be hardwired?????? Probably needs to be changed in GND also. */
113 if( energyType == MCGIDI_energyType_pr 113 if( energyType == MCGIDI_energyType_primaryGamma ) energy->primaryGammaMassFactor = energy->e_inCOMFactor; }
114 else { 114 else {
115 if( ( energyElement = xDataTOME_getOne 115 if( ( energyElement = xDataTOME_getOneElementByName( smr, element, "energy", 1 ) ) == NULL ) goto err;
116 if( ( nativeData = xDataTOM_getAttribu 116 if( ( nativeData = xDataTOM_getAttributesValueInElement( energyElement, "nativeData" ) ) == NULL ) goto err;
117 if( ( linearElement = xDataTOME_getOne 117 if( ( linearElement = xDataTOME_getOneElementByName( NULL, energyElement, "linear", 0 ) ) == NULL )
118 linearElement = xDataTOME_getOneEl 118 linearElement = xDataTOME_getOneElementByName( NULL, energyElement, "pointwise", 0 );
119 if( linearElement == NULL ) { 119 if( linearElement == NULL ) {
120 if( ( functional = xDataTOME_getOn 120 if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "generalEvaporation", 0 ) ) != NULL ) {
121 if( MCGIDI_energy_parseGeneral 121 if( MCGIDI_energy_parseGeneralEvaporationFromTOM( smr, functional, energy ) ) goto err; }
122 else if( ( functional = xDataTOME_ 122 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "simpleMaxwellianFission", 0 ) ) != NULL ) {
123 if( MCGIDI_energy_parseSimpleM 123 if( MCGIDI_energy_parseSimpleMaxwellianFissionFromTOM( smr, functional, energy ) ) goto err; }
124 else if( ( functional = xDataTOME_ 124 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "evaporation", 0 ) ) != NULL ) {
125 if( MCGIDI_energy_parseEvapora 125 if( MCGIDI_energy_parseEvaporationFromTOM( smr, functional, energy ) ) goto err; }
126 else if( ( functional = xDataTOME_ 126 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "Watt", 0 ) ) != NULL ) {
127 if( MCGIDI_energy_parseWattFro 127 if( MCGIDI_energy_parseWattFromTOM( smr, functional, energy ) ) goto err; }
128 else if( ( functional = xDataTOME_ 128 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "MadlandNix", 0 ) ) != NULL ) {
129 if( MCGIDI_energy_parseMadland 129 if( MCGIDI_energy_parseMadlandNixFromTOM( smr, functional, energy ) ) goto err; }
130 else if( ( functional = xDataTOME_ 130 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "NBodyPhaseSpace", 0 ) ) != NULL ) {
131 if( MCGIDI_energy_parseNBodyPh 131 if( MCGIDI_energy_parseNBodyPhaseSpaceFromTOM( smr, functional, energy, distribution ) ) goto err; }
132 else if( ( functional = xDataTOME_ 132 else if( ( functional = xDataTOME_getOneElementByName( NULL, energyElement, "weightedFunctionals", 0 ) ) != NULL ) {
133 if( MCGIDI_energy_parseWeighte 133 if( MCGIDI_energy_parseWeightedFunctionalsFromTOM( smr, functional, energy ) ) goto err; }
134 else { 134 else {
135 smr_setReportError2( smr, smr_ 135 smr_setReportError2( smr, smr_unknownID, 1, "unsupported energy type: nativeData = '%s'", nativeData );
136 goto err; 136 goto err;
137 } 137 }
138 frameElement = functional; } 138 frameElement = functional; }
139 else { 139 else {
140 char const *toUnits[3] = { "MeV", 140 char const *toUnits[3] = { "MeV", "MeV", "1/MeV" };
141 141
142 frameElement = linearElement; 142 frameElement = linearElement;
143 if( MCGIDI_fromTOM_pdfsOfXGivenW( 143 if( MCGIDI_fromTOM_pdfsOfXGivenW( smr, linearElement, &(energy->dists), norms, toUnits ) ) goto err;
144 energy->type = MCGIDI_energyType_l 144 energy->type = MCGIDI_energyType_linear;
145 } 145 }
146 if( ( energy->frame = MCGIDI_misc_getP 146 if( ( energy->frame = MCGIDI_misc_getProductFrame( smr, frameElement ) ) == xDataTOM_frame_invalid ) goto err;
147 } 147 }
148 distribution->energy = energy; 148 distribution->energy = energy;
149 149
150 return( 0 ); 150 return( 0 );
151 151
152 err: 152 err:
153 if( energy != NULL ) MCGIDI_energy_free( s 153 if( energy != NULL ) MCGIDI_energy_free( smr, energy );
154 return( 1 ); 154 return( 1 );
155 } 155 }
156 /* 156 /*
157 ********************************************** 157 ************************************************************
158 */ 158 */
159 static int MCGIDI_energy_parseWeightedFunction 159 static int MCGIDI_energy_parseWeightedFunctionalsFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
160 160
161 int i; 161 int i;
162 xDataTOM_element *child; 162 xDataTOM_element *child;
163 163
164 for( i = 0, child = xDataTOME_getFirstElem 164 for( i = 0, child = xDataTOME_getFirstElement( element ); child != NULL; i++, child = xDataTOME_getNextElement( child ) ) {
165 if( strcmp( child->name, "weighted" ) 165 if( strcmp( child->name, "weighted" ) ) goto err;
166 if( MCGIDI_energy_parseWeightFromTOM( 166 if( MCGIDI_energy_parseWeightFromTOM( smr, child, &(energy->weightedFunctionals.weightedFunctional[i]) ) ) goto err;
167 energy->weightedFunctionals.numberOfWe 167 energy->weightedFunctionals.numberOfWeights++;
168 } 168 }
169 energy->type = MCGIDI_energyType_weightedF 169 energy->type = MCGIDI_energyType_weightedFunctional;
170 return( 0 ); 170 return( 0 );
171 171
172 err: 172 err:
173 return( 1 ); 173 return( 1 );
174 } 174 }
175 /* 175 /*
176 ********************************************** 176 ************************************************************
177 */ 177 */
178 static int MCGIDI_energy_parseWeightFromTOM( s 178 static int MCGIDI_energy_parseWeightFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energyWeightedFunctional *weightedFunctional ) {
179 179
180 xDataTOM_element *child; 180 xDataTOM_element *child;
181 MCGIDI_energy *energy = NULL; 181 MCGIDI_energy *energy = NULL;
182 ptwXYPoints *weight = NULL; 182 ptwXYPoints *weight = NULL;
183 char const *toUnits[2] = { "MeV", "" }; 183 char const *toUnits[2] = { "MeV", "" };
184 184
185 if( ( energy = MCGIDI_energy_new( smr ) ) 185 if( ( energy = MCGIDI_energy_new( smr ) ) == NULL ) goto err;
186 for( child = xDataTOME_getFirstElement( el 186 for( child = xDataTOME_getFirstElement( element ); child != NULL; child = xDataTOME_getNextElement( child ) ) {
187 if( strcmp( child->name, "weight" ) == 187 if( strcmp( child->name, "weight" ) == 0 ) {
188 if( ( weight = MCGIDI_misc_dataFro 188 if( ( weight = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, child, toUnits ) ) == NULL ) goto err; }
189 else if( strcmp( child->name, "evapora 189 else if( strcmp( child->name, "evaporation" ) == 0 ) {
190 if( MCGIDI_energy_parseEvaporation 190 if( MCGIDI_energy_parseEvaporationFromTOM( smr, child, energy ) ) goto err; }
191 else { 191 else {
192 smr_setReportError2( smr, smr_unkn 192 smr_setReportError2( smr, smr_unknownID, 1, "unsupported energy type = '%s' in weighted functional", child->name );
193 goto err; 193 goto err;
194 } 194 }
195 } 195 }
196 weightedFunctional->weight = weight; 196 weightedFunctional->weight = weight;
197 weightedFunctional->energy = energy; 197 weightedFunctional->energy = energy;
198 return( 0 ); 198 return( 0 );
199 199
200 err: 200 err:
201 if( weight != NULL ) ptwXY_free( weight ); 201 if( weight != NULL ) ptwXY_free( weight );
202 if( energy != NULL ) MCGIDI_energy_free( s 202 if( energy != NULL ) MCGIDI_energy_free( smr, energy );
203 return( 1 ); 203 return( 1 );
204 } 204 }
205 /* 205 /*
206 ********************************************** 206 ************************************************************
207 */ 207 */
208 static int MCGIDI_energy_parseGeneralEvaporati 208 static int MCGIDI_energy_parseGeneralEvaporationFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
209 209
210 double norm; 210 double norm;
211 xDataTOM_element *thetaTOM, *gTOM; 211 xDataTOM_element *thetaTOM, *gTOM;
212 ptwXYPoints *theta = NULL, *g = NULL; 212 ptwXYPoints *theta = NULL, *g = NULL;
213 char const *toUnits[2] = { "MeV", "MeV" }; 213 char const *toUnits[2] = { "MeV", "MeV" };
214 214
215 if( ( thetaTOM = xDataTOME_getOneElementBy 215 if( ( thetaTOM = xDataTOME_getOneElementByName( smr, element, "theta", 1 ) ) == NULL ) goto err;
216 if( ( theta = MCGIDI_misc_dataFromElement2 216 if( ( theta = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, thetaTOM, toUnits ) ) == NULL ) goto err;
217 217
218 if( ( gTOM = xDataTOME_getOneElementByName 218 if( ( gTOM = xDataTOME_getOneElementByName( smr, element, "g", 1 ) ) == NULL ) goto err;
219 toUnits[0] = ""; 219 toUnits[0] = "";
220 toUnits[1] = ""; 220 toUnits[1] = "";
221 if( ( g = MCGIDI_misc_dataFromElement2ptwX 221 if( ( g = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, gTOM, toUnits ) ) == NULL ) goto err;
222 if( MCGIDI_fromTOM_pdfOfX( smr, g, &(energ 222 if( MCGIDI_fromTOM_pdfOfX( smr, g, &(energy->g), &norm ) ) goto err;
223 energy->gInterpolation = ptwXY_getInterpol 223 energy->gInterpolation = ptwXY_getInterpolation( g );
224 g = ptwXY_free( g ); 224 g = ptwXY_free( g );
225 if( std::fabs( 1. - norm ) > 0.001 ) print 225 if( std::fabs( 1. - norm ) > 0.001 ) printf( "bad norm = %e\n", norm );
226 226
227 energy->type = MCGIDI_energyType_generalEv 227 energy->type = MCGIDI_energyType_generalEvaporation;
228 energy->theta = theta; 228 energy->theta = theta;
229 return( 0 ); 229 return( 0 );
230 230
231 err: 231 err:
232 if( theta != NULL ) ptwXY_free( theta ); 232 if( theta != NULL ) ptwXY_free( theta );
233 if( g != NULL ) ptwXY_free( g ); 233 if( g != NULL ) ptwXY_free( g );
234 return( 1 ); 234 return( 1 );
235 } 235 }
236 /* 236 /*
237 ********************************************** 237 ************************************************************
238 */ 238 */
239 static int MCGIDI_energy_parseSimpleMaxwellian 239 static int MCGIDI_energy_parseSimpleMaxwellianFissionFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
240 240
241 char const *U, *toUnits[2] = { "MeV", "MeV 241 char const *U, *toUnits[2] = { "MeV", "MeV" };
242 xDataTOM_element *thetaTOM; 242 xDataTOM_element *thetaTOM;
243 243
244 if( ( U = xDataTOM_getAttributesValueInEle 244 if( ( U = xDataTOM_getAttributesValueInElement( element, "U" ) ) == NULL ) {
245 smr_setReportError2( smr, smr_unknownI 245 smr_setReportError2( smr, smr_unknownID, 1, "functional form '%s' missing 'U' attribute", element->name );
246 goto err; 246 goto err;
247 } 247 }
248 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 248 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, U, "MeV", &(energy->U) ) != 0 ) goto err;
249 if( ( thetaTOM = xDataTOME_getOneElementBy 249 if( ( thetaTOM = xDataTOME_getOneElementByName( smr, element, "theta", 1 ) ) == NULL ) goto err;
250 if( ( energy->theta = MCGIDI_misc_dataFrom 250 if( ( energy->theta = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, thetaTOM, toUnits ) ) == NULL ) goto err;
251 energy->type = MCGIDI_energyType_simpleMax 251 energy->type = MCGIDI_energyType_simpleMaxwellianFission;
252 return( 0 ); 252 return( 0 );
253 253
254 err: 254 err:
255 return( 1 ); 255 return( 1 );
256 } 256 }
257 /* 257 /*
258 ********************************************** 258 ************************************************************
259 */ 259 */
260 static int MCGIDI_energy_parseEvaporationFromT 260 static int MCGIDI_energy_parseEvaporationFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
261 261
262 char const *U, *toUnits[2] = { "MeV", "MeV 262 char const *U, *toUnits[2] = { "MeV", "MeV" };
263 xDataTOM_element *thetaTOM; 263 xDataTOM_element *thetaTOM;
264 264
265 if( ( U = xDataTOM_getAttributesValueInEle 265 if( ( U = xDataTOM_getAttributesValueInElement( element, "U" ) ) == NULL ) {
266 smr_setReportError2( smr, smr_unknownI 266 smr_setReportError2( smr, smr_unknownID, 1, "functional form '%s' missing 'U' attribute", element->name );
267 goto err; 267 goto err;
268 } 268 }
269 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 269 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, U, "MeV", &(energy->U) ) != 0 ) goto err;
270 if( ( thetaTOM = xDataTOME_getOneElementBy 270 if( ( thetaTOM = xDataTOME_getOneElementByName( smr, element, "theta", 1 ) ) == NULL ) goto err;
271 if( ( energy->theta = MCGIDI_misc_dataFrom 271 if( ( energy->theta = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, thetaTOM, toUnits ) ) == NULL ) goto err;
272 energy->type = MCGIDI_energyType_evaporati 272 energy->type = MCGIDI_energyType_evaporation;
273 return( 0 ); 273 return( 0 );
274 274
275 err: 275 err:
276 return( 1 ); 276 return( 1 );
277 } 277 }
278 /* 278 /*
279 ********************************************** 279 ************************************************************
280 */ 280 */
281 static int MCGIDI_energy_parseWattFromTOM( sta 281 static int MCGIDI_energy_parseWattFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_energy *energy ) {
282 282
283 char const *U, *toUnits[2] = { "MeV", "MeV 283 char const *U, *toUnits[2] = { "MeV", "MeV" };
284 xDataTOM_element *aOrBTOM; 284 xDataTOM_element *aOrBTOM;
285 285
286 if( ( U = xDataTOM_getAttributesValueInEle 286 if( ( U = xDataTOM_getAttributesValueInElement( element, "U" ) ) == NULL ) {
287 smr_setReportError2( smr, smr_unknownI 287 smr_setReportError2( smr, smr_unknownID, 1, "functional form '%s' missing 'U' attribute", element->name );
288 goto err; 288 goto err;
289 } 289 }
290 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 290 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, U, "MeV", &(energy->U) ) != 0 ) goto err;
291 291
292 if( ( aOrBTOM = xDataTOME_getOneElementByN 292 if( ( aOrBTOM = xDataTOME_getOneElementByName( smr, element, "a", 1 ) ) == NULL ) goto err;
293 if( ( energy->Watt_a = MCGIDI_misc_dataFro 293 if( ( energy->Watt_a = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, aOrBTOM, toUnits ) ) == NULL ) goto err;
294 294
295 toUnits[1] = "1/MeV"; 295 toUnits[1] = "1/MeV";
296 if( ( aOrBTOM = xDataTOME_getOneElementByN 296 if( ( aOrBTOM = xDataTOME_getOneElementByName( smr, element, "b", 1 ) ) == NULL ) goto err;
297 if( ( energy->Watt_b = MCGIDI_misc_dataFro 297 if( ( energy->Watt_b = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, aOrBTOM, toUnits ) ) == NULL ) goto err;
298 298
299 energy->type = MCGIDI_energyType_Watt; 299 energy->type = MCGIDI_energyType_Watt;
300 return( 0 ); 300 return( 0 );
301 301
302 err: 302 err:
303 return( 1 ); 303 return( 1 );
304 } 304 }
305 /* 305 /*
306 ********************************************** 306 ************************************************************
307 */ 307 */
308 static int MCGIDI_energy_parseMadlandNixFromTO 308 static int MCGIDI_energy_parseMadlandNixFromTOM( statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy ) {
309 309
310 int iE, length, nXs, i1, n; 310 int iE, length, nXs, i1, n;
311 double E=0., T_M=0., EFL=0., EFH=0., argLi 311 double E=0., T_M=0., EFL=0., EFH=0., argList[3] = { 0., 0., 0. },
312 xs[] = { 1e-5, 1e-3, 1e-1, 1e1, 1e3 312 xs[] = { 1e-5, 1e-3, 1e-1, 1e1, 1e3, 1e5, 3e7 }, norm;
313 ptwXYPoints *ptwXY_TM = NULL, *pdfXY = NUL 313 ptwXYPoints *ptwXY_TM = NULL, *pdfXY = NULL;
314 ptwXYPoint *point; 314 ptwXYPoint *point;
315 ptwXPoints *cdfX = NULL; 315 ptwXPoints *cdfX = NULL;
316 nfu_status status = nfu_Okay; 316 nfu_status status = nfu_Okay;
317 xDataTOM_element *TM_TOM; 317 xDataTOM_element *TM_TOM;
318 xDataTOM_XYs *XYs; 318 xDataTOM_XYs *XYs;
319 MCGIDI_pdfsOfXGivenW *dists = &(energy->di 319 MCGIDI_pdfsOfXGivenW *dists = &(energy->dists);
320 MCGIDI_pdfOfX *dist; 320 MCGIDI_pdfOfX *dist;
321 char const *EF, *TMUnits[2] = { "MeV", "Me 321 char const *EF, *TMUnits[2] = { "MeV", "MeV" };
322 322
323 nXs = sizeof( xs ) / sizeof( xs[0] ); 323 nXs = sizeof( xs ) / sizeof( xs[0] );
324 324
325 if( ( EF = xDataTOM_getAttributesValueInEl 325 if( ( EF = xDataTOM_getAttributesValueInElement( functional, "EFL" ) ) == NULL ) {
326 smr_setReportError2( smr, smr_unknownI 326 smr_setReportError2( smr, smr_unknownID, 1, "MadlandNix '%s' missing 'EFL' attribute", functional->name );
327 goto err; 327 goto err;
328 } 328 }
329 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 329 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, EF, TMUnits[0], &EFL ) != 0 ) goto err;
330 argList[0] = EFL; 330 argList[0] = EFL;
331 331
332 if( ( EF = xDataTOM_getAttributesValueInEl 332 if( ( EF = xDataTOM_getAttributesValueInElement( functional, "EFH" ) ) == NULL ) {
333 smr_setReportError2( smr, smr_unknownI 333 smr_setReportError2( smr, smr_unknownID, 1, "MadlandNix '%s' missing 'EFH' attribute", functional->name );
334 goto err; 334 goto err;
335 } 335 }
336 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( 336 if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, EF, TMUnits[0], &EFH ) != 0 ) goto err;
337 argList[1] = EFH; 337 argList[1] = EFH;
338 338
339 if( ( TM_TOM = xDataTOME_getOneElementByNa 339 if( ( TM_TOM = xDataTOME_getOneElementByName( smr, functional, "T_M", 1 ) ) == NULL ) goto err;
340 if( ( XYs = (xDataTOM_XYs *) xDataTOME_get 340 if( ( XYs = (xDataTOM_XYs *) xDataTOME_getXDataIfID( smr, TM_TOM, "XYs" ) ) == NULL ) goto err;
341 if( ( ptwXY_TM = MCGIDI_misc_dataFromXYs2p 341 if( ( ptwXY_TM = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, ptwXY_interpolationLinLin, TMUnits ) ) == NULL ) goto err;
342 342
343 length = (int) ptwXY_length( ptwXY_TM ); 343 length = (int) ptwXY_length( ptwXY_TM );
344 dists->interpolationWY = ptwXY_interpolati 344 dists->interpolationWY = ptwXY_interpolationLinLin;
345 dists->interpolationXY = ptwXY_interpolati 345 dists->interpolationXY = ptwXY_interpolationLinLin; /* Ignoring what the data says as it is probably wrong. */
346 if( ( dists->Ws = (double *) smr_malloc2( 346 if( ( dists->Ws = (double *) smr_malloc2( smr, length * sizeof( double ), 1, "dists->Ws" ) ) == NULL ) goto err;
347 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_ 347 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, length * sizeof( MCGIDI_pdfOfX ), 0, "dists->dist" ) ) == NULL ) goto err;
348 348
349 for( iE = 0; iE < length; iE++ ) { 349 for( iE = 0; iE < length; iE++ ) {
350 ptwXY_getXYPairAtIndex( ptwXY_TM, iE, 350 ptwXY_getXYPairAtIndex( ptwXY_TM, iE, &E, &T_M );
351 argList[2] = T_M; 351 argList[2] = T_M;
352 dist = &(dists->dist[iE]); 352 dist = &(dists->dist[iE]);
353 dists->Ws[iE] = E; 353 dists->Ws[iE] = E;
354 354
355 if( ( pdfXY = ptwXY_createFromFunction 355 if( ( pdfXY = ptwXY_createFromFunction( nXs, xs, (ptwXY_createFromFunction_callback) MCGIDI_energy_parseMadlandNixFromTOM_callback,
356 (void *) argList, 1e-3, 0, 12, &st 356 (void *) argList, 1e-3, 0, 12, &status ) ) == NULL ) goto err;
357 if( ( status = ptwXY_normalize( pdfXY 357 if( ( status = ptwXY_normalize( pdfXY ) ) != nfu_Okay ) {
358 smr_setReportError2( smr, smr_unkn 358 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_normalize err = %d: %s\n", status, nfu_statusMessage( status ) );
359 goto err; 359 goto err;
360 } 360 }
361 361
362 if( ptwXY_simpleCoalescePoints( pdfXY 362 if( ptwXY_simpleCoalescePoints( pdfXY ) != nfu_Okay ) goto err;
363 dist->numberOfXs = n = (int) ptwXY_len 363 dist->numberOfXs = n = (int) ptwXY_length( pdfXY );
364 364
365 if( ( dist->Xs = (double *) smr_malloc 365 if( ( dist->Xs = (double *) smr_malloc2( smr, 3 * n * sizeof( double ), 0, "dist->Xs" ) ) == NULL ) goto err;
366 dists->numberOfWs++; 366 dists->numberOfWs++;
367 dist->pdf = &(dist->Xs[n]); 367 dist->pdf = &(dist->Xs[n]);
368 dist->cdf = &(dist->pdf[n]); 368 dist->cdf = &(dist->pdf[n]);
369 369
370 for( i1 = 0; i1 < n; i1++ ) { 370 for( i1 = 0; i1 < n; i1++ ) {
371 point = ptwXY_getPointAtIndex_Unsa 371 point = ptwXY_getPointAtIndex_Unsafely( pdfXY, i1 );
372 dist->Xs[i1] = point->x; 372 dist->Xs[i1] = point->x;
373 dist->pdf[i1] = point->y; 373 dist->pdf[i1] = point->y;
374 } 374 }
375 375
376 if( ( cdfX = ptwXY_runningIntegral( pd 376 if( ( cdfX = ptwXY_runningIntegral( pdfXY, &status ) ) == NULL ) {
377 smr_setReportError2( smr, smr_unkn 377 smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_runningIntegral err = %d: %s\n", status, nfu_statusMessage( status ) );
378 goto err; 378 goto err;
379 } 379 }
380 380
381 norm = ptwX_getPointAtIndex_Unsafely( 381 norm = ptwX_getPointAtIndex_Unsafely( cdfX, n - 1 );
382 for( i1 = 0; i1 < n; i1++ ) dist->cdf[ 382 for( i1 = 0; i1 < n; i1++ ) dist->cdf[i1] = ptwX_getPointAtIndex_Unsafely( cdfX, i1 ) / norm;
383 for( i1 = 0; i1 < n; i1++ ) dist->pdf[ 383 for( i1 = 0; i1 < n; i1++ ) dist->pdf[i1] /= norm;
384 pdfXY = ptwXY_free( pdfXY ); 384 pdfXY = ptwXY_free( pdfXY );
385 cdfX = ptwX_free( cdfX ); 385 cdfX = ptwX_free( cdfX );
386 } 386 }
387 387
388 energy->type = MCGIDI_energyType_MadlandNi 388 energy->type = MCGIDI_energyType_MadlandNix;
389 389
390 ptwXY_free( ptwXY_TM ); 390 ptwXY_free( ptwXY_TM );
391 return( 0 ); 391 return( 0 );
392 392
393 err: 393 err:
394 if( ptwXY_TM != NULL ) ptwXY_free( ptwXY_T 394 if( ptwXY_TM != NULL ) ptwXY_free( ptwXY_TM );
395 if( pdfXY != NULL ) ptwXY_free( pdfXY ); 395 if( pdfXY != NULL ) ptwXY_free( pdfXY );
396 if( cdfX != NULL ) cdfX = ptwX_free( cdfX 396 if( cdfX != NULL ) cdfX = ptwX_free( cdfX );
397 397
398 return( 1 ); 398 return( 1 );
399 } 399 }
400 /* 400 /*
401 ********************************************** 401 ************************************************************
402 */ 402 */
403 static nfu_status MCGIDI_energy_parseMadlandNi 403 static nfu_status MCGIDI_energy_parseMadlandNixFromTOM_callback( double Ep, double *y, void *argList ) {
404 404
405 double *parameters = (double *) argList, E 405 double *parameters = (double *) argList, EFL, EFH, T_M;
406 nfu_status status = nfu_Okay; 406 nfu_status status = nfu_Okay;
407 407
408 EFL = parameters[0]; 408 EFL = parameters[0];
409 EFH = parameters[1]; 409 EFH = parameters[1];
410 T_M = parameters[2]; 410 T_M = parameters[2];
411 *y = MCGIDI_energy_parseMadlandNixFromTOM_ 411 *y = MCGIDI_energy_parseMadlandNixFromTOM_callback_g( Ep, EFL, T_M, &status );
412 if( status == nfu_Okay ) *y += MCGIDI_ener 412 if( status == nfu_Okay ) *y += MCGIDI_energy_parseMadlandNixFromTOM_callback_g( Ep, EFH, T_M, &status );
413 *y *= 0.5; 413 *y *= 0.5;
414 return( status ); 414 return( status );
415 } 415 }
416 /* 416 /*
417 ********************************************** 417 ************************************************************
418 */ 418 */
419 static double MCGIDI_energy_parseMadlandNixFro 419 static double MCGIDI_energy_parseMadlandNixFromTOM_callback_g( double Ep, double E_F, double T_M, nfu_status *status ) {
420 420
421 double u1, u2, E1, E2 = 0., gamma1 = 0., g 421 double u1, u2, E1, E2 = 0., gamma1 = 0., gamma2 = 0., signG = 1;
422 422
423 u1 = std::sqrt( Ep ) - std::sqrt( E_F ); 423 u1 = std::sqrt( Ep ) - std::sqrt( E_F );
424 u1 *= u1 / T_M; 424 u1 *= u1 / T_M;
425 u2 = std::sqrt( Ep ) + std::sqrt( E_F ); 425 u2 = std::sqrt( Ep ) + std::sqrt( E_F );
426 u2 *= u2 / T_M; 426 u2 *= u2 / T_M;
427 E1 = 0; /* u1^3/2 * E 427 E1 = 0; /* u1^3/2 * E1 is zero for u1 = 0. but E1 is infinity, whence, the next test. */
428 if( u1 != 0 ) E1 = nf_exponentialIntegral( 428 if( u1 != 0 ) E1 = nf_exponentialIntegral( 1, u1, status );
429 if( *status == nfu_Okay ) E2 = nf_exponent 429 if( *status == nfu_Okay ) E2 = nf_exponentialIntegral( 1, u2, status );
430 if( *status != nfu_Okay ) return( 0. ); 430 if( *status != nfu_Okay ) return( 0. );
431 if( u1 > 2. ) { 431 if( u1 > 2. ) {
432 signG = -1; 432 signG = -1;
433 gamma1 = nf_incompleteGammaFunctionCom 433 gamma1 = nf_incompleteGammaFunctionComplementary( 1.5, u1, status );
434 if( *status == nfu_Okay ) gamma2 = nf_ 434 if( *status == nfu_Okay ) gamma2 = nf_incompleteGammaFunctionComplementary( 1.5, u2, status ); }
435 else { 435 else {
436 gamma1 = nf_incompleteGammaFunction( 1 436 gamma1 = nf_incompleteGammaFunction( 1.5, u1, status );
437 if( *status == nfu_Okay ) gamma2 = nf_ 437 if( *status == nfu_Okay ) gamma2 = nf_incompleteGammaFunction( 1.5, u2, status );
438 } 438 }
439 if( *status != nfu_Okay ) return( 0. ); 439 if( *status != nfu_Okay ) return( 0. );
440 return( ( u2 * std::sqrt( u2 ) * E2 - u1 * 440 return( ( u2 * std::sqrt( u2 ) * E2 - u1 * std::sqrt( u1 ) * E1 + signG * ( gamma2 - gamma1 ) ) / ( 3 * std::sqrt( E_F * T_M ) ) );
441 } 441 }
442 /* 442 /*
443 ********************************************** 443 ************************************************************
444 */ 444 */
445 static int MCGIDI_energy_parseNBodyPhaseSpaceF 445 static int MCGIDI_energy_parseNBodyPhaseSpaceFromTOM( statusMessageReporting *smr, xDataTOM_element *functional, MCGIDI_energy *energy,
446 MCGIDI_distribution *distribution ) { 446 MCGIDI_distribution *distribution ) {
447 447
448 int argList[1]; 448 int argList[1];
449 double xs[2] = { 0.0, 1.0 }, productMass_M 449 double xs[2] = { 0.0, 1.0 }, productMass_MeV, norm;
450 ptwXYPoints *pdf = NULL; 450 ptwXYPoints *pdf = NULL;
451 nfu_status status; 451 nfu_status status;
452 char const *mass; 452 char const *mass;
453 453
454 if( xDataTOME_convertAttributeToInteger( N 454 if( xDataTOME_convertAttributeToInteger( NULL, functional, "numberOfProducts", &(energy->NBodyPhaseSpace.numberOfProducts) ) != 0 ) goto err;
455 if( ( mass = xDataTOM_getAttributesValueIn 455 if( ( mass = xDataTOM_getAttributesValueInElement( functional, "mass" ) ) == NULL ) {
456 smr_setReportError2( smr, smr_unknownI 456 smr_setReportError2( smr, smr_unknownID, 1, "functional form '%s' missing 'mass' attribute", functional->name );
457 goto err; 457 goto err;
458 } 458 }
459 if( MCGIDI_misc_PQUStringToDouble( smr, ma 459 if( MCGIDI_misc_PQUStringToDouble( smr, mass, "amu", MCGIDI_AMU2MeV, &(energy->NBodyPhaseSpace.mass) ) ) goto err;
460 argList[0] = energy->NBodyPhaseSpace.numbe 460 argList[0] = energy->NBodyPhaseSpace.numberOfProducts;
461 if( ( pdf = ptwXY_createFromFunction( 2, x 461 if( ( pdf = ptwXY_createFromFunction( 2, xs, MCGIDI_energy_NBodyPhaseSpacePDF_callback, (void *) argList, 1e-3, 0, 16, &status ) ) == NULL ) {
462 smr_setReportError2( smr, smr_unknownI 462 smr_setReportError2( smr, smr_unknownID, 1, "creating NBodyPhaseSpace pdf failed with ptwXY_createFromFunction error = %d (%s)",
463 status, nfu_statusMessage( status 463 status, nfu_statusMessage( status ) );
464 goto err; 464 goto err;
465 } 465 }
466 if( MCGIDI_fromTOM_pdfOfX( smr, pdf, &(ene 466 if( MCGIDI_fromTOM_pdfOfX( smr, pdf, &(energy->g), &norm ) ) goto err;
467 productMass_MeV = MCGIDI_product_getMass_M 467 productMass_MeV = MCGIDI_product_getMass_MeV( smr, distribution->product );
468 if( !smr_isOk( smr ) ) goto err; 468 if( !smr_isOk( smr ) ) goto err;
469 energy->NBodyPhaseSpace.massFactor = ( 1. 469 energy->NBodyPhaseSpace.massFactor = ( 1. - productMass_MeV / ( MCGIDI_AMU2MeV * energy->NBodyPhaseSpace.mass ) ); /* ??????? Hardwired MCGIDI_AMU2MeV */
470 energy->NBodyPhaseSpace.Q_MeV = MCGIDI_out 470 energy->NBodyPhaseSpace.Q_MeV = MCGIDI_outputChannel_getQ_MeV( smr, distribution->product->outputChannel, 0. );
471 if( !smr_isOk( smr ) ) goto err; 471 if( !smr_isOk( smr ) ) goto err;
472 472
473 ptwXY_free( pdf ); 473 ptwXY_free( pdf );
474 energy->type = MCGIDI_energyType_NBodyPhas 474 energy->type = MCGIDI_energyType_NBodyPhaseSpace;
475 475
476 return( 0 ); 476 return( 0 );
477 477
478 err: 478 err:
479 if( pdf != NULL ) ptwXY_free( pdf ); 479 if( pdf != NULL ) ptwXY_free( pdf );
480 return( 1 ); 480 return( 1 );
481 } 481 }
482 /* 482 /*
483 ********************************************** 483 ************************************************************
484 */ 484 */
485 static nfu_status MCGIDI_energy_NBodyPhaseSpac 485 static nfu_status MCGIDI_energy_NBodyPhaseSpacePDF_callback( double x, double *y, void *argList ) {
486 486
487 int numberOfProducts = *((int *) argList); 487 int numberOfProducts = *((int *) argList);
488 double e = 0.5 * ( 3 * numberOfProducts - 488 double e = 0.5 * ( 3 * numberOfProducts - 8 );
489 489
490 *y = std::sqrt( x ) * G4Pow::GetInstance() 490 *y = std::sqrt( x ) * G4Pow::GetInstance()->powA( 1.0 - x, e );
491 return( nfu_Okay ); 491 return( nfu_Okay );
492 } 492 }
493 /* 493 /*
494 ********************************************** 494 ************************************************************
495 */ 495 */
496 int MCGIDI_energy_sampleEnergy( statusMessageR 496 int MCGIDI_energy_sampleEnergy( statusMessageReporting *smr, MCGIDI_energy *energy, MCGIDI_quantitiesLookupModes &modes,
497 MCGIDI_decaySamplingInfo *decaySamplin 497 MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
498 /* 498 /*
499 * This function must be called before angula 499 * This function must be called before angular sampling as it sets the frame but does not test it.
500 */ 500 */
501 double theta, randomEp, Watt_a, Watt_b, e_ 501 double theta, randomEp, Watt_a, Watt_b, e_in = modes.getProjectileEnergy( );
502 MCGIDI_pdfsOfXGivenW_sampled sampled; 502 MCGIDI_pdfsOfXGivenW_sampled sampled;
503 503
504 decaySamplingInfo->frame = energy->frame; 504 decaySamplingInfo->frame = energy->frame;
505 switch( energy->type ) { 505 switch( energy->type ) {
506 case MCGIDI_energyType_primaryGamma : 506 case MCGIDI_energyType_primaryGamma :
507 decaySamplingInfo->Ep = energy->gammaE 507 decaySamplingInfo->Ep = energy->gammaEnergy_MeV + e_in * energy->primaryGammaMassFactor;
508 break; 508 break;
509 case MCGIDI_energyType_discreteGamma : 509 case MCGIDI_energyType_discreteGamma :
510 decaySamplingInfo->Ep = energy->gammaE 510 decaySamplingInfo->Ep = energy->gammaEnergy_MeV;
511 break; 511 break;
512 case MCGIDI_energyType_linear : 512 case MCGIDI_energyType_linear :
513 randomEp = decaySamplingInfo->rng( dec 513 randomEp = decaySamplingInfo->rng( decaySamplingInfo->rngState );
514 sampled.smr = smr; 514 sampled.smr = smr;
515 sampled.w = e_in; 515 sampled.w = e_in;
516 MCGIDI_sampling_sampleX_from_pdfsOfXGi 516 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(energy->dists), &sampled, randomEp );
517 decaySamplingInfo->Ep = sampled.x; 517 decaySamplingInfo->Ep = sampled.x;
518 break; 518 break;
519 case MCGIDI_energyType_generalEvaporation 519 case MCGIDI_energyType_generalEvaporation :
520 sampled.interpolationXY = energy->gInt 520 sampled.interpolationXY = energy->gInterpolation;
521 MCGIDI_sampling_sampleX_from_pdfOfX( & 521 MCGIDI_sampling_sampleX_from_pdfOfX( &(energy->g), &sampled, decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
522 theta = MCGIDI_sampling_ptwXY_getValue 522 theta = MCGIDI_sampling_ptwXY_getValueAtX( energy->theta, e_in );
523 decaySamplingInfo->Ep = theta * sample 523 decaySamplingInfo->Ep = theta * sampled.x;
524 break; 524 break;
525 case MCGIDI_energyType_simpleMaxwellianFis 525 case MCGIDI_energyType_simpleMaxwellianFission :
526 theta = MCGIDI_sampling_ptwXY_getValue 526 theta = MCGIDI_sampling_ptwXY_getValueAtX( energy->theta, e_in );
527 MCGIDI_energy_sampleSimpleMaxwellianFi 527 MCGIDI_energy_sampleSimpleMaxwellianFission( smr, ( e_in - energy->U ) / theta, decaySamplingInfo );
528 decaySamplingInfo->Ep *= theta; 528 decaySamplingInfo->Ep *= theta;
529 break; 529 break;
530 case MCGIDI_energyType_evaporation : 530 case MCGIDI_energyType_evaporation :
531 theta = MCGIDI_sampling_ptwXY_getValue 531 theta = MCGIDI_sampling_ptwXY_getValueAtX( energy->theta, e_in );
532 MCGIDI_energy_sampleEvaporation( smr, 532 MCGIDI_energy_sampleEvaporation( smr, ( e_in - energy->U ) / theta, decaySamplingInfo );
533 decaySamplingInfo->Ep *= theta; 533 decaySamplingInfo->Ep *= theta;
534 break; 534 break;
535 case MCGIDI_energyType_Watt : 535 case MCGIDI_energyType_Watt :
536 Watt_a = MCGIDI_sampling_ptwXY_getValu 536 Watt_a = MCGIDI_sampling_ptwXY_getValueAtX( energy->Watt_a, e_in );
537 Watt_b = MCGIDI_sampling_ptwXY_getValu 537 Watt_b = MCGIDI_sampling_ptwXY_getValueAtX( energy->Watt_b, e_in );
538 MCGIDI_energy_sampleWatt( smr, e_in - 538 MCGIDI_energy_sampleWatt( smr, e_in - energy->U, Watt_a, Watt_b, decaySamplingInfo );
539 break; 539 break;
540 case MCGIDI_energyType_MadlandNix : 540 case MCGIDI_energyType_MadlandNix :
541 MCGIDI_sampling_sampleX_from_pdfsOfXGi 541 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(energy->dists), &sampled, decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
542 decaySamplingInfo->Ep = sampled.x; 542 decaySamplingInfo->Ep = sampled.x;
543 break; 543 break;
544 case MCGIDI_energyType_NBodyPhaseSpace : 544 case MCGIDI_energyType_NBodyPhaseSpace :
545 MCGIDI_sampling_sampleX_from_pdfOfX( & 545 MCGIDI_sampling_sampleX_from_pdfOfX( &(energy->g), &sampled, decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
546 decaySamplingInfo->Ep = ( energy->e_in 546 decaySamplingInfo->Ep = ( energy->e_inCOMFactor * e_in + energy->NBodyPhaseSpace.Q_MeV ) * energy->NBodyPhaseSpace.massFactor * sampled.x;
547 break; 547 break;
548 case MCGIDI_energyType_weightedFunctional 548 case MCGIDI_energyType_weightedFunctional :
549 MCGIDI_energy_sampleWeightedFunctional 549 MCGIDI_energy_sampleWeightedFunctional( smr, energy, modes, decaySamplingInfo );
550 break; 550 break;
551 default : 551 default :
552 smr_setReportError2( smr, smr_unknownI 552 smr_setReportError2( smr, smr_unknownID, 1, "energy type = %d not supported", energy->type );
553 } 553 }
554 554
555 return( !smr_isOk( smr ) ); 555 return( !smr_isOk( smr ) );
556 } 556 }
557 /* 557 /*
558 ********************************************** 558 ************************************************************
559 */ 559 */
560 static int MCGIDI_energy_sampleSimpleMaxwellia 560 static int MCGIDI_energy_sampleSimpleMaxwellianFission( statusMessageReporting * /*smr*/, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
561 561
562 int i1; 562 int i1;
563 double a = e_in_U_theta, b, c, x, norm_a, 563 double a = e_in_U_theta, b, c, x, norm_a, xMin = 0., xMax = a, sqrt_x, sqrt_pi_2 = std::sqrt( M_PI ) / 2.;
564 564
565 sqrt_x = std::sqrt( a ); 565 sqrt_x = std::sqrt( a );
566 norm_a = sqrt_pi_2 * erf( sqrt_x ) - sqrt_ 566 norm_a = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x * G4Exp( -a );
567 b = norm_a * decaySamplingInfo->rng( decay 567 b = norm_a * decaySamplingInfo->rng( decaySamplingInfo->rngState );
568 for( i1 = 0; i1 < 16; i1++ ) { 568 for( i1 = 0; i1 < 16; i1++ ) {
569 x = 0.5 * ( xMin + xMax ); 569 x = 0.5 * ( xMin + xMax );
570 sqrt_x = std::sqrt( x ); 570 sqrt_x = std::sqrt( x );
571 c = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x 571 c = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x * G4Exp( -x );
572 if( b < c ) { 572 if( b < c ) {
573 xMax = x; } 573 xMax = x; }
574 else { 574 else {
575 xMin = x; 575 xMin = x;
576 } 576 }
577 } 577 }
578 /* To order e, the correct x is x + e 578 /* To order e, the correct x is x + e where e = 1 + ( 1 - b * exp( x ) ) / x. */
579 decaySamplingInfo->Ep = x; 579 decaySamplingInfo->Ep = x;
580 580
581 return( 0 ); 581 return( 0 );
582 } 582 }
583 /* 583 /*
584 ********************************************** 584 ************************************************************
585 */ 585 */
586 static int MCGIDI_energy_sampleEvaporation( st 586 static int MCGIDI_energy_sampleEvaporation( statusMessageReporting * /*smr*/, double e_in_U_theta, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
587 587
588 int i1; 588 int i1;
589 double a = e_in_U_theta, b, c, x, norm_a, 589 double a = e_in_U_theta, b, c, x, norm_a, xMin = 0., xMax = a;
590 590
591 norm_a = 1 - ( 1 + a ) * G4Exp( -a ); 591 norm_a = 1 - ( 1 + a ) * G4Exp( -a );
592 b = 1. - norm_a * decaySamplingInfo->rng( 592 b = 1. - norm_a * decaySamplingInfo->rng( decaySamplingInfo->rngState );
593 for( i1 = 0; i1 < 16; i1++ ) { 593 for( i1 = 0; i1 < 16; i1++ ) {
594 x = 0.5 * ( xMin + xMax ); 594 x = 0.5 * ( xMin + xMax );
595 c = ( 1 + x ) * G4Exp( -x ); 595 c = ( 1 + x ) * G4Exp( -x );
596 if( b > c ) { 596 if( b > c ) {
597 xMax = x; } 597 xMax = x; }
598 else { 598 else {
599 xMin = x; 599 xMin = x;
600 } 600 }
601 } 601 }
602 /* To order e, the correct x is x + e 602 /* To order e, the correct x is x + e where e = 1 + ( 1 - b * std::exp( x ) ) / x. */
603 decaySamplingInfo->Ep = x; 603 decaySamplingInfo->Ep = x;
604 604
605 return( 0 ); 605 return( 0 );
606 } 606 }
607 /* 607 /*
608 ********************************************** 608 ************************************************************
609 */ 609 */
610 static int MCGIDI_energy_sampleWatt( statusMes 610 static int MCGIDI_energy_sampleWatt( statusMessageReporting * /*smr*/, double e_in_U, double Watt_a, double Watt_b, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
611 /* 611 /*
612 * From MCAPM via Sample Watt Spectrum as in 612 * From MCAPM via Sample Watt Spectrum as in TART ( Kalos algorithm ).
613 */ 613 */
614 double WattMin = 0., WattMax = e_in_U, x, 614 double WattMin = 0., WattMax = e_in_U, x, y, z, energyOut = 0., rand1, rand2;
615 615
616 x = 1. + ( Watt_b / ( 8. * Watt_a ) ); 616 x = 1. + ( Watt_b / ( 8. * Watt_a ) );
617 y = ( x + std::sqrt( x * x - 1. ) ) / Watt 617 y = ( x + std::sqrt( x * x - 1. ) ) / Watt_a;
618 z = Watt_a * y - 1.; 618 z = Watt_a * y - 1.;
619 G4int icounter=0; 619 G4int icounter=0;
620 G4int icounter_max=1024; 620 G4int icounter_max=1024;
621 do 621 do
622 { 622 {
623 icounter++; 623 icounter++;
624 if ( icounter > icounter_max ) { 624 if ( icounter > icounter_max ) {
625 G4cout << "Loop-counter exceeded the 625 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
626 break; 626 break;
627 } 627 }
628 rand1 = -G4Log( decaySamplingInfo->rng 628 rand1 = -G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
629 rand2 = -G4Log( decaySamplingInfo->rng 629 rand2 = -G4Log( decaySamplingInfo->rng( decaySamplingInfo->rngState ) );
630 energyOut = y * rand1; 630 energyOut = y * rand1;
631 } 631 }
632 while( ( ( rand2 - z * ( rand1 + 1. ) ) * 632 while( ( ( rand2 - z * ( rand1 + 1. ) ) * ( rand2 - z * ( rand1 + 1. ) ) > Watt_b * y * rand1 ) || ( energyOut < WattMin ) || ( energyOut > WattMax ) ); // Loop checking, 11.06.2015, T. Koi
633 decaySamplingInfo->Ep = energyOut; 633 decaySamplingInfo->Ep = energyOut;
634 634
635 return( 0 ); 635 return( 0 );
636 } 636 }
637 /* 637 /*
638 ********************************************** 638 ************************************************************
639 */ 639 */
640 static int MCGIDI_energy_sampleWeightedFunctio 640 static int MCGIDI_energy_sampleWeightedFunctional( statusMessageReporting *smr, MCGIDI_energy *energy,
641 MCGIDI_quantitiesLookupModes &modes, M 641 MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
642 /* 642 /*
643 c This routine assumes that the weights sum 643 c This routine assumes that the weights sum to 1.
644 */ 644 */
645 int iW; 645 int iW;
646 double rW = decaySamplingInfo->rng( decayS 646 double rW = decaySamplingInfo->rng( decaySamplingInfo->rngState ), cumulativeW = 0., weight;
647 MCGIDI_energyWeightedFunctional *weightedF 647 MCGIDI_energyWeightedFunctional *weightedFunctional = NULL;
648 648
649 for( iW = 0; iW < energy->weightedFunction 649 for( iW = 0; iW < energy->weightedFunctionals.numberOfWeights; iW++ ) {
650 weightedFunctional = &(energy->weighte 650 weightedFunctional = &(energy->weightedFunctionals.weightedFunctional[iW]);
651 weight = MCGIDI_sampling_ptwXY_getValu 651 weight = MCGIDI_sampling_ptwXY_getValueAtX( weightedFunctional->weight, modes.getProjectileEnergy( ) );
652 cumulativeW += weight; 652 cumulativeW += weight;
653 if( cumulativeW >= rW ) break; 653 if( cumulativeW >= rW ) break;
654 } 654 }
655 return( MCGIDI_energy_sampleEnergy( smr, w 655 return( MCGIDI_energy_sampleEnergy( smr, weightedFunctional->energy, modes, decaySamplingInfo ) );
656 } 656 }
657 657
658 #if defined __cplusplus 658 #if defined __cplusplus
659 } 659 }
660 #endif 660 #endif
661 661
662 662