Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
tpia_misc.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 
37 #include "Randomize.hh"
38 
39 #include <stdio.h>
40 #include <stdlib.h>
41 #include <string.h>
42 #include <cmath>
43 #include <tpia_target.h>
44 #include <tpia_misc.h>
45 
46 #include <string>
47 
48 #if defined __cplusplus
49 namespace GIDI {
50 using namespace GIDI;
51 #endif
52 
53 struct ZSymbol {
54  int Z;
55  //char *symbol;
56  std::string symbol;
57 };
58 
59 static struct ZSymbol ZSymbols[] = { { 0, "n" }, { 1, "H" }, { 2, "He" }, { 3, "Li" }, { 4, "Be" }, { 5, "B" }, { 6, "C" },
60  { 7, "N" }, { 8, "O" }, { 9, "F" }, { 10, "Ne" }, { 11, "Na" }, { 12, "Mg" }, { 13, "Al" }, { 14, "Si" }, { 15, "P" },
61  { 16, "S" }, { 17, "Cl" }, { 18, "Ar" }, { 19, "K" }, { 20, "Ca" }, { 21, "Sc" }, { 22, "Ti" }, { 23, "V" }, { 24, "Cr" },
62  { 25, "Mn" }, { 26, "Fe" }, { 27, "Co" }, { 28, "Ni" }, { 29, "Cu" }, { 30, "Zn" }, { 31, "Ga" }, { 32, "Ge" }, { 33, "As" },
63  { 34, "Se" }, { 35, "Br" }, { 36, "Kr" }, { 37, "Rb" }, { 38, "Sr" }, { 39, "Y" }, { 40, "Zr" }, { 41, "Nb" }, { 42, "Mo" },
64  { 43, "Tc" }, { 44, "Ru" }, { 45, "Rh" }, { 46, "Pd" }, { 47, "Ag" }, { 48, "Cd" }, { 49, "In" }, { 50, "Sn" }, { 51, "Sb" },
65  { 52, "Te" }, { 53, "I" }, { 54, "Xe" }, { 55, "Cs" }, { 56, "Ba" }, { 57, "La" }, { 58, "Ce" }, { 59, "Pr" }, { 60, "Nd" },
66  { 61, "Pm" }, { 62, "Sm" }, { 63, "Eu" }, { 64, "Gd" }, { 65, "Tb" }, { 66, "Dy" }, { 67, "Ho" }, { 68, "Er" }, { 69, "Tm" },
67  { 70, "Yb" }, { 71, "Lu" }, { 72, "Hf" }, { 73, "Ta" }, { 74, "W" }, { 75, "Re" }, { 76, "Os" }, { 77, "Ir" }, { 78, "Pt" },
68  { 79, "Au" }, { 80, "Hg" }, { 81, "Tl" }, { 82, "Pb" }, { 83, "Bi" }, { 84, "Po" }, { 85, "At" }, { 86, "Rn" }, { 87, "Fr" },
69  { 88, "Ra" }, { 89, "Ac" }, { 90, "Th" }, { 91, "Pa" }, { 92, "U" }, { 93, "Np" }, { 94, "Pu" }, { 95, "Am" }, { 96, "Cm" },
70  { 97, "Bk" }, { 98, "Cf" }, { 99, "Es" }, { 100, "Fm" }, { 101, "Md" }, { 102, "No" }, { 103, "Lr" }, { 104, "Rf" }, { 105, "Db" },
71  { 106, "Sg" }, { 107, "Bh" }, { 108, "Hs" }, { 109, "Mt" } };
72 
73 /*
74 ************************************************************
75 */
77 
78  return( sizeof( ZSymbols ) / sizeof( struct ZSymbol ) );
79 }
80 /*
81 ************************************************************
82 */
83 const char *tpia_misc_ZToSymbol( int iZ ) {
84 
85  if( ( iZ < 0 ) || ( iZ >= tpia_misc_NumberOfZSymbols( ) ) ) return( NULL );
86  //return( ZSymbols[iZ].symbol );
87  return( ZSymbols[iZ].symbol.c_str() );
88 }
89 /*
90 ************************************************************
91 */
92 int tpia_misc_symbolToZ( const char *Z ) {
93 
94  int i, n = tpia_misc_NumberOfZSymbols( );
95 
96  for( i = 0; i < n; i++ ) {
97  //if( !strcmp( Z, ZSymbols[i].symbol ) ) return( ZSymbols[i].Z );
98  if( !strcmp( Z, ZSymbols[i].symbol.c_str() ) ) return( ZSymbols[i].Z );
99  }
100  return( -1 );
101 }
102 /*
103 ************************************************************
104 */
105 int tpia_miscNameToZAm( statusMessageReporting *smr, const char *name, int *Z, int *A, int *m ) {
106 
107  int i, n;
108  const char *p;
109  char s[1024] = "", *q, *e; /* Note 1) routine will fail when parts of a particle name can be longer than 1024. */
110 
111  n = sizeof( s ) - 1;
112 
113  *Z = *A = *m = 0;
114  if( !strncmp( "FissionProduct", name, 14 ) ) {
115  *Z = 99;
116  *A = 120;
117  return( 0 );
118  }
119  if( !strcmp( "gamma", name ) ) return( 0 );
120  for( p = name, q = s, i = 0; ( *p != '_' ) && ( i != n ); p++, q++, i++ ) *q = *p;
121  if( i == n ) { /* See note 1 above. */
122  smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Failed to find first '_' in particle name %s", name ); }
123  else {
124  *q = 0;
125  if( ( *Z = tpia_misc_symbolToZ( s ) ) < 0 ) {
126  smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Particle %s's symbol = '%s' not found", name, s ); }
127  else { /* Getting here implies that *p == '_'. */
128  for( p++, q = s; ( *p != '_' ) && ( *p != 0 ) && ( i != n ); p++, q++, i++ ) *q = *p;
129  if( i == n ) { /* See note 1 above. */
130  smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Failed to find second '_' in particle name %s", name ); }
131  else {
132  *q = 0;
133  if( strcmp( s, "natural" ) == 0 ) {
134  e = s;
135  while( *e ) e++; }
136  else {
137  *A = (int) strtol( s, &e, 10 );
138  }
139  if( *e != 0 ) {
140  smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Failed to convert A to integer in particle name %s", name ); }
141  else { /* Getting here implies that *p == '_' or 0. */
142  *m = 0;
143  if( *p == '_' ) {
144  p++;
145  if( *p != 'm' ) {
146  smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Particle name %s missing meta-stable label 'm'", name ); }
147  else {
148  p++;
149  *m = (int) strtol( p, &e, 10 );
150  if( *e != 0 ) smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Failed to convert m to integer in particle name %s", name );
151  }
152  }
153  }
154  }
155  }
156  }
157 
158  return( !smr_isOk( smr ) );
159 }
160 /*
161 ************************************************************
162 */
163 char *tpia_misc_pointerToAttributeIfAllOk( statusMessageReporting *smr, xData_element *element, const char *path, int required,
164  xData_attributionList *attributes, const char *name, const char *file, int line ) {
165 
166  char *value;
167 
168  if( !smr_isOk( smr ) ) return( NULL );
169  if( ( value = xData_getAttributesValue( attributes, name ) ) == NULL ) {
170  if( required ) {
171  if( element != NULL ) {
172  tpia_misc_setMessageError_Element( smr, NULL, element, file, line, 1, "element does not have attribute named %s", name ); }
173  else {
174  smr_setMessageError( smr, NULL, file, line, 1, "element does not have attribute named %s for file = %d", name, path );
175  }
176  }
177  }
178  return( value );
179 }
180 /*
181 ************************************************************
182 */
183 int tpia_misc_setMessageError_Element( statusMessageReporting *smr, void *userInterface, xData_element *element, const char *file, int line, int code,
184  const char *fmt, ... ) {
185 
186  int status = 0;
187  va_list args;
188  char *msg;
189 
190  va_start( args, fmt );
191  msg = smr_vallocateFormatMessage( fmt, &args );
192  va_end( args );
193  if( msg == NULL ) {
194  status = 1;
195  va_start( args, fmt );
196  smr_vsetMessageError( smr, userInterface, file, line, code, fmt, &args );
197  va_end( args ); }
198  else {
199  status = smr_setMessageError( smr, userInterface, file, line, code, "%s for element %s at line %d column %d", msg, element->fullName,
200  (int) element->docInfo.line, (int) element->docInfo.column );
201  free( msg );
202  }
203  return( status );
204 }
205 /*
206 ************************************************************
207 */
208 xData_Int tpia_misc_binarySearch( xData_Int n, double *ds, double d ) {
209 
210  xData_Int imin = 0, imid, imax = n - 1;
211 
212  if( d < ds[0] ) return( -2 );
213  if( d > ds[n-1] ) return( -1 );
214  while( 1 ) {
215  imid = ( imin + imax ) >> 1;
216  if( imid == imin ) break;
217  if( d < ds[imid] ) {
218  imax = imid; }
219  else {
220  imin = imid;
221  }
222  }
223  return( imin );
224 }
225 /*
226 ************************************************************
227 */
229 
230  xData_element *xDataElement;
231  double *d = NULL;
232  xData_Int length_;
233 
234  xData_addToAccessed( smr, element, 1 );
235  //if( ( xDataElement = xData_getOneElementByTagName( smr, element, "xData", 1 ) ) != NULL ) {
236  if( ( xDataElement = xData_getOneElementByTagName( smr, element, (char*)"xData", 1 ) ) != NULL ) {
237  xData_addToAccessed( smr, xDataElement, 1 );
238  if( xData_is_2d_xy( smr, &(xDataElement->xDataTypeInfo), 1 ) ) {
239  d = xData_2d_xy_allocateCopyData( smr, xDataElement, &length_ );
240  *length = (xData_Int) length_;
241  }
242  }
243  return( d );
244 }
245 /*
246 ************************************************************
247 */
248 double *tpia_misc_get2dxindex_y_data( statusMessageReporting *smr, xData_element *element, xData_Int *start, xData_Int *end, double *xValues ) {
249 
250  xData_element *xDataElement;
251  double *d = NULL;
252 
253  xData_addToAccessed( smr, element, 1 );
254  //if( ( xDataElement = xData_getOneElementByTagName( smr, element, "xData", 1 ) ) != NULL ) {
255  if( ( xDataElement = xData_getOneElementByTagName( smr, element, (char*) "xData", 1 ) ) != NULL ) {
256  xData_addToAccessed( smr, xDataElement, 1 );
257  if( xData_is_2d_xindex_y( smr, &(xDataElement->xDataTypeInfo), 1 ) ) {
258  if( start != NULL ) *start = xDataElement->xDataTypeInfo.start;
259  if( end != NULL ) *end = xDataElement->xDataTypeInfo.end;
260  d = xData_2d_xindex_y_toFilledYs( smr, xDataElement, xValues );
261  }
262  }
263  return( d );
264 }
265 /*
266 ************************************************************
267 */
269 
270  xData_Int i;
271  xData_element *xDataElement;
272  double *d = NULL;
273 
274  xData_addToAccessed( smr, element, 1 );
275  //if( ( xDataElement = xData_getOneElementByTagName( smr, element, "xData", 1 ) ) != NULL ) {
276  if( ( xDataElement = xData_getOneElementByTagName( smr, element, (char*) "xData", 1 ) ) != NULL ) {
277  xData_addToAccessed( smr, xDataElement, 1 );
278  if( ( d = xData_2d_xShared_yHistogram_copyData( smr, xDataElement, &i ) ) != NULL ) {
279  if( start != NULL ) *start = xDataElement->xDataTypeInfo.start;
280  if( end != NULL ) *end = xDataElement->xDataTypeInfo.end;
281  if( length != NULL ) *length = xDataElement->xDataTypeInfo.length;
282  }
283  }
284  return( d );
285 }
286 /*
287 ************************************************************
288 */
290 
291  if( ( group->data = tpia_misc_get2d_xShared_yHistogram_data( smr, element, &(group->start), &(group->end), &(group->length) ) ) == NULL ) return( 1 );
292  return( 0 );
293 }
294 /*
295 ************************************************************
296 */
297 //double tpia_misc_getPointwiseCrossSectionAtE( statusMessageReporting *smr, tpia_1dData *crossSection, double *energyGrid, xData_Int index, double e_in ) {
298 double tpia_misc_getPointwiseCrossSectionAtE( statusMessageReporting *, tpia_1dData *crossSection, double *energyGrid, xData_Int index, double e_in ) {
299 
300  double xsec = 0.0, e1, e2;
301 
302  if( ( index >= crossSection->start ) && ( index < crossSection->end ) ) {
303  e1 = energyGrid[index];
304  e2 = energyGrid[index + 1];
305  index -= crossSection->start;
306  if( e1 == e2 ) { /* Allow for future where step function may be allowed. */
307  xsec = 0.5 * ( crossSection->data[index] + crossSection->data[index + 1] ); }
308  else {
309  xsec = ( crossSection->data[index] * ( e2 - e_in ) + crossSection->data[index + 1] * ( e_in - e1 ) ) / ( e2 - e1 );
310  }
311  }
312  return( xsec );
313 }
314 /*
315 ************************************************************
316 */
318 
319  xData_element *element;
320 
321  xData_addToAccessed( smr, parent, 1 );
322  //if( ( element = xData_getOneElementByTagName( smr, parent, "equalProbableBins", 0 ) ) == NULL ) return( NULL );
323  if( ( element = xData_getOneElementByTagName( smr, parent, (char*) "equalProbableBins", 0 ) ) == NULL ) return( NULL );
324  if( xData_convertAttributeTo_xData_Int( smr, element, "nBins", nBins ) != 0 ) {
325  tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "missing or invalid nBins attribute" );
326  return( NULL );
327  }
328  return( tpia_misc_getEqualProbableBins( smr, element, "energy", *nBins, n ) );
329 }
330 /*
331 ************************************************************
332 */
334  xData_Int *n ) {
335 
336 
337  int i, j;
338  xData_Int index, size;
339  xData_elementList *list;
340  xData_element *element, *xData;
341  double *d;
342  tpia_EqualProbableBinSpectrum *epbs = NULL, *epb;
343 
344  xData_addToAccessed( smr, parent, 1 );
345  list = xData_getElementsByTagNameAndSort( smr, parent, name, NULL, NULL );
346  if( list->n == 0 ) {
347  tpia_misc_setMessageError_Element( smr, NULL, parent, __FILE__, __LINE__, 1, "bins does not contain any %s elements", name ); }
348  else {
349  *n = list->n;
350  size = list->n * ( sizeof( tpia_EqualProbableBinSpectrum ) + ( nBins + 1 ) * sizeof( double ) );
351  //if( ( epbs = xData_malloc2( smr, size, 0, "energies" ) ) != NULL ) {
352  if( ( epbs = (tpia_EqualProbableBinSpectrum*) xData_malloc2( smr, size, 0, "energies" ) ) != NULL ) {
353  d = (double *) &(epbs[list->n]);
354  for( i = 0, epb = epbs; i < list->n; i++, epb++ ) { /* Loop to test nBins and index are proper. */
355  element = list->items[i].element;
356  xData_addToAccessed( smr, element, 1 );
357  if( xData_convertAttributeTo_xData_Int( smr, element, "index", &index ) != 0 ) {
358  tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "missing or invalid index attribute" );
359  //epbs = xData_free( smr, epbs );
360  epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
361  break;
362  }
363  if( index != i ) {
364  tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "index = %lld is not incremental", index );
365  //epbs = xData_free( smr, epbs );
366  epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
367  break;
368  }
369  if( ( j = xData_convertAttributeToDouble( smr, element, "value", &(epb->value) ) ) != 0 ) {
370  if( j == 1 ) {
371  tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "element does not have value attribute" ); }
372  else {
373  tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "failed to convert value attribute to double" );
374  }
375  //epbs = xData_free( smr, epbs );
376  epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
377  break;
378  }
379  if( ( xData = xData_getElements_xDataElement( smr, element ) ) == NULL ) {
380  //epbs = xData_free( smr, epbs );
381  epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
382  break;
383  }
384  xData_addToAccessed( smr, xData, 1 );
385  epb->index = index;
386  epb->nBins = nBins;
387  epb->bins = d;
388  if( xData_1d_x_copyData( smr, xData, ( nBins + 1 ) * sizeof( double ), d ) != 0 ) {
389  //epbs = xData_free( smr, epbs );
390  epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
391  break;
392  }
393  d += nBins + 1;
394  }
395  }
396  }
397  xData_freeElementList( smr, list );
398  return( epbs );
399 }
400 /*
401 ************************************************************
402 */
403 double tpia_misc_drng( double (*rng)( void * ), void *rngState ) {
404 
405  double r;
406 
407  if( rng != NULL ) {
408  r = rng( rngState ); }
409  else {
410  //r = drand48( );
412 
413  }
414  return( r );
415 }
416 /*
417 ************************************************************
418 */
419 //int tpia_misc_sampleEqualProbableBin( statusMessageReporting *smr, tpia_decaySamplingInfo *decaySamplingInfo, double e_in, int nBins,
420 int tpia_misc_sampleEqualProbableBin( statusMessageReporting *, tpia_decaySamplingInfo *decaySamplingInfo, double e_in, int nBins,
421  tpia_EqualProbableBinSpectra *binned, double *value_ ) {
422 
423  int i, j, index1, index2, method = 0;
424  double fE = 1., r, value1, value2, value3, P12, P23, value = -2.;
425 
426  if( e_in <= binned->energies[0].value ) {
427  index1 = 0;
428  index2 = 0; }
429  else if( e_in >= binned->energies[binned->numberOfEs-1].value ) {
430  index1 = binned->numberOfEs - 1;
431  index2 = binned->numberOfEs - 1; }
432  else {
433  for( i = 0; i < binned->numberOfEs - 1; i++ ) {
434  if( e_in <= binned->energies[i].value ) break;
435  }
436  i--;
437  index1 = i;
438  index2 = i;
439  if( e_in != binned->energies[i].value ) {
440  index2++;
441  fE = ( e_in - binned->energies[i].value ) / ( binned->energies[i+1].value - binned->energies[i].value );
442  }
443  }
444  r = tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState );
445  j = (int) (r * nBins);
446  if( j >= nBins ) j = nBins - 1;
447  if( j < 0 ) j = 0;
448  r = tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState ); // Do not change r until after Point1 below.
450  method = 1;
451  if( ( ( j == 0 ) && ( r <= 0.5 ) ) || ( j == ( nBins - 1 ) && r > 0.5 ) ) method = 0;
452  }
453  if( method == 0 ) { /* Constant interpolaton. */
454  value1 = ( 1. - fE ) * binned->energies[index1].bins[j] + fE * binned->energies[index2].bins[j];
455  value2 = ( 1. - fE ) * binned->energies[index1].bins[j+1] + fE * binned->energies[index2].bins[j+1];
456  value = ( 1. - r ) * value1 + r * value2; }
457  else { /* Linear interpolation. */
458  if( r <= 0.5 ) j--; /* Point1: Above test insures that j is not 0 (nBins-1) if r <= 0.5 (> 0.5); */
459  value1 = ( 1. - fE ) * binned->energies[index1].bins[j] + fE * binned->energies[index2].bins[j];
460  value2 = ( 1. - fE ) * binned->energies[index1].bins[j+1] + fE * binned->energies[index2].bins[j+1];
461  value3 = ( 1. - fE ) * binned->energies[index1].bins[j+2] + fE * binned->energies[index2].bins[j+2];
462  P12 = 1. / ( value2 - value1 );
463  P23 = 1. / ( value3 - value2 );
464  r = tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState );
465  if( 0.25 * ( 1.0 + 2.0 * ( value2 - value1 ) / ( value3 - value1 ) ) > r ) {
466  P23 = 2. / ( value3 - value1 );
467  value3 = value2; }
468  else {
469  P12 = 2. / ( value3 - value1 );
470  value1 = value2;
471  }
472  r = tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState );
473  if( P23 != P12 ) r = ( -P12 + std::sqrt( P12 * P12 * ( 1. - r ) + r * P23 * P23 ) ) / ( P23 - P12 );
474  value = 0.5 * ( value1 + value2 + r * ( value3 - value1 ) );
475  }
476  *value_ = value;
477  return( 0 );
478 }
479 
480 #if defined __cplusplus
481 }
482 #endif