ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCGIDI_reaction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MCGIDI_reaction.cc
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 
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 );
25 static int MCGIDI_reaction_addReturnProduct( statusMessageReporting *smr, MCGIDI_productsInfo *productsInfo, int ID, MCGIDI_product *product,
26  MCGIDI_reaction *reaction, int transportable );
28 /*
29 ************************************************************
30 */
32 
33  MCGIDI_reaction *reaction;
34 
35  if( ( reaction = (MCGIDI_reaction *) smr_malloc2( smr, sizeof( MCGIDI_reaction ), 0, "reaction" ) ) == NULL ) return( NULL );
36  if( MCGIDI_reaction_initialize( smr, reaction ) ) reaction = MCGIDI_reaction_free( smr, reaction );
37  return( reaction );
38 }
39 /*
40 ************************************************************
41 */
43 
44  if( MCGIDI_reaction_initialize2( smr, reaction ) != 0 ) return( 1 );
45  reaction->transportabilities = new transportabilitiesMap( );
46  return( 0 );
47 }
48 /*
49 ************************************************************
50 */
52 
53  memset( reaction, 0, sizeof( MCGIDI_reaction ) );
54  xDataTOMAL_initial( smr, &(reaction->attributes) );
55  return( 0 );
56 }
57 /*
58 ************************************************************
59 */
61 
62  MCGIDI_reaction_release( smr, reaction );
63  smr_freeMemory( (void **) &reaction );
64  return( NULL );
65 }
66 /*
67 ************************************************************
68 */
70 
71  ptwXY_free( reaction->crossSection );
72  ptwX_free( reaction->crossSectionGrouped );
73  MCGIDI_outputChannel_release( smr, &(reaction->outputChannel) );
74  xDataTOMAL_release( &(reaction->attributes) );
75  smr_freeMemory( (void **) &(reaction->outputChannelStr) );
76  if( reaction->productsInfo.productInfo != NULL ) smr_freeMemory( (void **) &(reaction->productsInfo.productInfo) );
77  delete reaction->transportabilities;
78  MCGIDI_reaction_initialize2( smr, reaction );
79  return( 0 );
80 }
81 /*
82 ************************************************************
83 */
85  MCGIDI_POPs *pops, MCGIDI_reaction *reaction ) {
86 
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" };
91 
92  MCGIDI_reaction_initialize( smr, reaction );
93 
94  reaction->target = target;
96  if( xDataTOME_copyAttributionList( smr, &(reaction->attributes), element ) ) goto err;
97  if( xDataTOME_convertAttributeToInteger( smr, element, "ENDF_MT", &(reaction->ENDF_MT) ) ) goto err;
98  if( ( outputChannelStr = xDataTOM_getAttributesValueInElement( element, "outputChannel" ) ) == NULL ) goto err;
99  if( ( reaction->outputChannelStr = smr_allocateCopyString2( smr, outputChannelStr, "reaction->outputChannelStr" ) ) == NULL ) goto err;
100 
101  if( ( child = xDataTOME_getOneElementByName( smr, element, "crossSection", 1 ) ) == NULL ) goto err;
102  if( ( linear = xDataTOME_getOneElementByName( smr, child, "linear", 0 ) ) == NULL ) {
103  if( ( linear = xDataTOME_getOneElementByName( smr, child, "pointwise", 1 ) ) == NULL ) goto err;
104  }
105  if( xDataTOME_getInterpolation( smr, linear, 0, &independent, &dependent, &qualifier ) ) 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 );
108  goto err;
109  }
110  if( ( reaction->crossSection = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, linear, crossSectionUnits ) ) == NULL ) goto err;
111  reaction->domainValuesPresent = 1;
112  reaction->EMin = ptwXY_getXMin( reaction->crossSection );
113  reaction->EMax = ptwXY_getXMax( reaction->crossSection );
114 
115  if( ( outputChannel = xDataTOME_getOneElementByName( smr, element, "outputChannel", 1 ) ) == NULL ) goto err;
116  if( MCGIDI_outputChannel_parseFromTOM( smr, outputChannel, pops, &(reaction->outputChannel), reaction, NULL ) ) goto err;
117 
118  if( MCGIDI_reaction_ParseReactionTypeAndDetermineProducts( smr, pops, reaction ) != 0 ) goto err;
119 
120  return( 0 );
121 
122 err:
123  MCGIDI_reaction_release( smr, reaction );
124  return( 1 );
125 }
126 /*
127 ************************************************************
128 */
130 
131  MCGIDI_outputChannel *outputChannel = &(reaction->outputChannel);
132  int MT;
133  int particlesChanges[nParticleChanges], numberOfChanges;
134  double finalQ = 0.;
135 
136  if( MCGIDI_reaction_ParseDetermineReactionProducts( smr, pops, outputChannel, &(reaction->productsInfo), reaction, &finalQ, 0 ) != 0 ) return( 1 );
137  reaction->finalQ = finalQ;
138  MT = MCGIDI_reaction_getENDF_MTNumber( reaction );
139  switch( MT ) {
140  case 2 :
142  break;
143  case 18 : case 19 : case 20 : case 21 : case 38 :
145  break;
146  case 102 :
148  break;
149  case 5 :
151  break;
152  default :
153  numberOfChanges = MCGIDI_reaction_particleChanges( reaction->target->projectilePOP, reaction->target->targetPOP, &(reaction->productsInfo),
154  nParticleChanges, particlesChanges );
155 
157  if( numberOfChanges == 0 ) {
159  else {
161  }
162 
163 /*
164  Currently, these are not handled properly:
165  MCGIDI_reactionType_nuclearLevelTransition_e
166  MCGIDI_reactionType_atomic_e
167 */
168  break;
169  }
170 
171  MCGIDI_reaction_setENDL_CSNumbers( smr, reaction );
172  return( 0 );
173 }
174 /*
175 ************************************************************
176 */
177 static int MCGIDI_reaction_particleChanges( MCGIDI_POP *projectile, MCGIDI_POP *target, MCGIDI_productsInfo *productsInfo, int n1, int *particlesChanges ) {
178 
179  int projectileGlobalIndex = projectile->globalPoPsIndex, targetGlobalIndex = target->globalPoPsIndex, i1, i2 = 0;
180  int gammaIndex = PoPs_particleIndex( "gamma" );
181 
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;
185  }
186 
187  for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ ) if( targetGlobalIndex == productsInfo->productInfo[i1].globalPoPsIndex ) break;
188  if( i1 == productsInfo->numberOfProducts ) particlesChanges[i2++] = targetGlobalIndex;
189 
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;
196  }
197  return( i2 );
198 }
199 /*
200 ************************************************************
201 */
203  MCGIDI_productsInfo *productsInfo, MCGIDI_reaction *reaction, double *finalQ, int level ) {
204 /*
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 as both are decay to the groud state (unless a meta-stable is encountered).
207 * Some examples for projectile 'n' and target 'U238' are:
208 * outputChannel products returned during sampling.
209 * 'n + U238' n, U238
210 * 'n + U238 + gamma' n, U238, gamma
211 * 'n + U238_c' n, U238 (even if no gammas are give, return ground state of residual.
212 * 'n + (U238_e3 -> U238 + gamma)' n, U238, gamma
213 */
214  int iProduct, nProducts = MCGIDI_outputChannel_numberOfProducts( outputChannel ), globalPoPsIndex, productIsTrackable;
215  int twoBodyProductsWithData = 0;
216  MCGIDI_product *product;
217  MCGIDI_POP *residual;
218 
219  if( ( level == 0 ) && ( outputChannel->genre == MCGIDI_channelGenre_twoBody_e ) ) {
220  for( iProduct = 0; iProduct < nProducts; iProduct++ ) {
221  product = MCGIDI_outputChannel_getProductAtIndex( smr, outputChannel, 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;
226  }
227  }
228  }
229  if( twoBodyProductsWithData < 0 ) twoBodyProductsWithData = 0;
230  *finalQ += MCGIDI_outputChannel_getQ_MeV( smr, outputChannel, 0 );
231  for( iProduct = 0; iProduct < nProducts; iProduct++ ) {
232  productIsTrackable = twoBodyProductsWithData;
233  product = MCGIDI_outputChannel_getProductAtIndex( smr, outputChannel, iProduct );
234  globalPoPsIndex = product->pop->globalPoPsIndex;
236  productIsTrackable = 1;
237  if( globalPoPsIndex < 0 ) {
238  if( product->distribution.angular != NULL ) {
239  if( product->distribution.angular->type == MCGIDI_angularType_recoil ) productIsTrackable = 0;
240  }
241  if( productIsTrackable ) {
242  int len = (int) strlen( product->pop->name );
243 
244  if( len > 2 ) { /* Special case for continuum reactions with data for residual (e.g., n + U233 -> n + U233_c). */
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;
249  }
250  }
251  if( globalPoPsIndex < 0 ) {
252  smr_setReportError2( smr, smr_unknownID, 1, "product determination for '%s' cannot be determined", product->pop->name );
253  return( 1 );
254  }
255  }
256  }
257  }
258  if( productIsTrackable ) {
259  if( MCGIDI_reaction_addReturnProduct( smr, productsInfo, globalPoPsIndex, product, reaction, 1 ) != 0 ) return( 1 ); }
260  else {
262  if( MCGIDI_reaction_ParseDetermineReactionProducts( smr, pops, &(product->decayChannel), productsInfo, reaction, finalQ, level + 1 ) != 0 ) return( 1 ); }
263  else {
264  *finalQ += product->pop->level_MeV;
265  for( residual = product->pop; residual->globalPoPsIndex < 0; residual = residual->parent ) ;
266  if( MCGIDI_reaction_addReturnProduct( smr, productsInfo, residual->globalPoPsIndex, product, reaction, 0 ) != 0 ) return( 1 );
267  if( product->pop->numberOfGammaBranchs != 0 ) {
268  int gammaIndex = PoPs_particleIndex( "gamma" );
269  if( MCGIDI_reaction_addReturnProduct( smr, productsInfo, gammaIndex, NULL, reaction, 1 ) != 0 ) return( 1 );
270  }
271  }
272  }
273  }
274  return( 0 );
275 }
276 /*
277 ************************************************************
278 */
280  MCGIDI_reaction *reaction, int transportable ) {
281 
282  int i1;
283  enum MCGIDI_productMultiplicityType productMultiplicityType;
284 
285  MCGIDI_misc_updateTransportabilitiesMap2( reaction->transportabilities, ID, transportable );
286  for( i1 = 0; i1 < productsInfo->numberOfProducts; i1++ ) {
287  if( productsInfo->productInfo[i1].globalPoPsIndex == ID ) break;
288  }
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 );
294  }
295  productsInfo->numberOfProducts++;
296  productsInfo->productInfo[i1].globalPoPsIndex = ID;
298  productsInfo->productInfo[i1].multiplicity = 0;
299  productsInfo->productInfo[i1].transportable = transportable;
300  }
301  if( product == NULL ) {
302  productMultiplicityType = MCGIDI_productMultiplicityType_gammaBranching_e; }
303  else {
304  if( ( product->multiplicityVsEnergy != NULL ) || ( product->piecewiseMultiplicities != NULL ) ) {
305  productMultiplicityType = MCGIDI_productMultiplicityType_energyDependent_e; }
306  else {
307  productsInfo->productInfo[i1].multiplicity += product->multiplicity;
308  productMultiplicityType = MCGIDI_productMultiplicityType_integer_e;
309  }
310  }
312  ( productsInfo->productInfo[i1].productMultiplicityType == productMultiplicityType ) ) {
313  productsInfo->productInfo[i1].productMultiplicityType = productMultiplicityType; }
314  else {
316  }
317  return( 0 );
318 }
319 /*
320 ************************************************************
321 */
323 
324  return( reaction->reactionType );
325 }
326 /*
327 ************************************************************
328 */
330 
331  return( reaction->target );
332 }
333 /*
334 ************************************************************
335 */
337 
338  return( MCGIDI_target_heated_getProjectileMass_MeV( smr, reaction->target ) );
339 }
340 /*
341 ************************************************************
342 */
344 
345  return( MCGIDI_target_heated_getTargetMass_MeV( smr, reaction->target ) );
346 }
347 /*
348 ************************************************************
349 */
350 int MCGIDI_reaction_getDomain( statusMessageReporting * /*smr*/, MCGIDI_reaction *reaction, double *EMin, double *EMax ) {
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 ) return( -1 );
359  *EMin = reaction->EMin;
360  *EMax = reaction->EMax;
361  return( 0 );
362 }
363 /*
364 ************************************************************
365 */
366 int MCGIDI_reaction_fixDomains( statusMessageReporting * /*smr*/, MCGIDI_reaction *reaction, double EMin, double EMax, nfu_status *status ) {
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. ) ) return( 0 );
373 
374  *status = ptwXY_dullEdges( reaction->crossSection, lowerEps, upperEps, 1 );
375  return( *status != nfu_Okay );
376 }
377 /*
378 ************************************************************
379 */
381  bool sampling ) {
382 
383  double e_in = modes.getProjectileEnergy( ), xsec;
384 
386  if( e_in < reaction->EMin ) e_in = reaction->EMin;
387  if( e_in > reaction->EMax ) e_in = reaction->EMax;
388  ptwXY_getValueAtX( reaction->crossSection, e_in, &xsec ); }
390  int index = modes.getGroupIndex( );
391  double *xSecP = ptwX_getPointAtIndex( reaction->crossSectionGrouped, index );
392 
393  if( xSecP != NULL ) {
394  xsec = *xSecP;
395  if( sampling && ( index == reaction->thresholdGroupIndex ) ) xsec += reaction->thresholdGroupedDeltaCrossSection; }
396  else {
397  xsec = 0.;
398  smr_setReportError2( smr, smr_unknownID, 1, "Invalid cross section group index %d", index );
399  } }
400  else {
401  xsec = 0.;
402  }
403  return( xsec );
404 }
405 /*
406 ************************************************************
407 */
409 
410  return( reaction->finalQ );
411 }
412 /*
413 ************************************************************
414 */
416 
417  return( reaction->ENDF_MT );
418 }
419 /*
420 ************************************************************
421 */
423 
424  if( S != NULL ) *S = reaction->ENDL_S;
425  return( reaction->ENDL_C );
426 }
427 /*
428 ************************************************************
429 */
431 
432  int MT = MCGIDI_reaction_getENDF_MTNumber( reaction );
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 };
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 - 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 ) {
464  if( MT < 650 ) {
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;
482  }
483  }
484  return( 0 );
485 }
486 /*
487 ************************************************************
488 */
490 
491  return( &(reaction->productsInfo) );
492 }
493 /*
494 ************************************************************
495 */
497  GIDI_settings_particle const *projectileSettings, double temperature_MeV, ptwXPoints *totalGroupedCrossSection ) {
498 
499  if( totalGroupedCrossSection != NULL ) {
500  nfu_status status_nf;
501  GIDI_settings_group group( projectileSettings->getGroup( ) );
502 
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 );
506 
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 ) {
512  /* factor 2 for linear reject in bin but above threshold. */
514  ( 2 * ( group[reaction->thresholdGroupIndex+1] - group[reaction->thresholdGroupIndex] ) / reaction->thresholdGroupDomain - 1 );
515  }
516  }
517  }
518  return( 0 );
519 }
520 
521 /*
522 *********************** productsInfo ***********************
523 */
524 /*
525 ************************************************************
526 */
528 
529  return( productsInfo->numberOfProducts );
530 }
531 /*
532 ************************************************************
533 */
535 
536  if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) ) return( -1 );
537  return( productsInfo->productInfo[index].globalPoPsIndex );
538 }
539 /*
540 ************************************************************
541 */
543 
544  if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) ) return( MCGIDI_productMultiplicityType_invalid_e );
545  return( productsInfo->productInfo[index].productMultiplicityType );
546 }
547 /*
548 ************************************************************
549 */
551 
552  if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) ) return( -1 );
553  return( productsInfo->productInfo[index].multiplicity );
554 }
555 /*
556 ************************************************************
557 */
559 
560  if( ( index < 0 ) || ( index >= productsInfo->numberOfProducts ) ) return( -1 );
561  return( productsInfo->productInfo[index].transportable );
562 }
563 
564 #if defined __cplusplus
565 }
566 #endif
567