Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
tpia_product.cc
Go to the documentation of this file.
1 /*
2 # <<BEGIN-copyright>>
3 # Copyright (c) 2010, Lawrence Livermore National Security, LLC.
4 # Produced at the Lawrence Livermore National Laboratory
5 # Written by Bret R. Beck, beck6@llnl.gov.
6 # CODE-461393
7 # All rights reserved.
8 #
9 # This file is part of GIDI. For details, see nuclear.llnl.gov.
10 # Please also read the "Additional BSD Notice" at nuclear.llnl.gov.
11 #
12 # Redistribution and use in source and binary forms, with or without modification,
13 # are permitted provided that the following conditions are met:
14 #
15 # 1) Redistributions of source code must retain the above copyright notice,
16 # this list of conditions and the disclaimer below.
17 # 2) Redistributions in binary form must reproduce the above copyright notice,
18 # this list of conditions and the disclaimer (as noted below) in the
19 # documentation and/or other materials provided with the distribution.
20 # 3) Neither the name of the LLNS/LLNL nor the names of its contributors may be
21 # used to endorse or promote products derived from this software without
22 # specific prior written permission.
23 #
24 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
25 # EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
26 # OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
27 # SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR
28 # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
29 # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
30 # OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
31 # AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
32 # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
33 # EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34 # <<END-copyright>>
35 */
36 #include <string.h>
37 #include <ctype.h>
38 
39 #include <gString.h>
40 #include <tpia_target.h>
41 #include <tpia_misc.h>
42 
43 #if defined __cplusplus
44 namespace GIDI {
45 using namespace GIDI;
46 #endif
47 
48 static const int tpia_b_unknown = 0,
49  tpia_b_twoBody_angular = tpia_m_angular,
50  tpia_b_twoBody_formFactor = 0, /* ??? */
51  tpia_b_NBody_Legendre = tpia_m_Legendre,
52  tpia_b_NBody_angular_energy = tpia_m_angular | tpia_m_angular_energy,
53  tpia_b_NBody_uncorrelate_Legendre = tpia_m_angular | tpia_m_Legendre,
54  tpia_b_NBody_pairProduction = 0; /* ??? */
55 
56 const char *tpia_productGenre_unknown = "unknown",
57  *tpia_productGenre_twoBody_angular = "twoBody_angular",
58  *tpia_productGenre_twoBody_formFactor = "twoBody_formFactor",
59  *tpia_productGenre_NBody_Legendre = "NBody_Legendre",
60  *tpia_productGenre_NBody_angular_energy = "NBody_angular_energy",
61  *tpia_productGenre_NBody_uncorrelate_Legendre = "NBody_uncorrelate_Legendre",
62  *tpia_productGenre_NBody_pairProduction = "NBody_pairProduction";
63 
64 static int _tpia_product_getProductOutgoingData( statusMessageReporting *smr, xData_element *productElement, tpia_product *product );
65 static int _tpia_product_checkRequiredData( statusMessageReporting *smr, int allowMany, int m, xData_element *productElement, tpia_product *product, char *str );
66 static int _tpia_product_getDepositionEnergy( statusMessageReporting *smr, xData_element *depositionEnergy, tpia_product *product );
67 static int _tpia_product_getMultiplicityFromElement( statusMessageReporting *smr, xData_element *data, tpia_product *product );
68 /*
69 ************************************************************
70 */
72 
73  tpia_product *product;
74 
75  //if( ( product = xData_malloc2( smr, sizeof( tpia_product ), 0, "product" ) ) == NULL ) return( NULL );
76  if( ( product = (tpia_product*) xData_malloc2( smr, sizeof( tpia_product ), 0, "product" ) ) == NULL ) return( NULL );
77  if( tpia_product_initialize( smr, product ) ) product = tpia_product_free( smr, product );
78  return( product );
79 }
80 /*
81 ************************************************************
82 */
84 
85  memset( product, 0, sizeof( tpia_product ) );
86  if( tpia_angular_initialize( smr, &(product->angular) ) ) return( 1 );
87  if( tpia_Legendre_initialize( smr, &(product->Legendre) ) ) return( 1 );
88  return( 0 );
89 }
90 /*
91 ************************************************************
92 */
94  xData_element *productElement ) {
95 
96  tpia_product *product;
97 
98  if( ( product = tpia_product_create( smr ) ) == NULL ) return( NULL );
99  if( tpia_product_getFromElement( smr, channel, parentProduct, productElement, product ) != 0 ) product = tpia_product_free( smr, product );
100  return( product );
101 }
102 /*
103 ************************************************************
104 */
106 
107  tpia_product_release( smr, product );
108  xData_free( smr, product );
109  return( NULL );
110 }
111 /*
112 ************************************************************
113 */
115 
116  tpia_multiplicity *multiplicity, *multiplicity_next;
117  tpia_product *decayProduct, *nextProduct;
118 
119  xData_releaseAttributionList( smr, &(product->attributes) );
120  //product->depositionEnergyGrouped.data = xData_free( smr, product->depositionEnergyGrouped.data );
121  product->depositionEnergyGrouped.data = (double*) xData_free( smr, product->depositionEnergyGrouped.data );
122 
123  if( product->multiplicityVsEnergy != NULL ) tpia_multiplicity_free( smr, product->multiplicityVsEnergy );
124  for( multiplicity = product->delayedNeutronMultiplicityVsEnergy; multiplicity != NULL; multiplicity = multiplicity_next ) {
125  multiplicity_next = multiplicity->next;
126  tpia_multiplicity_free( smr, multiplicity );
127  }
128  tpia_angular_release( smr, &(product->angular) );
129  tpia_Legendre_release( smr, &(product->Legendre ) );
130  tpia_angularEnergy_release( smr, &(product->angularEnergy) );
131  for( decayProduct = product->decayChannel.products; decayProduct != NULL; decayProduct = nextProduct ) {
132  nextProduct = decayProduct->next;
133  tpia_product_free( smr, decayProduct );
134  }
135  product->decayChannel.numberOfProducts = 0;
136  product->decayChannel.products = NULL;
137  return( 0 );
138 }
139 /*
140 ************************************************************
141 */
142 int tpia_product_getFromElement( statusMessageReporting *smr, tpia_channel *channel, tpia_product *parentProduct, xData_element *productElement,
143  tpia_product *product ) {
144 
145  char const *productGenre;
146  char *name, *multiplicity, *e;
147 
148  xData_addToAccessed( smr, productElement, 1 );
149  product->channel = channel;
150  product->parentProduct = parentProduct;
151  if( xData_copyAttributionList( smr, &(product->attributes), &(productElement->attributes) ) != 0 ) return( 0 );
152  name = tpia_misc_pointerToAttributeIfAllOk2( smr, productElement, 1, &(product->attributes), "particle" );
153  if( name != NULL ) {
154  product->productID = tpia_particle_getInternalID( smr, name );
155  multiplicity = tpia_misc_pointerToAttributeIfAllOk2( smr, productElement, 1, &(product->attributes), "multiplicity" );
156  if( multiplicity != NULL ) {
157  if( strcmp( multiplicity, "energyDependent" ) && strcmp( multiplicity, "partialProduction" ) ) { /* Must be an integer. */
158  product->multiplicity = strtol( multiplicity, &e, 10 );
159  while( isspace( *e ) ) e++;
160  if( *e != 0 ) tpia_misc_setMessageError_Element( smr, NULL, productElement, __FILE__, __LINE__, 1, "bad multiplicity = %s", multiplicity );
161  }
162  }
163  }
164  if( ( productGenre = tpia_misc_pointerToAttributeIfAllOk2( smr, productElement, 1, &(product->attributes), "genre" ) ) != NULL ) {
165  if( strcmp( productGenre, tpia_productGenre_unknown ) == 0 ) {
166  product->b_dataRequired = 0;
167  product->genre = tpia_productGenre_unknown; }
168  else if( strcmp( productGenre, tpia_productGenre_twoBody_angular ) == 0 ) {
169  product->b_dataRequired = tpia_b_twoBody_angular;
171  else if( strcmp( productGenre, tpia_productGenre_twoBody_formFactor ) == 0 ) {
172  product->b_dataRequired = tpia_b_twoBody_formFactor;
174  else if( strcmp( productGenre, tpia_productGenre_NBody_Legendre ) == 0 ) {
175  product->b_dataRequired = tpia_b_NBody_Legendre;
177  else if( strcmp( productGenre, tpia_productGenre_NBody_angular_energy ) == 0 ) {
178  product->b_dataRequired = tpia_b_NBody_angular_energy;
180  else if( strcmp( productGenre, tpia_productGenre_NBody_uncorrelate_Legendre ) == 0 ) {
181  product->b_dataRequired = tpia_b_NBody_uncorrelate_Legendre;
183  else if( strcmp( productGenre, tpia_productGenre_NBody_pairProduction ) == 0 ) {
184  product->b_dataRequired = tpia_b_NBody_pairProduction;
186  else {
187  tpia_misc_setMessageError_Element( smr, NULL, productElement, __FILE__, __LINE__, 1, "unsupported product genre = %s", productGenre );
188  }
189  if( smr_isOk( smr ) ) _tpia_product_getProductOutgoingData( smr, productElement, product );
190  }
191  return( !smr_isOk( smr ) );
192 }
193 /*
194 ************************************************************
195 */
196 static int _tpia_product_getProductOutgoingData( statusMessageReporting *smr, xData_element *productElement, tpia_product *product ) {
197 
199  int allowMany = 0;
200 
201  for( data = xData_getFirstElement( productElement ); data != NULL; data = xData_getNextElement( data ) ) {
202  if( strcmp( data->name, "depositionEnergy" ) == 0 ) {
203  //if( _tpia_product_checkRequiredData( smr, allowMany, tpia_m_depositionEnergy, productElement, product, "deposition energy" ) ) return( 1 );
204  if( _tpia_product_checkRequiredData( smr, allowMany, tpia_m_depositionEnergy, productElement, product, (char*)"deposition energy" ) ) return( 1 );
205  if( _tpia_product_getDepositionEnergy( smr, data, product ) != 0 ) return( 1 ); }
206  else if( strcmp( data->name, "multiplicity" ) == 0 ) {
207  allowMany = ( product->channel->fission != NULL ) && ( strcmp( product->productID->name, "n_1" ) == 0 );
208  //if( _tpia_product_checkRequiredData( smr, allowMany, tpia_m_multiplicity, productElement, product, "multiplicity" ) ) return( 1 );
209  if( _tpia_product_checkRequiredData( smr, allowMany, tpia_m_multiplicity, productElement, product, (char*) "multiplicity" ) ) return( 1 );
210  if( _tpia_product_getMultiplicityFromElement( smr, data, product ) != 0 ) return( 1 ); }
211  else if( strcmp( data->name, "angular" ) == 0 ) {
212  //if( _tpia_product_checkRequiredData( smr, allowMany, tpia_m_angular, productElement, product, "angular" ) ) return( 1 );
213  if( _tpia_product_checkRequiredData( smr, allowMany, tpia_m_angular, productElement, product, (char*) "angular" ) ) return( 1 );
214  if( tpia_angular_getFromElement( smr, data, &(product->angular) ) != 0 ) return( 1 ); }
215  else if( strcmp( data->name, "Legendre" ) == 0 ) {
216  //if( _tpia_product_checkRequiredData( smr, allowMany, tpia_m_Legendre, productElement, product, "Legendre" ) ) return( 1 );
217  if( _tpia_product_checkRequiredData( smr, allowMany, tpia_m_Legendre, productElement, product, (char*) "Legendre" ) ) return( 1 );
218  if( tpia_Legendre_getFromElement( smr, data, &(product->Legendre) ) != 0 ) return( 1 ); }
219  else if( strcmp( data->name, "angularEnergy" ) == 0 ) {
220  if( _tpia_product_checkRequiredData( smr, allowMany, tpia_m_angular_energy, productElement, product, (char*) "angularEnergy" ) ) return( 1 );
221  if( tpia_angularEnergy_getFromElement( smr, data, &(product->angularEnergy) ) != 0 ) return( 1 ); }
222  else if( strcmp( data->name, "decayChannel" ) == 0 ) {
223  xData_addToAccessed( smr, data, 1 );
224  if( tpia_product_getDecayChannelFromElement( smr, data, product->channel, product, &(product->decayChannel.products) ) ) return( 1 ); }
225  else {
226  printf( " %s\n", data->name );
227  }
228  }
229  if( ( product->b_dataPresent >> tpia_m_commonShift ) != ( product->b_dataRequired >> tpia_m_commonShift ) ) {
230  gString gStr;
231  int missing = ~product->b_dataPresent & product->b_dataRequired;
232  char const *str = "";
233  if( gString_initialize( NULL, &gStr, 100, 100 ) == 0 ) {
234  if( missing & tpia_m_angular ) gString_addTo( NULL, &gStr, "angular " );
235  if( missing & tpia_m_formFactor ) gString_addTo( NULL, &gStr, "formFactor " );
236  if( missing & tpia_m_Legendre ) gString_addTo( NULL, &gStr, "Legendre " );
237  if( missing & tpia_m_angular_energy ) gString_addTo( NULL, &gStr, "angular_energy " );
238  str = gString_string( NULL, &gStr );
239  }
240  tpia_misc_setMessageError_Element( smr, NULL, productElement, __FILE__, __LINE__, 1, "missing data %s for product %s", str,
241  product->productID->name );
242  gString_release( NULL, &gStr );
243  return( 1 );
244  }
245  return( 0 );
246 }
247 /*
248 ************************************************************
249 */
250 static int _tpia_product_checkRequiredData(statusMessageReporting *smr, int allowMany, int m, xData_element *productElement, tpia_product *product, char *str) {
251 
252  if( !allowMany && ( product->b_dataPresent & m ) ) {
253  tpia_misc_setMessageError_Element( smr, NULL, productElement, __FILE__, __LINE__, 1, "multiple %s", str );
254  return( 1 );
255  }
257  if( ( product->b_dataRequired & m ) == 0 ) {
258  tpia_misc_setMessageError_Element( smr, NULL, productElement, __FILE__, __LINE__, 1, "extra product data %s", str );
259  return( 1 );
260  }
261  }
262  product->b_dataPresent += m;
263  return( 0 );
264 }
265 /*
266 ************************************************************
267 */
269  tpia_product **priorProductNext ) {
270 
271  xData_elementList *list;
272  tpia_product *product;
273  int i, status = 0;
274 
275  list = xData_getElementsByTagName( smr, parentElement, "product" );
276  for( i = 0; i < list->n; i++ ) {
277  if( ( product = tpia_product_createGetFromElement( smr, channel, parentProduct, list->items[i].element ) ) == NULL ) {
278  status = 1;
279  break;
280  }
281  if( parentProduct == NULL ) {
284  channel->decayChannel.numberOfProducts++; }
285  else {
286  channel->decayChannel.m1_fullMass_MeV = parentProduct->productID->fullMass_MeV;
287  channel->decayChannel.m2_fullMass_MeV = 0.;
288  parentProduct->decayChannel.numberOfProducts++;
289  }
290  *priorProductNext = product;
291  priorProductNext = &(product->next);
292  }
293  xData_freeElementList( smr, list );
294  return( status );
295 }
296 /*
297 ************************************************************
298 */
299 static int _tpia_product_getDepositionEnergy( statusMessageReporting *smr, xData_element *depositionEnergy, tpia_product *product ) {
300 
302 
303  xData_addToAccessed( smr, depositionEnergy, 1 );
304  for( data = xData_getFirstElement( depositionEnergy ); data != NULL; data = xData_getNextElement( data ) ) {
305  if( strcmp( data->name, "grouped" ) == 0 ) {
306  if( tpia_misc_get2d_xShared_yHistogram_data_Grouped( smr, data, &(product->depositionEnergyGrouped) ) ) return( 1 ); }
307  else {
308  tpia_misc_setMessageError_Element( smr, NULL, depositionEnergy, __FILE__, __LINE__, 1, "unsupported deposition energy type = %s", data->name );
309  return( 1 );
310  }
311  }
312  return( 0 );
313 }
314 /*
315 ************************************************************
316 */
317 static int _tpia_product_getMultiplicityFromElement( statusMessageReporting *smr, xData_element *data, tpia_product *product ) {
318 
319  tpia_multiplicity *multiplicity, *prior, *current;
320  const char *timeScale;
321  int isDelayedNeutrons;
322  double dTimeScale;
323 
324  if( tpia_multiplicity_getTimeScaleFromElement( smr, data, &timeScale, &isDelayedNeutrons, &dTimeScale ) ) return( 1 );
325  if( ( isDelayedNeutrons == 0 ) && ( product->multiplicityVsEnergy != NULL ) ) {
326  tpia_misc_setMessageError_Element( smr, NULL, data, __FILE__, __LINE__, 1, "extra product multiplicity data" );
327  return( 1 );
328  }
329  if( ( multiplicity = tpia_multiplicity_createGetFromElement( smr, data, product->channel->target->nGroups ) ) == NULL ) return( 1 );
330  if( isDelayedNeutrons == 0 ) {
331  product->multiplicityVsEnergy = multiplicity; }
332  else {
333  if( product->delayedNeutronMultiplicityVsEnergy == NULL ) {
334  product->delayedNeutronMultiplicityVsEnergy = multiplicity; }
335  else {
336  if( product->delayedNeutronMultiplicityVsEnergy->timeScale > multiplicity->timeScale ) {
337  multiplicity->next = product->delayedNeutronMultiplicityVsEnergy;
338  product->delayedNeutronMultiplicityVsEnergy = multiplicity; }
339  else {
340  for( current = product->delayedNeutronMultiplicityVsEnergy->next, prior = product->delayedNeutronMultiplicityVsEnergy; current != NULL;
341  current = current->next ) {
342  if( current->timeScale > multiplicity->timeScale ) {
343  multiplicity->next = current;
344  prior->next = multiplicity;
345  break;
346  }
347  prior = current;
348  }
349  if( current == NULL ) prior->next = multiplicity;
350  }
351  }
352  }
353  return( 0 );
354 }
355 /*
356 ************************************************************
357 */
358 //long tpia_product_dataRequired( statusMessageReporting *smr, tpia_product *product ) {
360 
361  return( product->b_dataRequired );
362 }
363 /*
364 ************************************************************
365 */
367 
368  return( tpia_decayChannel_getFirstProduct( &(product->decayChannel) ) );
369 }
370 /*
371 ************************************************************
372 */
373 //tpia_product *tpia_product_getProductByIndex( statusMessageReporting *smr, tpia_product *product, int index ) {
375 
376  int i = 0;
377  tpia_product *p;
378 
379  if( index < 0 ) return( NULL );
380  for( p = tpia_product_getFirstProduct( product ); ( p != NULL ) && ( i < index ); p = tpia_decayChannel_getNextProduct( p ), i++ ) ;
381  return( p );
382 }
383 /*
384 ************************************************************
385 */
386 //int tpia_product_doesDecay( statusMessageReporting *smr, tpia_product *product ) {
388 
389  return( product->decayChannel.products != NULL );
390 }
391 /*
392 ************************************************************
393 */
394 //int tpia_product_numberOfProducts( statusMessageReporting *smr, tpia_product *product ) {
396 
397  return( product->decayChannel.numberOfProducts );
398 }
399 /*
400 ************************************************************
401 */
402 //int tpia_product_isDataPresent( statusMessageReporting *smr, tpia_product *product, int b_data ) {
404 
405  return( product->b_dataPresent && b_data );
406 }
407 /*
408 ************************************************************
409 */
410 //int tpia_product_sampleMultiplicity( statusMessageReporting *smr, tpia_product *product, double e_in, double r ) {
411 int tpia_product_sampleMultiplicity( statusMessageReporting *, tpia_product *product, double e_in, double r ) {
412 
413  int i, multiplicity;
414  tpia_multiplicity *multiplicityVsEnergy = product->multiplicityVsEnergy;
415  double *p = multiplicityVsEnergy->pointwise, dMult;
416 
417  if( e_in <= p[0] ) {
418  dMult = p[1]; }
419  else if( e_in >= p[2 * ( multiplicityVsEnergy->numberOfPointwise - 1 )] ) {
420  dMult = p[2 * multiplicityVsEnergy->numberOfPointwise - 1]; }
421  else {
422  for( i = 0; i < multiplicityVsEnergy->numberOfPointwise - 1; i++, p += 2 ) if( e_in < p[2] ) break;
423  dMult = ( e_in - p[0] ) / ( p[2] - p[0] );
424  dMult = dMult * p[3] + ( 1. - dMult ) * p[1];
425  }
426  multiplicity = (int) dMult;
427  if( r < ( dMult - multiplicity ) ) multiplicity++;
428 
429  return( multiplicity );
430 }
431 
432 #if defined __cplusplus
433 }
434 #endif