10 #include "MCGIDI_misc.h"
11 #include "MCGIDI_private.h"
13 #if defined __cplusplus
18 #define nParticleChanges 6
21 static int MCGIDI_reaction_particleChanges( MCGIDI_POP *projectile, MCGIDI_POP *target, MCGIDI_productsInfo *productsInfo,
int n1,
int *particlesChanges );
24 MCGIDI_productsInfo *productsInfo, MCGIDI_reaction *reaction,
double *finalQ,
int level );
26 MCGIDI_reaction *reaction,
int transportable );
33 MCGIDI_reaction *reaction;
35 if( ( reaction = (MCGIDI_reaction *) smr_malloc2( smr,
sizeof( MCGIDI_reaction ), 0,
"reaction" ) ) == NULL )
return( NULL );
45 reaction->transportabilities =
new transportabilitiesMap( );
53 memset( reaction, 0,
sizeof( MCGIDI_reaction ) );
72 ptwX_free( reaction->crossSectionGrouped );
76 if( reaction->productsInfo.productInfo != NULL )
smr_freeMemory( (
void **) &(reaction->productsInfo.productInfo) );
77 delete reaction->transportabilities;
85 MCGIDI_POPs *pops, MCGIDI_reaction *reaction ) {
87 xDataTOM_element *child, *linear, *outputChannel;
88 enum xDataTOM_interpolationFlag independent, dependent;
89 enum xDataTOM_interpolationQualifier qualifier;
90 char const *outputChannelStr, *crossSectionUnits[2] = {
"MeV",
"b" };
94 reaction->target = target;
95 reaction->reactionType = MCGIDI_reactionType_unknown_e;
99 if( ( reaction->outputChannelStr = smr_allocateCopyString2( smr, outputChannelStr,
"reaction->outputChannelStr" ) ) == NULL )
goto err;
106 if( ( independent != xDataTOM_interpolationFlag_linear ) || ( dependent != xDataTOM_interpolationFlag_linear ) ) {
107 smr_setReportError2( smr, smr_unknownID, 1,
"cross section interpolation (%d,%d) is not linear-linear", independent, dependent );
111 reaction->domainValuesPresent = 1;
131 MCGIDI_outputChannel *outputChannel = &(reaction->outputChannel);
137 reaction->finalQ = finalQ;
141 reaction->reactionType = MCGIDI_reactionType_elastic_e;
143 case 18 :
case 19 :
case 20 :
case 21 :
case 38 :
144 reaction->reactionType = MCGIDI_reactionType_fission_e;
147 reaction->reactionType = MCGIDI_reactionType_capture_e;
150 reaction->reactionType = MCGIDI_reactionType_sumOfRemainingOutputChannels_e;
156 reaction->reactionType = MCGIDI_reactionType_unknown_e;
157 if( numberOfChanges == 0 ) {
158 reaction->reactionType = MCGIDI_reactionType_scattering_e; }
160 reaction->reactionType = MCGIDI_reactionType_nuclearIsomerTransmutation_e;
179 int projectileGlobalIndex = projectile->globalPoPsIndex, targetGlobalIndex = target->globalPoPsIndex, i1, i2 = 0;
182 if( projectileGlobalIndex != gammaIndex ) {
183 for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ )
if( projectileGlobalIndex == productsInfo->productInfo[i1].globalPoPsIndex )
break;
184 if( i1 == productsInfo->numberOfProducts ) particlesChanges[i2++] = projectileGlobalIndex;
187 for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ )
if( targetGlobalIndex == productsInfo->productInfo[i1].globalPoPsIndex )
break;
188 if( i1 == productsInfo->numberOfProducts ) particlesChanges[i2++] = targetGlobalIndex;
190 for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ ) {
191 if( i2 == n1 )
break;
192 if( projectileGlobalIndex == productsInfo->productInfo[i1].globalPoPsIndex )
continue;
193 if( targetGlobalIndex == productsInfo->productInfo[i1].globalPoPsIndex )
continue;
194 if( gammaIndex == productsInfo->productInfo[i1].globalPoPsIndex )
continue;
195 particlesChanges[i2++] = productsInfo->productInfo[i1].globalPoPsIndex;
203 MCGIDI_productsInfo *productsInfo, MCGIDI_reaction *reaction,
double *finalQ,
int level ) {
215 int twoBodyProductsWithData = 0;
216 MCGIDI_product *product;
217 MCGIDI_POP *residual;
219 if( ( level == 0 ) && ( outputChannel->genre == MCGIDI_channelGenre_twoBody_e ) ) {
220 for( iProduct = 0; iProduct < nProducts; iProduct++ ) {
222 if( product->pop->globalPoPsIndex < 0 ) {
223 twoBodyProductsWithData = -1; }
224 else if( product->distribution.type == MCGIDI_distributionType_angular_e ) {
225 if( twoBodyProductsWithData >= 0 ) twoBodyProductsWithData = 1;
229 if( twoBodyProductsWithData < 0 ) twoBodyProductsWithData = 0;
231 for( iProduct = 0; iProduct < nProducts; iProduct++ ) {
232 productIsTrackable = twoBodyProductsWithData;
234 globalPoPsIndex = product->pop->globalPoPsIndex;
235 if( ( product->distribution.type != MCGIDI_distributionType_none_e ) && ( product->distribution.type != MCGIDI_distributionType_unknown_e ) ) {
236 productIsTrackable = 1;
237 if( globalPoPsIndex < 0 ) {
238 if( product->distribution.angular != NULL ) {
239 if( product->distribution.angular->type == MCGIDI_angularType_recoil ) productIsTrackable = 0;
241 if( productIsTrackable ) {
242 int len = (int) strlen( product->pop->name );
245 if( ( product->pop->name[len-2] ==
'_' ) && ( product->pop->name[len-1] ==
'c' ) ) {
246 for( residual = product->pop; residual->globalPoPsIndex < 0; residual = residual->parent ) ;
247 productIsTrackable = 1;
248 globalPoPsIndex = residual->globalPoPsIndex;
251 if( globalPoPsIndex < 0 ) {
252 smr_setReportError2( smr, smr_unknownID, 1,
"product determination for '%s' cannot be determined", product->pop->name );
258 if( productIsTrackable ) {
261 if( product->decayChannel.genre != MCGIDI_channelGenre_undefined_e ) {
264 *finalQ += product->pop->level_MeV;
265 for( residual = product->pop; residual->globalPoPsIndex < 0; residual = residual->parent ) ;
267 if( product->pop->numberOfGammaBranchs != 0 ) {
280 MCGIDI_reaction *reaction,
int transportable ) {
283 enum MCGIDI_productMultiplicityType productMultiplicityType;
286 for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ ) {
287 if( productsInfo->productInfo[i1].globalPoPsIndex == ID )
break;
289 if( i1 == productsInfo->numberOfProducts ) {
290 if( productsInfo->numberOfProducts == productsInfo->numberOfAllocatedProducts ) {
291 productsInfo->numberOfAllocatedProducts += 4;
292 if( ( productsInfo->productInfo = (MCGIDI_productInfo *) smr_realloc2( smr, productsInfo->productInfo,
293 productsInfo->numberOfAllocatedProducts *
sizeof( MCGIDI_productInfo ),
"productsInfo->productInfo" ) ) == NULL )
return( 1 );
295 productsInfo->numberOfProducts++;
296 productsInfo->productInfo[i1].globalPoPsIndex = ID;
297 productsInfo->productInfo[i1].productMultiplicityType = MCGIDI_productMultiplicityType_unknown_e;
298 productsInfo->productInfo[i1].multiplicity = 0;
299 productsInfo->productInfo[i1].transportable = transportable;
301 if( product == NULL ) {
302 productMultiplicityType = MCGIDI_productMultiplicityType_gammaBranching_e; }
304 if( ( product->multiplicityVsEnergy != NULL ) || ( product->piecewiseMultiplicities != NULL ) ) {
305 productMultiplicityType = MCGIDI_productMultiplicityType_energyDependent_e; }
307 productsInfo->productInfo[i1].multiplicity += product->multiplicity;
308 productMultiplicityType = MCGIDI_productMultiplicityType_integer_e;
311 if( ( productsInfo->productInfo[i1].productMultiplicityType == MCGIDI_productMultiplicityType_unknown_e ) ||
312 ( productsInfo->productInfo[i1].productMultiplicityType == productMultiplicityType ) ) {
313 productsInfo->productInfo[i1].productMultiplicityType = productMultiplicityType; }
315 productsInfo->productInfo[i1].productMultiplicityType = MCGIDI_productMultiplicityType_mixed_e;
324 return( reaction->reactionType );
331 return( reaction->target );
358 if( !reaction->domainValuesPresent )
return( -1 );
359 *EMin = reaction->EMin;
360 *EMax = reaction->EMax;
368 double lowerEps = 1e-14, upperEps = -1e-14;
370 if( reaction->EMin == EMin ) lowerEps = 0.;
371 if( reaction->EMax == EMax ) upperEps = 0.;
372 if( ( lowerEps == 0. ) && ( upperEps == 0. ) )
return( 0 );
374 *status =
ptwXY_dullEdges( reaction->crossSection, lowerEps, upperEps, 1 );
375 return( *status != nfu_Okay );
383 double e_in = modes.getProjectileEnergy( ), xsec;
385 if( modes.getCrossSectionMode( ) == MCGIDI_quantityLookupMode_pointwise ) {
386 if( e_in < reaction->EMin ) e_in = reaction->EMin;
387 if( e_in > reaction->EMax ) e_in = reaction->EMax;
389 else if( modes.getCrossSectionMode( ) == MCGIDI_quantityLookupMode_grouped ) {
390 int index = modes.getGroupIndex( );
393 if( xSecP != NULL ) {
395 if( sampling && ( index == reaction->thresholdGroupIndex ) ) xsec += reaction->thresholdGroupedDeltaCrossSection; }
398 smr_setReportError2( smr, smr_unknownID, 1,
"Invalid cross section group index %d", index );
410 return( reaction->finalQ );
417 return( reaction->ENDF_MT );
424 if( S != NULL ) *S = reaction->ENDL_S;
425 return( reaction->ENDL_C );
433 int MT1_50ToC[] = { 1, 10, -3, -4, -5, 0, 0, 0, 0, -10,
434 32, 0, 0, 0, 0, 12, 13, 15, 15, 15,
435 15, 26, 36, 33, -25, 0, -27, 20, 27, -30,
436 0, 22, 24, 25, -35, -36, 14, 15, 0, 0,
437 29, 16, 0, 17, 34, 0, 0, 0, 0 };
438 int MT100_200ToC[] = { -101, 46, 40, 41, 42, 44, 45, 37, -109, 0,
439 18, 48, -113, -114, 19, 39, 47, 0, 0, 0,
440 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
441 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
442 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
443 0, -152, -153, -154, 43, -156, -157, 23, 31, -160,
444 -161, -162, -163, -164, -165, -166, -167, -168, -169, -170,
445 -171, -172, -173, -174, -175, -176, -177, -178, -179, -180,
446 -181, -182, -183, -184, -185, -186, -187, -188, 28, -190,
447 -191, -192, 38, -194, -195, -196, -197, -198, -199, -200 };
449 reaction->ENDL_C = 0;
450 reaction->ENDL_S = 0;
451 if( MT <= 0 )
return( 1 );
452 if( MT > 891 )
return( 1 );
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 - 101]; }
460 else if( ( MT == 452 ) || ( MT == 455 ) || ( MT == 456 ) || ( MT == 458 ) ) {
461 reaction->ENDL_C = 15;
462 if( MT == 455 ) reaction->ENDL_S = 7; }
463 else if( MT >= 600 ) {
465 reaction->ENDL_C = 40;
466 if( MT != 649 ) reaction->ENDL_S = 1; }
467 else if( MT < 700 ) {
468 reaction->ENDL_C = 41;
469 if( MT != 699 ) reaction->ENDL_S = 1; }
470 else if( MT < 750 ) {
471 reaction->ENDL_C = 42;
472 if( MT != 749 ) reaction->ENDL_S = 1; }
473 else if( MT < 800 ) {
474 reaction->ENDL_C = 44;
475 if( MT != 799 ) reaction->ENDL_S = 1; }
476 else if( MT < 850 ) {
477 reaction->ENDL_C = 45;
478 if( MT != 849 ) reaction->ENDL_S = 1; }
479 else if( ( MT >= 875 ) && ( MT <= 891 ) ) {
480 reaction->ENDL_C = 12;
481 if( MT != 891 ) reaction->ENDL_S = 1;
491 return( &(reaction->productsInfo) );
497 GIDI_settings_particle const *projectileSettings,
double temperature_MeV, ptwXPoints *totalGroupedCrossSection ) {
499 if( totalGroupedCrossSection != NULL ) {
500 nfu_status status_nf;
503 if( reaction->crossSectionGrouped != NULL ) reaction->crossSectionGrouped =
ptwX_free( reaction->crossSectionGrouped );
504 if( ( reaction->crossSectionGrouped = projectileSettings->
groupFunction( smr, reaction->crossSection, temperature_MeV, 0 ) ) == NULL )
return( 1 );
505 if( ( status_nf =
ptwX_add_ptwX( totalGroupedCrossSection, reaction->crossSectionGrouped ) ) != nfu_Okay )
return( 1 );
507 reaction->thresholdGroupDomain = reaction->thresholdGroupedDeltaCrossSection = 0.;
508 reaction->thresholdGroupIndex = group.getGroupIndexFromEnergy( reaction->EMin,
false );
509 if( reaction->thresholdGroupIndex > -1 ) {
510 reaction->thresholdGroupDomain = group[reaction->thresholdGroupIndex+1] - reaction->EMin;
511 if( reaction->thresholdGroupDomain > 0 ) {
513 reaction->thresholdGroupedDeltaCrossSection = *
ptwX_getPointAtIndex( reaction->crossSectionGrouped, reaction->thresholdGroupIndex ) *
514 ( 2 * ( group[reaction->thresholdGroupIndex+1] - group[reaction->thresholdGroupIndex] ) / reaction->thresholdGroupDomain - 1 );
529 return( productsInfo->numberOfProducts );
536 if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) )
return( -1 );
537 return( productsInfo->productInfo[index].globalPoPsIndex );
544 if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) )
return( MCGIDI_productMultiplicityType_invalid_e );
545 return( productsInfo->productInfo[index].productMultiplicityType );
552 if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) )
return( -1 );
553 return( productsInfo->productInfo[index].multiplicity );
560 if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) )
return( -1 );
561 return( productsInfo->productInfo[index].transportable );
564 #if defined __cplusplus
xDataTOM_element * xDataTOME_getOneElementByName(statusMessageReporting *smr, xDataTOM_element *element, char const *name, int required)
int MCGIDI_productsInfo_getTransportableAtIndex(MCGIDI_productsInfo *productsInfo, int index)
int MCGIDI_reaction_fixDomains(statusMessageReporting *, MCGIDI_reaction *reaction, double EMin, double EMax, nfu_status *status)
int MCGIDI_productsInfo_getIntegerMultiplicityAtIndex(MCGIDI_productsInfo *productsInfo, int index)
GIDI::ptwXPoints * groupFunction(GIDI::statusMessageReporting *smr, GIDI::ptwXYPoints *ptwXY1, double temperature, int order) const
int MCGIDI_outputChannel_release(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel)
nfu_status ptwX_add_ptwX(ptwXPoints *ptwX1, ptwXPoints *ptwX2)
int MCGIDI_productsInfo_getPoPsIndexAtIndex(MCGIDI_productsInfo *productsInfo, int index)
double MCGIDI_reaction_getCrossSectionAtE(statusMessageReporting *smr, MCGIDI_reaction *reaction, MCGIDI_quantitiesLookupModes &modes, bool sampling)
int MCGIDI_reaction_getDomain(statusMessageReporting *, MCGIDI_reaction *reaction, double *EMin, double *EMax)
int PoPs_particleIndex(char const *name)
double MCGIDI_outputChannel_getQ_MeV(statusMessageReporting *, MCGIDI_outputChannel *outputChannel, double)
int MCGIDI_reaction_release(statusMessageReporting *smr, MCGIDI_reaction *reaction)
enum MCGIDI_reactionType MCGIDI_reaction_getReactionType(statusMessageReporting *, MCGIDI_reaction *reaction)
double MCGIDI_reaction_getTargetMass_MeV(statusMessageReporting *smr, MCGIDI_reaction *reaction)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
GIDI_settings_group getGroup(void) const
void xDataTOMAL_initial(statusMessageReporting *, xDataTOM_attributionList *attributes)
double ptwXY_getXMax(ptwXYPoints *ptwXY)
MCGIDI_reaction * MCGIDI_reaction_new(statusMessageReporting *smr)
int MCGIDI_reaction_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_target_heated *target, MCGIDI_POPs *pops, MCGIDI_reaction *reaction)
int MCGIDI_reaction_initialize(statusMessageReporting *smr, MCGIDI_reaction *reaction)
static int MCGIDI_reaction_particleChanges(MCGIDI_POP *projectile, MCGIDI_POP *target, MCGIDI_productsInfo *productsInfo, int n1, int *particlesChanges)
double MCGIDI_target_heated_getProjectileMass_MeV(statusMessageReporting *, MCGIDI_target_heated *target)
int MCGIDI_reaction_recast(statusMessageReporting *smr, MCGIDI_reaction *reaction, GIDI_settings &, GIDI_settings_particle const *projectileSettings, double temperature_MeV, ptwXPoints *totalGroupedCrossSection)
int xDataTOME_convertAttributeToInteger(statusMessageReporting *smr, xDataTOM_element *element, char const *name, int *n)
MCGIDI_reaction * MCGIDI_reaction_free(statusMessageReporting *smr, MCGIDI_reaction *reaction)
void MCGIDI_misc_updateTransportabilitiesMap2(transportabilitiesMap *transportabilities, int PoPID, int transportable)
int MCGIDI_outputChannel_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_POPs *pops, MCGIDI_outputChannel *outputChannel, MCGIDI_reaction *reaction, MCGIDI_product *parent)
MCGIDI_target_heated * MCGIDI_reaction_getTargetHeated(statusMessageReporting *, MCGIDI_reaction *reaction)
void * smr_freeMemory(void **p)
MCGIDI_product * MCGIDI_outputChannel_getProductAtIndex(statusMessageReporting *smr, MCGIDI_outputChannel *outputChannel, int i)
enum MCGIDI_productMultiplicityType MCGIDI_productsInfo_getMultiplicityTypeAtIndex(MCGIDI_productsInfo *productsInfo, int index)
int xDataTOME_getInterpolation(statusMessageReporting *smr, xDataTOM_element *element, int index, enum xDataTOM_interpolationFlag *independent, enum xDataTOM_interpolationFlag *dependent, enum xDataTOM_interpolationQualifier *qualifier)
int xDataTOME_copyAttributionList(statusMessageReporting *smr, xDataTOM_attributionList *desc, xDataTOM_element *element)
static int MCGIDI_reaction_ParseReactionTypeAndDetermineProducts(statusMessageReporting *smr, MCGIDI_POPs *pops, MCGIDI_reaction *reaction)
void xDataTOMAL_release(xDataTOM_attributionList *attributes)
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
char const * xDataTOM_getAttributesValueInElement(xDataTOM_element *element, char const *name)
double * ptwX_getPointAtIndex(ptwXPoints *ptwX, int64_t index)
double MCGIDI_target_heated_getTargetMass_MeV(statusMessageReporting *, MCGIDI_target_heated *target)
double MCGIDI_reaction_getProjectileMass_MeV(statusMessageReporting *smr, MCGIDI_reaction *reaction)
double MCGIDI_reaction_getFinalQ(statusMessageReporting *, MCGIDI_reaction *reaction, MCGIDI_quantitiesLookupModes &)
static int MCGIDI_reaction_ParseDetermineReactionProducts(statusMessageReporting *smr, MCGIDI_POPs *pops, MCGIDI_outputChannel *outputChannel, MCGIDI_productsInfo *productsInfo, MCGIDI_reaction *reaction, double *finalQ, int level)
static int MCGIDI_reaction_initialize2(statusMessageReporting *smr, MCGIDI_reaction *reaction)
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
static int MCGIDI_reaction_addReturnProduct(statusMessageReporting *smr, MCGIDI_productsInfo *productsInfo, int ID, MCGIDI_product *product, MCGIDI_reaction *reaction, int transportable)
int MCGIDI_productsInfo_getNumberOfUniqueProducts(MCGIDI_productsInfo *productsInfo)
ptwXYPoints * MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf(statusMessageReporting *smr, xDataTOM_element *linear, char const *toUnits[2])
double ptwXY_getXMin(ptwXYPoints *ptwXY)
int MCGIDI_reaction_getENDL_CSNumbers(MCGIDI_reaction *reaction, int *S)
int MCGIDI_reaction_getENDF_MTNumber(MCGIDI_reaction *reaction)
int MCGIDI_outputChannel_numberOfProducts(MCGIDI_outputChannel *outputChannel)
nfu_status ptwXY_dullEdges(ptwXYPoints *ptwXY, double lowerEps, double upperEps, int positiveXOnly)
MCGIDI_productsInfo * MCGIDI_reaction_getProductsInfo(MCGIDI_reaction *reaction)
static int MCGIDI_reaction_setENDL_CSNumbers(statusMessageReporting *smr, MCGIDI_reaction *reaction)