ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCGIDI_product.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MCGIDI_product.cc
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 #include <string.h>
6 #include <cmath>
7 
8 #include "MCGIDI.h"
9 #include "MCGIDI_misc.h"
10 #include "MCGIDI_fromTOM.h"
11 
12 #if defined __cplusplus
13 namespace GIDI {
14 using namespace GIDI;
15 #endif
16 
17 typedef struct polynomialCallbackArgs_s {
18  int length;
19  double energyFactor;
20  double *coefficients;
22 
26  ptwXYPoints **multiplicityVsEnergy, ptwXYPoints **norms );
27 static double MCGIDI_product_evaluatePolynomial( double x, polynomialCallbackArgs *args );
28 /*
29 ************************************************************
30 */
32 
33  MCGIDI_product *product;
34 
35  if( ( product = (MCGIDI_product *) smr_malloc2( smr, sizeof( MCGIDI_product ), 0, "product" ) ) == NULL ) return( NULL );
36  if( MCGIDI_product_initialize( smr, product ) ) product = MCGIDI_product_free( smr, product );
37  return( product );
38 }
39 /*
40 ************************************************************
41 */
43 
44  memset( product, 0, sizeof( MCGIDI_product ) );
45  product->delayedNeutronIndex = -1;
46  return( 0 );
47 }
48 /*
49 ************************************************************
50 */
52 
53  MCGIDI_product_release( smr, product );
54  smr_freeMemory( (void **) &product );
55  return( NULL );
56 }
57 /*
58 ************************************************************
59 */
61 
62  int i;
63 
64  if( product->label != NULL ) smr_freeMemory( (void **) &(product->label) );
65 
66  if( product->multiplicityVsEnergy != NULL ) ptwXY_free( product->multiplicityVsEnergy );
67  if( product->piecewiseMultiplicities != NULL ) {
68  for( i = 0; i < product->numberOfPiecewiseMultiplicities; i++ ) ptwXY_free( product->piecewiseMultiplicities[i] );
69  smr_freeMemory( (void **) &(product->piecewiseMultiplicities) );
70  }
71  if( product->norms != NULL ) ptwXY_free( product->norms );
72 
73  MCGIDI_distribution_release( smr, &(product->distribution) );
74  MCGIDI_outputChannel_release( smr, &(product->decayChannel) );
75 
76  MCGIDI_product_initialize( smr, product );
77  return( 0 );
78 }
79 /*
80 ************************************************************
81 */
83  MCGIDI_POPs *pops, MCGIDI_product *product, int *delayedNeutronIndex ) {
84 
85  char const *name, *label, *delayedNeutron, *multiplicityStr, *multiplicityUnits[2] = { "MeV", "" };
86  xDataTOM_element *multiplicity, *multiplicityTOM, *decayChannelElement;
87  nfu_status status;
88  ptwXYPoints *multiplicityVsEnergy = NULL, *norms1 = NULL, *norms2 = NULL;
89 
90  MCGIDI_product_initialize( smr, product );
91 
92  product->outputChannel = outputChannel;
93  if( ( name = xDataTOM_getAttributesValueInElement( element, "name" ) ) == NULL ) goto err;
94  if( ( product->pop = MCGIDI_POPs_findParticle( pops, name ) ) == NULL ) {
95  smr_setReportError2( smr, smr_unknownID, 1, "product '%s' not found in pops", name );
96  goto err;
97  }
98  if( ( label = xDataTOM_getAttributesValueInElement( element, "label" ) ) != NULL ) {
99  if( ( product->label = smr_allocateCopyString2( smr, label, "product->label" ) ) == NULL ) goto err;
100  }
101 
102  if( ( delayedNeutron = xDataTOM_getAttributesValueInElement( element, "emissionMode" ) ) != NULL ) {
103  if( strcmp( delayedNeutron, "delayed" ) == 0 ) {
104  if( ( delayedNeutron = xDataTOM_getAttributesValueInElement( element, "decayRate" ) ) == NULL ) {
105  goto err;
106  }
107  if( MCGIDI_misc_PQUStringToDoubleInUnitOf( smr, delayedNeutron, "1/s", &(product->delayedNeutronRate) ) != 0 ) goto err;
108  product->delayedNeutronIndex = *delayedNeutronIndex;
109  (*delayedNeutronIndex)++;
110  }
111  }
112 
113  if( ( multiplicityStr = xDataTOM_getAttributesValueInElement( element, "multiplicity" ) ) == NULL ) goto err;
114  if( xDataTOME_convertAttributeToInteger( NULL, element, "multiplicity", &(product->multiplicity) ) ) {
115  if( strcmp( multiplicityStr, "energyDependent" ) ) {
116  smr_setReportError2( smr, smr_unknownID, 1, "invalid multiplicity '%s' for product '%s'", multiplicityStr, name );
117  goto err;
118  }
119  if( ( multiplicity = xDataTOME_getOneElementByName( smr, element, "multiplicity", 1 ) ) == NULL ) goto err;
120  if( ( multiplicityTOM = xDataTOME_getOneElementByName( NULL, multiplicity, "weightedReference", 0 ) ) != NULL ) {
121  if( MCGIDI_product_parseWeightedReferenceMultiplicityFromTOM( smr, multiplicityTOM, product, &multiplicityVsEnergy, &norms1 ) ) goto err; }
122  else if( ( multiplicityTOM = xDataTOME_getOneElementByName( NULL, multiplicity, "piecewise", 0 ) ) != NULL ) {
123  if( MCGIDI_product_parsePiecewiseMultiplicity( smr, multiplicityTOM, product ) ) goto err; }
124  else if( ( multiplicityTOM = xDataTOME_getOneElementByName( NULL, multiplicity, "polynomial", 0 ) ) != NULL ) {
125  if( ( multiplicityVsEnergy = MCGIDI_product_parsePolynomialMultiplicity( smr, multiplicityTOM, product ) ) == NULL ) goto err; }
126  else {
127 /* ??????? Need to check interpolation. */
128  if( ( multiplicityTOM = xDataTOME_getOneElementByName( smr, multiplicity, "pointwise", 1 ) ) == NULL ) goto err;
129  if( ( multiplicityVsEnergy = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, multiplicityTOM, multiplicityUnits ) ) == NULL ) goto err;
130  }
131  }
132 
133  if( strcmp( product->pop->name, "gamma" ) == 0 ) {
134  if( ( norms2 = ptwXY_new( ptwXY_interpolationLinLin, NULL, 2., 1e-3, 200, 10, &status, 0 ) ) == NULL ) {
135  smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_new err = %d: %s\n", status, nfu_statusMessage( status ) );
136  goto err;
137  }
138  }
139  if( MCGIDI_distribution_parseFromTOM( smr, element, product, pops, norms2 ) ) goto err;
140  if( norms2 != NULL ) {
141  if( ptwXY_length( norms2 ) < 2 ) {
142  norms2 = ptwXY_free( norms2 ); }
143  else {
144  if( ptwXY_simpleCoalescePoints( norms2 ) != nfu_Okay ) goto err;
145  if( ( ptwXY_getYMin( norms2 ) > 0.99 ) && ( ptwXY_getYMax( norms2 ) < 1.01 ) ) norms2 = ptwXY_free( norms2 );
146  }
147  }
148  if( ( norms1 != NULL ) && ( norms2 != NULL ) ) {
149  smr_setReportError2p( smr, smr_unknownID, 1, "norm1 and norm2 are both not NULL" );
150  goto err;
151  }
152 
153  product->multiplicityVsEnergy = multiplicityVsEnergy;
154  product->norms = norms1;
155  if( norms2 != NULL ) product->norms = norms2;
156 
157  if( ( decayChannelElement = xDataTOME_getOneElementByName( NULL, element, "decayChannel", 0 ) ) != NULL ) {
158  if( MCGIDI_outputChannel_parseFromTOM( smr, decayChannelElement, pops, &(product->decayChannel), NULL, product ) ) goto err;
159  }
160 
161  return( 0 );
162 
163 err:
164  if( multiplicityVsEnergy != NULL ) ptwXY_free( multiplicityVsEnergy );
165  if( norms1 != NULL ) ptwXY_free( norms1 );
166  if( norms2 != NULL ) ptwXY_free( norms2 );
167  MCGIDI_product_release( smr, product );
168  return( 1 );
169 }
170 /*
171 ************************************************************
172 */
174 
175  int i;
176  xDataTOM_XYs *XYs;
177  xDataTOM_regionsXYs *regionsXYs = (xDataTOM_regionsXYs *) element->xDataInfo.data;
178  ptwXYPoints *multiplicityVsEnergy;
179  char const *multiplicityUnits[2] = { "MeV", "" };
180 
181  if( ( product->piecewiseMultiplicities = (ptwXYPoints **) smr_malloc2( smr, regionsXYs->length * sizeof( ptwXYPoints * ), 1, "piecewiseMultiplicities" ) ) == NULL ) return( 1 );
182 
183  for( i = 0; i < regionsXYs->length; i++ ) {
184 /* ??????? Need to check interpolation. */
185  XYs = &(regionsXYs->XYs[i]);
186  if( ( multiplicityVsEnergy = MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( smr, XYs, ptwXY_interpolationLinLin, multiplicityUnits )
187  ) == NULL ) return( 1 );
188  product->piecewiseMultiplicities[i] = multiplicityVsEnergy;
189  product->numberOfPiecewiseMultiplicities++;
190  }
191 
192  return( 0 );
193 }
194 /*
195 ************************************************************
196 */
198 
199  int length;
200  double *coefficients;
201  char const *energyUnit;
202  ptwXYPoints *ptwXY = NULL;
203  nfu_status status;
204  double EMin, EMax;
206 
207  if( MCGIDI_product_getDomain( smr, product, &EMin, &EMax ) ) goto err;
208 
209  length = xDataTOM_polynomial_getDataFromXDataInfo( (xDataTOM_xDataInfo *) &(element->xDataInfo), &coefficients );
210  if( ( ptwXY = ptwXY_new( ptwXY_interpolationLinLin, NULL, 2., 1e-3, length, 10, &status, 0 ) ) == NULL ) {
211  smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_new err = %d: %s\n", status, nfu_statusMessage( status ) );
212  goto err;
213  }
214 
215  if( ( energyUnit = xDataTOM_axes_getUnit( smr, &(element->xDataInfo.axes), 0 ) ) == NULL ) goto err;
216  args.energyFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyUnit, "MeV" );
217  if( !smr_isOk( smr ) ) goto err;
218 
219  args.length = length;
220  args.coefficients = coefficients;
221  ptwXY_setValueAtX( ptwXY, EMin, MCGIDI_product_evaluatePolynomial( EMin, &args ) );
222  ptwXY_setValueAtX( ptwXY, EMax, MCGIDI_product_evaluatePolynomial( EMax, &args ) );
223  if( length > 2 ) { /* ?????????????? This needs work. */
224  int i, n = 4 * length;
225  double E = EMin, dE = ( EMax - EMin ) / n;
226 
227  for( i = 1; i < n; i++ ) {
228  E += dE;
229  ptwXY_setValueAtX( ptwXY, E, MCGIDI_product_evaluatePolynomial( E, &args ) );
230  }
231  }
232  return( ptwXY );
233 
234 err:
235  if( ptwXY != NULL ) ptwXY_free( ptwXY );
236  return( NULL );
237 }
238 /*
239 ************************************************************
240 */
242 
243  int i;
244  double value = 0.;
245 
246  x /= args->energyFactor;
247  for( i = args->length; i > 0; i-- ) value = value * x + args->coefficients[i-1];
248 
249  return( value );
250 }
251 /*
252 ************************************************************
253 */
255  ptwXYPoints **multiplicityVsEnergy, ptwXYPoints **norms ) {
256 
257  char const *link, *energyInMWUnits[2] = { "MeV", "" };
258  xDataTOM_element *reference, *productTOM, *multiplicity, *weights, *pointwise;
259 
260  if( ( reference = xDataTOME_getOneElementByName( smr, element, "reference", 1 ) ) == NULL ) goto err;
261  if( ( link = xDataTOM_getAttributesValueInElement( reference, "xlink:href" ) ) == NULL ) goto err;
262  if( ( productTOM = xDataTOM_getLinksElement( smr, reference, link ) ) == NULL ) goto err;
263  if( ( multiplicity = xDataTOME_getOneElementByName( smr, productTOM, "multiplicity", 1 ) ) == NULL ) goto err;
264  /* Currently, only pointwise supported. */
265  if( ( pointwise = xDataTOME_getOneElementByName( smr, multiplicity, "pointwise", 1 ) ) == NULL ) goto err;
266  if( ( *multiplicityVsEnergy = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, pointwise, energyInMWUnits ) ) == NULL ) goto err;
267 
268  if( ( weights = xDataTOME_getOneElementByName( smr, element, "weights", 1 ) ) == NULL ) goto err;
269  if( ( pointwise = xDataTOME_getOneElementByName( smr, weights, "pointwise", 1 ) ) == NULL ) goto err;
270  if( ( *norms = MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( smr, pointwise, energyInMWUnits ) ) == NULL ) goto err;
271 
272  return( 0 );
273 
274 err:
275  if( *multiplicityVsEnergy != NULL ) *multiplicityVsEnergy = ptwXY_free( *multiplicityVsEnergy );
276  if( *norms != NULL ) *norms = ptwXY_free( *norms );
277  return( 1 );
278 }
279 /*
280 ************************************************************
281 */
282 int MCGIDI_product_getDomain( statusMessageReporting *smr, MCGIDI_product *product, double *EMin, double *EMax ) {
283 
284  return( MCGIDI_outputChannel_getDomain( smr, product->outputChannel, EMin, EMax ) );
285 }
286 /*
287 ************************************************************
288 */
289 int MCGIDI_product_setTwoBodyMasses( statusMessageReporting *smr, MCGIDI_product *product, double projectileMass_MeV, double targetMass_MeV,
290  double productMass_MeV, double residualMass_MeV ) {
291 
292  return( MCGIDI_angular_setTwoBodyMasses( smr, product->distribution.angular, projectileMass_MeV, targetMass_MeV, productMass_MeV, residualMass_MeV ) );
293 }
294 /*
295 ************************************************************
296 */
298 
299  return( MCGIDI_POP_getMass_MeV( product->pop ) );
300 }
301 /*
302 ************************************************************
303 */
305 
306  return( MCGIDI_outputChannel_getTargetHeated( smr, product->outputChannel ) );
307 }
308 /*
309 ************************************************************
310 */
312 
314 }
315 /*
316 ************************************************************
317 */
319 
320  return( MCGIDI_outputChannel_getTargetMass_MeV( smr, product->outputChannel ) );
321 }
322 /*
323 ************************************************************
324 */
325 int MCGIDI_product_sampleMultiplicity( statusMessageReporting * /*smr*/, MCGIDI_product *product, double e_in, double r ) {
326 
327  int i, multiplicity;
328  double y, norm = 1.0;
329  ptwXYPoints *ptwXY = product->multiplicityVsEnergy;
330 
331  if( product->piecewiseMultiplicities != NULL ) {
332  for( i = 0; i < product->numberOfPiecewiseMultiplicities - 1; i++ ) {
333  if( e_in < ptwXY_getXMax( product->piecewiseMultiplicities[i] ) ) break;
334  }
335  ptwXY = product->piecewiseMultiplicities[i];
336  }
337  y = MCGIDI_sampling_ptwXY_getValueAtX( ptwXY, e_in );
338  if( product->norms != NULL ) norm = MCGIDI_sampling_ptwXY_getValueAtX( product->norms, e_in );
339  y *= norm;
340  multiplicity = (int) y;
341  if( r < ( y - multiplicity ) ) multiplicity++;
342 
343  return( multiplicity );
344 }
345 /*
346 ************************************************************
347 */
349  MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
350 
352  smr_setReportError2( smr, smr_unknownID, 1, "product distribution is not angular: type = %d", product->distribution.type );
353  return( 1 );
354  }
355  return( MCGIDI_angular_sampleMu( smr, product->distribution.angular, modes, decaySamplingInfo ) );
356 }
357 
358 
359 /*
360 ************************************************************
361 */
362 int MCGIDI_sampledProducts_initialize( statusMessageReporting *smr, MCGIDI_sampledProductsDatas *sampledProductsDatas, int incrementSize ) {
363 
364  if( incrementSize < 10 ) incrementSize = 10;
365  sampledProductsDatas->numberOfProducts = 0;
366  sampledProductsDatas->numberAllocated = 0;
367  sampledProductsDatas->incrementSize = incrementSize;
368  sampledProductsDatas->products = NULL;
369  return( MCGIDI_sampledProducts_remalloc( smr, sampledProductsDatas ) );
370 }
371 /*
372 ************************************************************
373 */
375 
376  smr_freeMemory( (void **) &(sampledProductsDatas->products) );
377  return( 0 );
378 }
379 /*
380 ************************************************************
381 */
383 
384  int size = sampledProductsDatas->numberAllocated + sampledProductsDatas->incrementSize;
385 
386  if( ( sampledProductsDatas->products = (MCGIDI_sampledProductsData *) smr_realloc2( smr, sampledProductsDatas->products,
387  size * sizeof( MCGIDI_sampledProductsData ), "products" ) ) != NULL ) {
388  sampledProductsDatas->numberAllocated = size;
389  return( 0 );
390  }
391  sampledProductsDatas->numberOfProducts = 0;
392  sampledProductsDatas->numberAllocated = 0;
393  return( 1 );
394 }
395 /*
396 ************************************************************
397 */
399 
400  if( sampledProductsDatas->numberOfProducts == sampledProductsDatas->numberAllocated ) {
401  if( ( MCGIDI_sampledProducts_remalloc( smr, sampledProductsDatas ) ) != 0 ) return( 1 );
402  }
403  sampledProductsDatas->products[sampledProductsDatas->numberOfProducts] = *sampledProductsData;
404  sampledProductsDatas->numberOfProducts++;
405  return( 0 );
406 }
407 /*
408 ************************************************************
409 */
411 
412  return( sampledProductsDatas->numberOfProducts );
413 }
414 /*
415 ************************************************************
416 */
418 
419  if( index < 0 ) return( NULL );
420  if( index >= sampledProductsDatas->numberOfProducts ) return( NULL );
421  return( &(sampledProductsDatas->products[index]) );
422 }
423 
424 #if defined __cplusplus
425 }
426 #endif
427