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