| 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 <PoPs.h> 8 #include <PoPs.h>
9 #include "MCGIDI.h" 9 #include "MCGIDI.h"
10 #include "MCGIDI_misc.h" 10 #include "MCGIDI_misc.h"
11 #include "MCGIDI_private.h" 11 #include "MCGIDI_private.h"
12 12
13 #if defined __cplusplus 13 #if defined __cplusplus
14 namespace GIDI { 14 namespace GIDI {
15 using namespace GIDI; 15 using namespace GIDI;
16 #endif 16 #endif
17 17
18 #define nParticleChanges 6 18 #define nParticleChanges 6
19 19
20 static int MCGIDI_reaction_initialize2( status 20 static int MCGIDI_reaction_initialize2( statusMessageReporting *smr, MCGIDI_reaction *reaction );
21 static int MCGIDI_reaction_particleChanges( MC 21 static int MCGIDI_reaction_particleChanges( MCGIDI_POP *projectile, MCGIDI_POP *target, MCGIDI_productsInfo *productsInfo, int n1, int *particlesChanges );
22 static int MCGIDI_reaction_ParseReactionTypeAn 22 static int MCGIDI_reaction_ParseReactionTypeAndDetermineProducts( statusMessageReporting *smr, MCGIDI_POPs *pops, MCGIDI_reaction *reaction );
23 static int MCGIDI_reaction_ParseDetermineReact 23 static int MCGIDI_reaction_ParseDetermineReactionProducts( statusMessageReporting *smr, MCGIDI_POPs *pops, MCGIDI_outputChannel *outputChannel,
24 MCGIDI_productsInfo *productsInfo, MCGIDI_ 24 MCGIDI_productsInfo *productsInfo, MCGIDI_reaction *reaction, double *finalQ, int level );
25 static int MCGIDI_reaction_addReturnProduct( s 25 static int MCGIDI_reaction_addReturnProduct( statusMessageReporting *smr, MCGIDI_productsInfo *productsInfo, int ID, MCGIDI_product *product,
26 MCGIDI_reaction *reaction, int transpo 26 MCGIDI_reaction *reaction, int transportable );
27 static int MCGIDI_reaction_setENDL_CSNumbers( 27 static int MCGIDI_reaction_setENDL_CSNumbers( statusMessageReporting *smr, MCGIDI_reaction *reaction );
28 /* 28 /*
29 ********************************************** 29 ************************************************************
30 */ 30 */
31 MCGIDI_reaction *MCGIDI_reaction_new( statusMe 31 MCGIDI_reaction *MCGIDI_reaction_new( statusMessageReporting *smr ) {
32 32
33 MCGIDI_reaction *reaction; 33 MCGIDI_reaction *reaction;
34 34
35 if( ( reaction = (MCGIDI_reaction *) smr_m 35 if( ( reaction = (MCGIDI_reaction *) smr_malloc2( smr, sizeof( MCGIDI_reaction ), 0, "reaction" ) ) == NULL ) return( NULL );
36 if( MCGIDI_reaction_initialize( smr, react 36 if( MCGIDI_reaction_initialize( smr, reaction ) ) reaction = MCGIDI_reaction_free( smr, reaction );
37 return( reaction ); 37 return( reaction );
38 } 38 }
39 /* 39 /*
40 ********************************************** 40 ************************************************************
41 */ 41 */
42 int MCGIDI_reaction_initialize( statusMessageR 42 int MCGIDI_reaction_initialize( statusMessageReporting *smr, MCGIDI_reaction *reaction ) {
43 43
44 if( MCGIDI_reaction_initialize2( smr, reac 44 if( MCGIDI_reaction_initialize2( smr, reaction ) != 0 ) return( 1 );
45 reaction->transportabilities = new transpo 45 reaction->transportabilities = new transportabilitiesMap( );
46 return( 0 ); 46 return( 0 );
47 } 47 }
48 /* 48 /*
49 ********************************************** 49 ************************************************************
50 */ 50 */
51 static int MCGIDI_reaction_initialize2( status 51 static int MCGIDI_reaction_initialize2( statusMessageReporting *smr, MCGIDI_reaction *reaction ) {
52 52
53 memset( reaction, 0, sizeof( MCGIDI_reacti 53 memset( reaction, 0, sizeof( MCGIDI_reaction ) );
54 xDataTOMAL_initial( smr, &(reaction->attri 54 xDataTOMAL_initial( smr, &(reaction->attributes) );
55 return( 0 ); 55 return( 0 );
56 } 56 }
57 /* 57 /*
58 ********************************************** 58 ************************************************************
59 */ 59 */
60 MCGIDI_reaction *MCGIDI_reaction_free( statusM 60 MCGIDI_reaction *MCGIDI_reaction_free( statusMessageReporting *smr, MCGIDI_reaction *reaction ) {
61 61
62 MCGIDI_reaction_release( smr, reaction ); 62 MCGIDI_reaction_release( smr, reaction );
63 smr_freeMemory( (void **) &reaction ); 63 smr_freeMemory( (void **) &reaction );
64 return( NULL ); 64 return( NULL );
65 } 65 }
66 /* 66 /*
67 ********************************************** 67 ************************************************************
68 */ 68 */
69 int MCGIDI_reaction_release( statusMessageRepo 69 int MCGIDI_reaction_release( statusMessageReporting *smr, MCGIDI_reaction *reaction ) {
70 70
71 ptwXY_free( reaction->crossSection ); 71 ptwXY_free( reaction->crossSection );
72 ptwX_free( reaction->crossSectionGrouped ) 72 ptwX_free( reaction->crossSectionGrouped );
73 MCGIDI_outputChannel_release( smr, &(react 73 MCGIDI_outputChannel_release( smr, &(reaction->outputChannel) );
74 xDataTOMAL_release( &(reaction->attributes 74 xDataTOMAL_release( &(reaction->attributes) );
75 smr_freeMemory( (void **) &(reaction->outp 75 smr_freeMemory( (void **) &(reaction->outputChannelStr) );
76 if( reaction->productsInfo.productInfo != 76 if( reaction->productsInfo.productInfo != NULL ) smr_freeMemory( (void **) &(reaction->productsInfo.productInfo) );
77 delete reaction->transportabilities; 77 delete reaction->transportabilities;
78 MCGIDI_reaction_initialize2( smr, reaction 78 MCGIDI_reaction_initialize2( smr, reaction );
79 return( 0 ); 79 return( 0 );
80 } 80 }
81 /* 81 /*
82 ********************************************** 82 ************************************************************
83 */ 83 */
84 int MCGIDI_reaction_parseFromTOM( statusMessag 84 int MCGIDI_reaction_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_target_heated *target,
85 MCGIDI_POPs *pops, MCGIDI_reaction *re 85 MCGIDI_POPs *pops, MCGIDI_reaction *reaction ) {
86 86
87 xDataTOM_element *child, *linear, *outputC 87 xDataTOM_element *child, *linear, *outputChannel;
88 enum xDataTOM_interpolationFlag independen 88 enum xDataTOM_interpolationFlag independent, dependent;
89 enum xDataTOM_interpolationQualifier quali 89 enum xDataTOM_interpolationQualifier qualifier;
90 char const *outputChannelStr, *crossSectio 90 char const *outputChannelStr, *crossSectionUnits[2] = { "MeV", "b" };
91 91
92 MCGIDI_reaction_initialize( smr, reaction 92 MCGIDI_reaction_initialize( smr, reaction );
93 93
94 reaction->target = target; 94 reaction->target = target;
95 reaction->reactionType = MCGIDI_reactionTy 95 reaction->reactionType = MCGIDI_reactionType_unknown_e;
96 if( xDataTOME_copyAttributionList( smr, &( 96 if( xDataTOME_copyAttributionList( smr, &(reaction->attributes), element ) ) goto err;
97 if( xDataTOME_convertAttributeToInteger( s 97 if( xDataTOME_convertAttributeToInteger( smr, element, "ENDF_MT", &(reaction->ENDF_MT) ) ) goto err;
98 if( ( outputChannelStr = xDataTOM_getAttri 98 if( ( outputChannelStr = xDataTOM_getAttributesValueInElement( element, "outputChannel" ) ) == NULL ) goto err;
99 if( ( reaction->outputChannelStr = smr_all 99 if( ( reaction->outputChannelStr = smr_allocateCopyString2( smr, outputChannelStr, "reaction->outputChannelStr" ) ) == NULL ) goto err;
100 100
101 if( ( child = xDataTOME_getOneElementByNam 101 if( ( child = xDataTOME_getOneElementByName( smr, element, "crossSection", 1 ) ) == NULL ) goto err;
102 if( ( linear = xDataTOME_getOneElementByNa 102 if( ( linear = xDataTOME_getOneElementByName( smr, child, "linear", 0 ) ) == NULL ) {
103 if( ( linear = xDataTOME_getOneElement 103 if( ( linear = xDataTOME_getOneElementByName( smr, child, "pointwise", 1 ) ) == NULL ) goto err;
104 } 104 }
105 if( xDataTOME_getInterpolation( smr, linea 105 if( xDataTOME_getInterpolation( smr, linear, 0, &independent, &dependent, &qualifier ) ) goto err;
106 if( ( independent != xDataTOM_interpolatio 106 if( ( independent != xDataTOM_interpolationFlag_linear ) || ( dependent != xDataTOM_interpolationFlag_linear ) ) {
107 smr_setReportError2( smr, smr_unknownI 107 smr_setReportError2( smr, smr_unknownID, 1, "cross section interpolation (%d,%d) is not linear-linear", independent, dependent );
108 goto err; 108 goto err;
109 } 109 }
110 if( ( reaction->crossSection = MCGIDI_misc 110 if( ( reaction->crossSection = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, linear, crossSectionUnits ) ) == NULL ) goto err;
111 reaction->domainValuesPresent = 1; 111 reaction->domainValuesPresent = 1;
112 reaction->EMin = ptwXY_getXMin( reaction-> 112 reaction->EMin = ptwXY_getXMin( reaction->crossSection );
113 reaction->EMax = ptwXY_getXMax( reaction-> 113 reaction->EMax = ptwXY_getXMax( reaction->crossSection );
114 114
115 if( ( outputChannel = xDataTOME_getOneElem 115 if( ( outputChannel = xDataTOME_getOneElementByName( smr, element, "outputChannel", 1 ) ) == NULL ) goto err;
116 if( MCGIDI_outputChannel_parseFromTOM( smr 116 if( MCGIDI_outputChannel_parseFromTOM( smr, outputChannel, pops, &(reaction->outputChannel), reaction, NULL ) ) goto err;
117 117
118 if( MCGIDI_reaction_ParseReactionTypeAndDe 118 if( MCGIDI_reaction_ParseReactionTypeAndDetermineProducts( smr, pops, reaction ) != 0 ) goto err;
119 119
120 return( 0 ); 120 return( 0 );
121 121
122 err: 122 err:
123 MCGIDI_reaction_release( smr, reaction ); 123 MCGIDI_reaction_release( smr, reaction );
124 return( 1 ); 124 return( 1 );
125 } 125 }
126 /* 126 /*
127 ********************************************** 127 ************************************************************
128 */ 128 */
129 static int MCGIDI_reaction_ParseReactionTypeAn 129 static int MCGIDI_reaction_ParseReactionTypeAndDetermineProducts( statusMessageReporting *smr, MCGIDI_POPs *pops, MCGIDI_reaction *reaction ) {
130 130
131 MCGIDI_outputChannel *outputChannel = &(re 131 MCGIDI_outputChannel *outputChannel = &(reaction->outputChannel);
132 int MT; 132 int MT;
133 int particlesChanges[nParticleChanges], nu 133 int particlesChanges[nParticleChanges], numberOfChanges;
134 double finalQ = 0.; 134 double finalQ = 0.;
135 135
136 if( MCGIDI_reaction_ParseDetermineReaction 136 if( MCGIDI_reaction_ParseDetermineReactionProducts( smr, pops, outputChannel, &(reaction->productsInfo), reaction, &finalQ, 0 ) != 0 ) return( 1 );
137 reaction->finalQ = finalQ; 137 reaction->finalQ = finalQ;
138 MT = MCGIDI_reaction_getENDF_MTNumber( rea 138 MT = MCGIDI_reaction_getENDF_MTNumber( reaction );
139 switch( MT ) { 139 switch( MT ) {
140 case 2 : 140 case 2 :
141 reaction->reactionType = MCGIDI_reacti 141 reaction->reactionType = MCGIDI_reactionType_elastic_e;
142 break; 142 break;
143 case 18 : case 19 : case 20 : case 21 : ca 143 case 18 : case 19 : case 20 : case 21 : case 38 :
144 reaction->reactionType = MCGIDI_reacti 144 reaction->reactionType = MCGIDI_reactionType_fission_e;
145 break; 145 break;
146 case 102 : 146 case 102 :
147 reaction->reactionType = MCGIDI_reacti 147 reaction->reactionType = MCGIDI_reactionType_capture_e;
148 break; 148 break;
149 case 5 : 149 case 5 :
150 reaction->reactionType = MCGIDI_reacti 150 reaction->reactionType = MCGIDI_reactionType_sumOfRemainingOutputChannels_e;
151 break; 151 break;
152 default : 152 default :
153 numberOfChanges = MCGIDI_reaction_part 153 numberOfChanges = MCGIDI_reaction_particleChanges( reaction->target->projectilePOP, reaction->target->targetPOP, &(reaction->productsInfo),
154 nParticleChanges, particlesChanges 154 nParticleChanges, particlesChanges );
155 155
156 reaction->reactionType = MCGIDI_reacti 156 reaction->reactionType = MCGIDI_reactionType_unknown_e;
157 if( numberOfChanges == 0 ) { 157 if( numberOfChanges == 0 ) {
158 reaction->reactionType = MCGIDI_re 158 reaction->reactionType = MCGIDI_reactionType_scattering_e; }
159 else { 159 else {
160 reaction->reactionType = MCGIDI_re 160 reaction->reactionType = MCGIDI_reactionType_nuclearIsomerTransmutation_e;
161 } 161 }
162 162
163 /* 163 /*
164 Currently, these are not handled properly: 164 Currently, these are not handled properly:
165 MCGIDI_reactionType_nuclearLevelTransi 165 MCGIDI_reactionType_nuclearLevelTransition_e
166 MCGIDI_reactionType_atomic_e 166 MCGIDI_reactionType_atomic_e
167 */ 167 */
168 break; 168 break;
169 } 169 }
170 170
171 MCGIDI_reaction_setENDL_CSNumbers( smr, re 171 MCGIDI_reaction_setENDL_CSNumbers( smr, reaction );
172 return( 0 ); 172 return( 0 );
173 } 173 }
174 /* 174 /*
175 ********************************************** 175 ************************************************************
176 */ 176 */
177 static int MCGIDI_reaction_particleChanges( MC 177 static int MCGIDI_reaction_particleChanges( MCGIDI_POP *projectile, MCGIDI_POP *target, MCGIDI_productsInfo *productsInfo, int n1, int *particlesChanges ) {
178 178
179 int projectileGlobalIndex = projectile->gl 179 int projectileGlobalIndex = projectile->globalPoPsIndex, targetGlobalIndex = target->globalPoPsIndex, i1, i2 = 0;
180 int gammaIndex = PoPs_particleIndex( "gamm 180 int gammaIndex = PoPs_particleIndex( "gamma" );
181 181
182 if( projectileGlobalIndex != gammaIndex ) 182 if( projectileGlobalIndex != gammaIndex ) {
183 for( i1 = 0; i1 < productsInfo->number 183 for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ ) if( projectileGlobalIndex == productsInfo->productInfo[i1].globalPoPsIndex ) break;
184 if( i1 == productsInfo->numberOfProduc 184 if( i1 == productsInfo->numberOfProducts ) particlesChanges[i2++] = projectileGlobalIndex;
185 } 185 }
186 186
187 for( i1 = 0; i1 < productsInfo->numberOfPr 187 for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ ) if( targetGlobalIndex == productsInfo->productInfo[i1].globalPoPsIndex ) break;
188 if( i1 == productsInfo->numberOfProducts ) 188 if( i1 == productsInfo->numberOfProducts ) particlesChanges[i2++] = targetGlobalIndex;
189 189
190 for( i1 = 0; i1 < productsInfo->numberOfPr 190 for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ ) {
191 if( i2 == n1 ) break; 191 if( i2 == n1 ) break;
192 if( /*(*/ projectileGlobalIndex == pro 192 if( /*(*/ projectileGlobalIndex == productsInfo->productInfo[i1].globalPoPsIndex /*)*/ ) continue;
193 if( /*(*/ targetGlobalIndex == product 193 if( /*(*/ targetGlobalIndex == productsInfo->productInfo[i1].globalPoPsIndex /*)*/ ) continue;
194 if( /*(*/ gammaIndex == productsInfo-> 194 if( /*(*/ gammaIndex == productsInfo->productInfo[i1].globalPoPsIndex /*)*/ ) continue;
195 particlesChanges[i2++] = productsInfo- 195 particlesChanges[i2++] = productsInfo->productInfo[i1].globalPoPsIndex;
196 } 196 }
197 return( i2 ); 197 return( i2 );
198 } 198 }
199 /* 199 /*
200 ********************************************** 200 ************************************************************
201 */ 201 */
202 static int MCGIDI_reaction_ParseDetermineReact 202 static int MCGIDI_reaction_ParseDetermineReactionProducts( statusMessageReporting *smr, MCGIDI_POPs *pops, MCGIDI_outputChannel *outputChannel,
203 MCGIDI_productsInfo *productsInfo, MCGIDI_ 203 MCGIDI_productsInfo *productsInfo, MCGIDI_reaction *reaction, double *finalQ, int level ) {
204 /* 204 /*
205 * This function determines all products that 205 * This function determines all products that can be returned during sampling for this outputChannel. Note, products like 'U238_c' and
206 * 'U238_e3' are not returned during sampling 206 * 'U238_e3' are not returned during sampling as both are decay to the groud state (unless a meta-stable is encountered).
207 * Some examples for projectile 'n' and targe 207 * Some examples for projectile 'n' and target 'U238' are:
208 * outputChannel pr 208 * outputChannel products returned during sampling.
209 * 'n + U238' n, 209 * 'n + U238' n, U238
210 * 'n + U238 + gamma' n, 210 * 'n + U238 + gamma' n, U238, gamma
211 * 'n + U238_c' n, 211 * 'n + U238_c' n, U238 (even if no gammas are give, return ground state of residual.
212 * 'n + (U238_e3 -> U238 + gamma)' n, 212 * 'n + (U238_e3 -> U238 + gamma)' n, U238, gamma
213 */ 213 */
214 int iProduct, nProducts = MCGIDI_outputCha 214 int iProduct, nProducts = MCGIDI_outputChannel_numberOfProducts( outputChannel ), globalPoPsIndex, productIsTrackable;
215 int twoBodyProductsWithData = 0; 215 int twoBodyProductsWithData = 0;
216 MCGIDI_product *product; 216 MCGIDI_product *product;
217 MCGIDI_POP *residual; 217 MCGIDI_POP *residual;
218 218
219 if( ( level == 0 ) && ( outputChannel->gen 219 if( ( level == 0 ) && ( outputChannel->genre == MCGIDI_channelGenre_twoBody_e ) ) {
220 for( iProduct = 0; iProduct < nProduct 220 for( iProduct = 0; iProduct < nProducts; iProduct++ ) {
221 product = MCGIDI_outputChannel_get 221 product = MCGIDI_outputChannel_getProductAtIndex( smr, outputChannel, iProduct );
222 if( product->pop->globalPoPsIndex 222 if( product->pop->globalPoPsIndex < 0 ) {
223 twoBodyProductsWithData = -1; 223 twoBodyProductsWithData = -1; }
224 else if( product->distribution.typ 224 else if( product->distribution.type == MCGIDI_distributionType_angular_e ) {
225 if( twoBodyProductsWithData >= 225 if( twoBodyProductsWithData >= 0 ) twoBodyProductsWithData = 1;
226 } 226 }
227 } 227 }
228 } 228 }
229 if( twoBodyProductsWithData < 0 ) twoBodyP 229 if( twoBodyProductsWithData < 0 ) twoBodyProductsWithData = 0;
230 *finalQ += MCGIDI_outputChannel_getQ_MeV( 230 *finalQ += MCGIDI_outputChannel_getQ_MeV( smr, outputChannel, 0 );
231 for( iProduct = 0; iProduct < nProducts; i 231 for( iProduct = 0; iProduct < nProducts; iProduct++ ) {
232 productIsTrackable = twoBodyProductsWi 232 productIsTrackable = twoBodyProductsWithData;
233 product = MCGIDI_outputChannel_getProd 233 product = MCGIDI_outputChannel_getProductAtIndex( smr, outputChannel, iProduct );
234 globalPoPsIndex = product->pop->global 234 globalPoPsIndex = product->pop->globalPoPsIndex;
235 if( ( product->distribution.type != MC 235 if( ( product->distribution.type != MCGIDI_distributionType_none_e ) && ( product->distribution.type != MCGIDI_distributionType_unknown_e ) ) {
236 productIsTrackable = 1; 236 productIsTrackable = 1;
237 if( globalPoPsIndex < 0 ) { 237 if( globalPoPsIndex < 0 ) {
238 if( product->distribution.angu 238 if( product->distribution.angular != NULL ) {
239 if( product->distribution. 239 if( product->distribution.angular->type == MCGIDI_angularType_recoil ) productIsTrackable = 0;
240 } 240 }
241 if( productIsTrackable ) { 241 if( productIsTrackable ) {
242 int len = (int) strlen( pr 242 int len = (int) strlen( product->pop->name );
243 243
244 if( len > 2 ) { /* Spe 244 if( len > 2 ) { /* Special case for continuum reactions with data for residual (e.g., n + U233 -> n + U233_c). */
245 if( ( product->pop->na 245 if( ( product->pop->name[len-2] == '_' ) && ( product->pop->name[len-1] == 'c' ) ) {
246 for( residual = pr 246 for( residual = product->pop; residual->globalPoPsIndex < 0; residual = residual->parent ) ;
247 productIsTrackable 247 productIsTrackable = 1;
248 globalPoPsIndex = 248 globalPoPsIndex = residual->globalPoPsIndex;
249 } 249 }
250 } 250 }
251 if( globalPoPsIndex < 0 ) 251 if( globalPoPsIndex < 0 ) {
252 smr_setReportError2( s 252 smr_setReportError2( smr, smr_unknownID, 1, "product determination for '%s' cannot be determined", product->pop->name );
253 return( 1 ); 253 return( 1 );
254 } 254 }
255 } 255 }
256 } 256 }
257 } 257 }
258 if( productIsTrackable ) { 258 if( productIsTrackable ) {
259 if( MCGIDI_reaction_addReturnProdu 259 if( MCGIDI_reaction_addReturnProduct( smr, productsInfo, globalPoPsIndex, product, reaction, 1 ) != 0 ) return( 1 ); }
260 else { 260 else {
261 if( product->decayChannel.genre != 261 if( product->decayChannel.genre != MCGIDI_channelGenre_undefined_e ) {
262 if( MCGIDI_reaction_ParseDeter 262 if( MCGIDI_reaction_ParseDetermineReactionProducts( smr, pops, &(product->decayChannel), productsInfo, reaction, finalQ, level + 1 ) != 0 ) return( 1 ); }
263 else { 263 else {
264 *finalQ += product->pop->level 264 *finalQ += product->pop->level_MeV;
265 for( residual = product->pop; 265 for( residual = product->pop; residual->globalPoPsIndex < 0; residual = residual->parent ) ;
266 if( MCGIDI_reaction_addReturnP 266 if( MCGIDI_reaction_addReturnProduct( smr, productsInfo, residual->globalPoPsIndex, product, reaction, 0 ) != 0 ) return( 1 );
267 if( product->pop->numberOfGamm 267 if( product->pop->numberOfGammaBranchs != 0 ) {
268 int gammaIndex = PoPs_part 268 int gammaIndex = PoPs_particleIndex( "gamma" );
269 if( MCGIDI_reaction_addRet 269 if( MCGIDI_reaction_addReturnProduct( smr, productsInfo, gammaIndex, NULL, reaction, 1 ) != 0 ) return( 1 );
270 } 270 }
271 } 271 }
272 } 272 }
273 } 273 }
274 return( 0 ); 274 return( 0 );
275 } 275 }
276 /* 276 /*
277 ********************************************** 277 ************************************************************
278 */ 278 */
279 static int MCGIDI_reaction_addReturnProduct( s 279 static int MCGIDI_reaction_addReturnProduct( statusMessageReporting *smr, MCGIDI_productsInfo *productsInfo, int ID, MCGIDI_product *product,
280 MCGIDI_reaction *reaction, int transpo 280 MCGIDI_reaction *reaction, int transportable ) {
281 281
282 int i1; 282 int i1;
283 enum MCGIDI_productMultiplicityType produc 283 enum MCGIDI_productMultiplicityType productMultiplicityType;
284 284
285 MCGIDI_misc_updateTransportabilitiesMap2( 285 MCGIDI_misc_updateTransportabilitiesMap2( reaction->transportabilities, ID, transportable );
286 for( i1 = 0; i1 < productsInfo->numberOfPr 286 for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ ) {
287 if( productsInfo->productInfo[i1].glob 287 if( productsInfo->productInfo[i1].globalPoPsIndex == ID ) break;
288 } 288 }
289 if( i1 == productsInfo->numberOfProducts ) 289 if( i1 == productsInfo->numberOfProducts ) {
290 if( productsInfo->numberOfProducts == 290 if( productsInfo->numberOfProducts == productsInfo->numberOfAllocatedProducts ) {
291 productsInfo->numberOfAllocatedPro 291 productsInfo->numberOfAllocatedProducts += 4;
292 if( ( productsInfo->productInfo = 292 if( ( productsInfo->productInfo = (MCGIDI_productInfo *) smr_realloc2( smr, productsInfo->productInfo,
293 productsInfo->numberOfAllocate 293 productsInfo->numberOfAllocatedProducts * sizeof( MCGIDI_productInfo ), "productsInfo->productInfo" ) ) == NULL ) return( 1 );
294 } 294 }
295 productsInfo->numberOfProducts++; 295 productsInfo->numberOfProducts++;
296 productsInfo->productInfo[i1].globalPo 296 productsInfo->productInfo[i1].globalPoPsIndex = ID;
297 productsInfo->productInfo[i1].productM 297 productsInfo->productInfo[i1].productMultiplicityType = MCGIDI_productMultiplicityType_unknown_e;
298 productsInfo->productInfo[i1].multipli 298 productsInfo->productInfo[i1].multiplicity = 0;
299 productsInfo->productInfo[i1].transpor 299 productsInfo->productInfo[i1].transportable = transportable;
300 } 300 }
301 if( product == NULL ) { 301 if( product == NULL ) {
302 productMultiplicityType = MCGIDI_produ 302 productMultiplicityType = MCGIDI_productMultiplicityType_gammaBranching_e; }
303 else { 303 else {
304 if( ( product->multiplicityVsEnergy != 304 if( ( product->multiplicityVsEnergy != NULL ) || ( product->piecewiseMultiplicities != NULL ) ) {
305 productMultiplicityType = MCGIDI_p 305 productMultiplicityType = MCGIDI_productMultiplicityType_energyDependent_e; }
306 else { 306 else {
307 productsInfo->productInfo[i1].mult 307 productsInfo->productInfo[i1].multiplicity += product->multiplicity;
308 productMultiplicityType = MCGIDI_p 308 productMultiplicityType = MCGIDI_productMultiplicityType_integer_e;
309 } 309 }
310 } 310 }
311 if( ( productsInfo->productInfo[i1].produc 311 if( ( productsInfo->productInfo[i1].productMultiplicityType == MCGIDI_productMultiplicityType_unknown_e ) ||
312 ( productsInfo->productInfo[i1].produc 312 ( productsInfo->productInfo[i1].productMultiplicityType == productMultiplicityType ) ) {
313 productsInfo->productInfo[i1].productM 313 productsInfo->productInfo[i1].productMultiplicityType = productMultiplicityType; }
314 else { 314 else {
315 productsInfo->productInfo[i1].productM 315 productsInfo->productInfo[i1].productMultiplicityType = MCGIDI_productMultiplicityType_mixed_e;
316 } 316 }
317 return( 0 ); 317 return( 0 );
318 } 318 }
319 /* 319 /*
320 ********************************************** 320 ************************************************************
321 */ 321 */
322 enum MCGIDI_reactionType MCGIDI_reaction_getRe 322 enum MCGIDI_reactionType MCGIDI_reaction_getReactionType( statusMessageReporting * /*smr*/, MCGIDI_reaction *reaction ) {
323 323
324 return( reaction->reactionType ); 324 return( reaction->reactionType );
325 } 325 }
326 /* 326 /*
327 ********************************************** 327 ************************************************************
328 */ 328 */
329 MCGIDI_target_heated *MCGIDI_reaction_getTarge 329 MCGIDI_target_heated *MCGIDI_reaction_getTargetHeated( statusMessageReporting * /*smr*/, MCGIDI_reaction *reaction ) {
330 330
331 return( reaction->target ); 331 return( reaction->target );
332 } 332 }
333 /* 333 /*
334 ********************************************** 334 ************************************************************
335 */ 335 */
336 double MCGIDI_reaction_getProjectileMass_MeV( 336 double MCGIDI_reaction_getProjectileMass_MeV( statusMessageReporting *smr, MCGIDI_reaction *reaction ) {
337 337
338 return( MCGIDI_target_heated_getProjectile 338 return( MCGIDI_target_heated_getProjectileMass_MeV( smr, reaction->target ) );
339 } 339 }
340 /* 340 /*
341 ********************************************** 341 ************************************************************
342 */ 342 */
343 double MCGIDI_reaction_getTargetMass_MeV( stat 343 double MCGIDI_reaction_getTargetMass_MeV( statusMessageReporting *smr, MCGIDI_reaction *reaction ) {
344 344
345 return( MCGIDI_target_heated_getTargetMass 345 return( MCGIDI_target_heated_getTargetMass_MeV( smr, reaction->target ) );
346 } 346 }
347 /* 347 /*
348 ********************************************** 348 ************************************************************
349 */ 349 */
350 int MCGIDI_reaction_getDomain( statusMessageRe 350 int MCGIDI_reaction_getDomain( statusMessageReporting * /*smr*/, MCGIDI_reaction *reaction, double *EMin, double *EMax ) {
351 /* 351 /*
352 * Return value 352 * Return value
353 * < 0 No cross section data. 353 * < 0 No cross section data.
354 * == 0 Okay and EMin and EMax set. 354 * == 0 Okay and EMin and EMax set.
355 * > 0 error, EMin and EMax undefined 355 * > 0 error, EMin and EMax undefined.
356 */ 356 */
357 357
358 if( !reaction->domainValuesPresent ) retur 358 if( !reaction->domainValuesPresent ) return( -1 );
359 *EMin = reaction->EMin; 359 *EMin = reaction->EMin;
360 *EMax = reaction->EMax; 360 *EMax = reaction->EMax;
361 return( 0 ); 361 return( 0 );
362 } 362 }
363 /* 363 /*
364 ********************************************** 364 ************************************************************
365 */ 365 */
366 int MCGIDI_reaction_fixDomains( statusMessageR 366 int MCGIDI_reaction_fixDomains( statusMessageReporting * /*smr*/, MCGIDI_reaction *reaction, double EMin, double EMax, nfu_status *status ) {
367 367
368 double lowerEps = 1e-14, upperEps = -1e-14 368 double lowerEps = 1e-14, upperEps = -1e-14;
369 369
370 if( reaction->EMin == EMin ) lowerEps = 0. 370 if( reaction->EMin == EMin ) lowerEps = 0.;
371 if( reaction->EMax == EMax ) upperEps = 0. 371 if( reaction->EMax == EMax ) upperEps = 0.;
372 if( ( lowerEps == 0. ) && ( upperEps == 0. 372 if( ( lowerEps == 0. ) && ( upperEps == 0. ) ) return( 0 );
373 373
374 *status = ptwXY_dullEdges( reaction->cross 374 *status = ptwXY_dullEdges( reaction->crossSection, lowerEps, upperEps, 1 );
375 return( *status != nfu_Okay ); 375 return( *status != nfu_Okay );
376 } 376 }
377 /* 377 /*
378 ********************************************** 378 ************************************************************
379 */ 379 */
380 double MCGIDI_reaction_getCrossSectionAtE( sta 380 double MCGIDI_reaction_getCrossSectionAtE( statusMessageReporting *smr, MCGIDI_reaction *reaction, MCGIDI_quantitiesLookupModes &modes,
381 bool sampling ) { 381 bool sampling ) {
382 382
383 double e_in = modes.getProjectileEnergy( ) 383 double e_in = modes.getProjectileEnergy( ), xsec;
384 384
385 if( modes.getCrossSectionMode( ) == MCGIDI 385 if( modes.getCrossSectionMode( ) == MCGIDI_quantityLookupMode_pointwise ) {
386 if( e_in < reaction->EMin ) e_in = rea 386 if( e_in < reaction->EMin ) e_in = reaction->EMin;
387 if( e_in > reaction->EMax ) e_in = rea 387 if( e_in > reaction->EMax ) e_in = reaction->EMax;
388 ptwXY_getValueAtX( reaction->crossSect 388 ptwXY_getValueAtX( reaction->crossSection, e_in, &xsec ); }
389 else if( modes.getCrossSectionMode( ) == M 389 else if( modes.getCrossSectionMode( ) == MCGIDI_quantityLookupMode_grouped ) {
390 int index = modes.getGroupIndex( ); 390 int index = modes.getGroupIndex( );
391 double *xSecP = ptwX_getPointAtIndex( 391 double *xSecP = ptwX_getPointAtIndex( reaction->crossSectionGrouped, index );
392 392
393 if( xSecP != NULL ) { 393 if( xSecP != NULL ) {
394 xsec = *xSecP; 394 xsec = *xSecP;
395 if( sampling && ( index == reactio 395 if( sampling && ( index == reaction->thresholdGroupIndex ) ) xsec += reaction->thresholdGroupedDeltaCrossSection; }
396 else { 396 else {
397 xsec = 0.; 397 xsec = 0.;
398 smr_setReportError2( smr, smr_unkn 398 smr_setReportError2( smr, smr_unknownID, 1, "Invalid cross section group index %d", index );
399 } } 399 } }
400 else { 400 else {
401 xsec = 0.; 401 xsec = 0.;
402 } 402 }
403 return( xsec ); 403 return( xsec );
404 } 404 }
405 /* 405 /*
406 ********************************************** 406 ************************************************************
407 */ 407 */
408 double MCGIDI_reaction_getFinalQ( statusMessag 408 double MCGIDI_reaction_getFinalQ( statusMessageReporting * /*smr*/, MCGIDI_reaction *reaction, MCGIDI_quantitiesLookupModes &/*modes*/ ) {
409 409
410 return( reaction->finalQ ); 410 return( reaction->finalQ );
411 } 411 }
412 /* 412 /*
413 ********************************************** 413 ************************************************************
414 */ 414 */
415 int MCGIDI_reaction_getENDF_MTNumber( MCGIDI_r 415 int MCGIDI_reaction_getENDF_MTNumber( MCGIDI_reaction *reaction ) {
416 416
417 return( reaction->ENDF_MT ); 417 return( reaction->ENDF_MT );
418 } 418 }
419 /* 419 /*
420 ********************************************** 420 ************************************************************
421 */ 421 */
422 int MCGIDI_reaction_getENDL_CSNumbers( MCGIDI_ 422 int MCGIDI_reaction_getENDL_CSNumbers( MCGIDI_reaction *reaction, int *S ) {
423 423
424 if( S != NULL ) *S = reaction->ENDL_S; 424 if( S != NULL ) *S = reaction->ENDL_S;
425 return( reaction->ENDL_C ); 425 return( reaction->ENDL_C );
426 } 426 }
427 /* 427 /*
428 ********************************************** 428 ************************************************************
429 */ 429 */
430 static int MCGIDI_reaction_setENDL_CSNumbers( 430 static int MCGIDI_reaction_setENDL_CSNumbers( statusMessageReporting * /*smr*/, MCGIDI_reaction *reaction ) {
431 431
432 int MT = MCGIDI_reaction_getENDF_MTNumber( 432 int MT = MCGIDI_reaction_getENDF_MTNumber( reaction );
433 int MT1_50ToC[] = { 1, 10, -3, -4, 433 int MT1_50ToC[] = { 1, 10, -3, -4, -5, 0, 0, 0, 0, -10,
434 32, 0, 0, 0, 434 32, 0, 0, 0, 0, 12, 13, 15, 15, 15,
435 15, 26, 36, 33, - 435 15, 26, 36, 33, -25, 0, -27, 20, 27, -30,
436 0, 22, 24, 25, - 436 0, 22, 24, 25, -35, -36, 14, 15, 0, 0,
437 29, 16, 0, 17, 437 29, 16, 0, 17, 34, 0, 0, 0, 0 };
438 int MT100_200ToC[] = { -101, 46, 40, 438 int MT100_200ToC[] = { -101, 46, 40, 41, 42, 44, 45, 37, -109, 0,
439 18, 48, -113, - 439 18, 48, -113, -114, 19, 39, 47, 0, 0, 0,
440 0, 0, 0, 440 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
441 0, 0, 0, 441 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
442 0, 0, 0, 442 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
443 0, -152, -153, - 443 0, -152, -153, -154, 43, -156, -157, 23, 31, -160,
444 -161, -162, -163, - 444 -161, -162, -163, -164, -165, -166, -167, -168, -169, -170,
445 -171, -172, -173, - 445 -171, -172, -173, -174, -175, -176, -177, -178, -179, -180,
446 -181, -182, -183, - 446 -181, -182, -183, -184, -185, -186, -187, -188, 28, -190,
447 -191, -192, 38, - 447 -191, -192, 38, -194, -195, -196, -197, -198, -199, -200 };
448 448
449 reaction->ENDL_C = 0; 449 reaction->ENDL_C = 0;
450 reaction->ENDL_S = 0; 450 reaction->ENDL_S = 0;
451 if( MT <= 0 ) return( 1 ); 451 if( MT <= 0 ) return( 1 );
452 if( MT > 891 ) return( 1 ); 452 if( MT > 891 ) return( 1 );
453 if( MT < 50 ) { 453 if( MT < 50 ) {
454 reaction->ENDL_C = MT1_50ToC[MT - 1]; 454 reaction->ENDL_C = MT1_50ToC[MT - 1]; }
455 else if( MT <= 91 ) { 455 else if( MT <= 91 ) {
456 reaction->ENDL_C = 11; 456 reaction->ENDL_C = 11;
457 if( MT != 91 ) reaction->ENDL_S = 1; } 457 if( MT != 91 ) reaction->ENDL_S = 1; }
458 else if( ( MT > 100 ) && ( MT <= 200 ) ) { 458 else if( ( MT > 100 ) && ( MT <= 200 ) ) {
459 reaction->ENDL_C = MT100_200ToC[MT - 1 459 reaction->ENDL_C = MT100_200ToC[MT - 101]; }
460 else if( ( MT == 452 ) || ( MT == 455 ) || 460 else if( ( MT == 452 ) || ( MT == 455 ) || ( MT == 456 ) || ( MT == 458 ) ) {
461 reaction->ENDL_C = 15; 461 reaction->ENDL_C = 15;
462 if( MT == 455 ) reaction->ENDL_S = 7; 462 if( MT == 455 ) reaction->ENDL_S = 7; }
463 else if( MT >= 600 ) { 463 else if( MT >= 600 ) {
464 if( MT < 650 ) { 464 if( MT < 650 ) {
465 reaction->ENDL_C = 40; 465 reaction->ENDL_C = 40;
466 if( MT != 649 ) reaction->ENDL_S = 466 if( MT != 649 ) reaction->ENDL_S = 1; }
467 else if( MT < 700 ) { 467 else if( MT < 700 ) {
468 reaction->ENDL_C = 41; 468 reaction->ENDL_C = 41;
469 if( MT != 699 ) reaction->ENDL_S = 469 if( MT != 699 ) reaction->ENDL_S = 1; }
470 else if( MT < 750 ) { 470 else if( MT < 750 ) {
471 reaction->ENDL_C = 42; 471 reaction->ENDL_C = 42;
472 if( MT != 749 ) reaction->ENDL_S = 472 if( MT != 749 ) reaction->ENDL_S = 1; }
473 else if( MT < 800 ) { 473 else if( MT < 800 ) {
474 reaction->ENDL_C = 44; 474 reaction->ENDL_C = 44;
475 if( MT != 799 ) reaction->ENDL_S = 475 if( MT != 799 ) reaction->ENDL_S = 1; }
476 else if( MT < 850 ) { 476 else if( MT < 850 ) {
477 reaction->ENDL_C = 45; 477 reaction->ENDL_C = 45;
478 if( MT != 849 ) reaction->ENDL_S = 478 if( MT != 849 ) reaction->ENDL_S = 1; }
479 else if( ( MT >= 875 ) && ( MT <= 891 479 else if( ( MT >= 875 ) && ( MT <= 891 ) ) {
480 reaction->ENDL_C = 12; 480 reaction->ENDL_C = 12;
481 if( MT != 891 ) reaction->ENDL_S = 481 if( MT != 891 ) reaction->ENDL_S = 1;
482 } 482 }
483 } 483 }
484 return( 0 ); 484 return( 0 );
485 } 485 }
486 /* 486 /*
487 ********************************************** 487 ************************************************************
488 */ 488 */
489 MCGIDI_productsInfo *MCGIDI_reaction_getProduc 489 MCGIDI_productsInfo *MCGIDI_reaction_getProductsInfo( MCGIDI_reaction *reaction ) {
490 490
491 return( &(reaction->productsInfo) ); 491 return( &(reaction->productsInfo) );
492 } 492 }
493 /* 493 /*
494 ********************************************** 494 ************************************************************
495 */ 495 */
496 int MCGIDI_reaction_recast( statusMessageRepor 496 int MCGIDI_reaction_recast( statusMessageReporting *smr, MCGIDI_reaction *reaction, GIDI_settings & /*settings*/,
497 GIDI_settings_particle const *projecti 497 GIDI_settings_particle const *projectileSettings, double temperature_MeV, ptwXPoints *totalGroupedCrossSection ) {
498 498
499 if( totalGroupedCrossSection != NULL ) { 499 if( totalGroupedCrossSection != NULL ) {
500 nfu_status status_nf; 500 nfu_status status_nf;
501 GIDI_settings_group group( projectileS 501 GIDI_settings_group group( projectileSettings->getGroup( ) );
502 502
503 if( reaction->crossSectionGrouped != N 503 if( reaction->crossSectionGrouped != NULL ) reaction->crossSectionGrouped = ptwX_free( reaction->crossSectionGrouped );
504 if( ( reaction->crossSectionGrouped = 504 if( ( reaction->crossSectionGrouped = projectileSettings->groupFunction( smr, reaction->crossSection, temperature_MeV, 0 ) ) == NULL ) return( 1 );
505 if( ( status_nf = ptwX_add_ptwX( total 505 if( ( status_nf = ptwX_add_ptwX( totalGroupedCrossSection, reaction->crossSectionGrouped ) ) != nfu_Okay ) return( 1 );
506 506
507 reaction->thresholdGroupDomain = react 507 reaction->thresholdGroupDomain = reaction->thresholdGroupedDeltaCrossSection = 0.;
508 reaction->thresholdGroupIndex = group. 508 reaction->thresholdGroupIndex = group.getGroupIndexFromEnergy( reaction->EMin, false );
509 if( reaction->thresholdGroupIndex > -1 509 if( reaction->thresholdGroupIndex > -1 ) {
510 reaction->thresholdGroupDomain = g 510 reaction->thresholdGroupDomain = group[reaction->thresholdGroupIndex+1] - reaction->EMin;
511 if( reaction->thresholdGroupDomain 511 if( reaction->thresholdGroupDomain > 0 ) {
512 512 /* factor 2 for linear reject in bin but above threshold. */
513 reaction->thresholdGroupedDelt 513 reaction->thresholdGroupedDeltaCrossSection = *ptwX_getPointAtIndex( reaction->crossSectionGrouped, reaction->thresholdGroupIndex ) *
514 ( 2 * ( group[reaction 514 ( 2 * ( group[reaction->thresholdGroupIndex+1] - group[reaction->thresholdGroupIndex] ) / reaction->thresholdGroupDomain - 1 );
515 } 515 }
516 } 516 }
517 } 517 }
518 return( 0 ); 518 return( 0 );
519 } 519 }
520 520
521 /* 521 /*
522 *********************** productsInfo ********* 522 *********************** productsInfo ***********************
523 */ 523 */
524 /* 524 /*
525 ********************************************** 525 ************************************************************
526 */ 526 */
527 int MCGIDI_productsInfo_getNumberOfUniqueProdu 527 int MCGIDI_productsInfo_getNumberOfUniqueProducts( MCGIDI_productsInfo *productsInfo ) {
528 528
529 return( productsInfo->numberOfProducts ); 529 return( productsInfo->numberOfProducts );
530 } 530 }
531 /* 531 /*
532 ********************************************** 532 ************************************************************
533 */ 533 */
534 int MCGIDI_productsInfo_getPoPsIndexAtIndex( M 534 int MCGIDI_productsInfo_getPoPsIndexAtIndex( MCGIDI_productsInfo *productsInfo, int index ) {
535 535
536 if( ( index < 0 ) || ( index >= productsIn 536 if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) ) return( -1 );
537 return( productsInfo->productInfo[index].g 537 return( productsInfo->productInfo[index].globalPoPsIndex );
538 } 538 }
539 /* 539 /*
540 ********************************************** 540 ************************************************************
541 */ 541 */
542 enum MCGIDI_productMultiplicityType MCGIDI_pro 542 enum MCGIDI_productMultiplicityType MCGIDI_productsInfo_getMultiplicityTypeAtIndex( MCGIDI_productsInfo *productsInfo, int index ) {
543 543
544 if( ( index < 0 ) || ( index >= productsIn 544 if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) ) return( MCGIDI_productMultiplicityType_invalid_e );
545 return( productsInfo->productInfo[index].p 545 return( productsInfo->productInfo[index].productMultiplicityType );
546 } 546 }
547 /* 547 /*
548 ********************************************** 548 ************************************************************
549 */ 549 */
550 int MCGIDI_productsInfo_getIntegerMultiplicity 550 int MCGIDI_productsInfo_getIntegerMultiplicityAtIndex( MCGIDI_productsInfo *productsInfo, int index ) {
551 551
552 if( ( index < 0 ) || ( index >= productsIn 552 if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) ) return( -1 );
553 return( productsInfo->productInfo[index].m 553 return( productsInfo->productInfo[index].multiplicity );
554 } 554 }
555 /* 555 /*
556 ********************************************** 556 ************************************************************
557 */ 557 */
558 int MCGIDI_productsInfo_getTransportableAtInde 558 int MCGIDI_productsInfo_getTransportableAtIndex( MCGIDI_productsInfo *productsInfo, int index ) {
559 559
560 if( ( index < 0 ) || ( index >= productsIn 560 if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) ) return( -1 );
561 return( productsInfo->productInfo[index].t 561 return( productsInfo->productInfo[index].transportable );
562 } 562 }
563 563
564 #if defined __cplusplus 564 #if defined __cplusplus
565 } 565 }
566 #endif 566 #endif
567 567
568 568