ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MCGIDI_misc.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MCGIDI_misc.cc
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <string.h>
8 #include <cmath>
9 #include <ctype.h>
10 
11 #ifdef WIN32
12 #include <direct.h>
13 #else
14 #include <unistd.h>
15 #endif
16 
17 #include <ptwXY.h>
19 
20 #include "MCGIDI.h"
21 #include "MCGIDI_misc.h"
22 #include "MCGIDI_fromTOM.h"
23 
24 #if defined __cplusplus
25 namespace GIDI {
26 using namespace GIDI;
27 #endif
28 
29 struct ZSymbol {
30  int Z;
31  char const *symbol;
32 };
33 
34 static struct ZSymbol ZSymbols[] = { { 0, "n" }, { 1, "H" }, { 2, "He" }, { 3, "Li" }, { 4, "Be" }, { 5, "B" }, { 6, "C" },
35  { 7, "N" }, { 8, "O" }, { 9, "F" }, { 10, "Ne" }, { 11, "Na" }, { 12, "Mg" }, { 13, "Al" }, { 14, "Si" }, { 15, "P" },
36  { 16, "S" }, { 17, "Cl" }, { 18, "Ar" }, { 19, "K" }, { 20, "Ca" }, { 21, "Sc" }, { 22, "Ti" }, { 23, "V" }, { 24, "Cr" },
37  { 25, "Mn" }, { 26, "Fe" }, { 27, "Co" }, { 28, "Ni" }, { 29, "Cu" }, { 30, "Zn" }, { 31, "Ga" }, { 32, "Ge" }, { 33, "As" },
38  { 34, "Se" }, { 35, "Br" }, { 36, "Kr" }, { 37, "Rb" }, { 38, "Sr" }, { 39, "Y" }, { 40, "Zr" }, { 41, "Nb" }, { 42, "Mo" },
39  { 43, "Tc" }, { 44, "Ru" }, { 45, "Rh" }, { 46, "Pd" }, { 47, "Ag" }, { 48, "Cd" }, { 49, "In" }, { 50, "Sn" }, { 51, "Sb" },
40  { 52, "Te" }, { 53, "I" }, { 54, "Xe" }, { 55, "Cs" }, { 56, "Ba" }, { 57, "La" }, { 58, "Ce" }, { 59, "Pr" }, { 60, "Nd" },
41  { 61, "Pm" }, { 62, "Sm" }, { 63, "Eu" }, { 64, "Gd" }, { 65, "Tb" }, { 66, "Dy" }, { 67, "Ho" }, { 68, "Er" }, { 69, "Tm" },
42  { 70, "Yb" }, { 71, "Lu" }, { 72, "Hf" }, { 73, "Ta" }, { 74, "W" }, { 75, "Re" }, { 76, "Os" }, { 77, "Ir" }, { 78, "Pt" },
43  { 79, "Au" }, { 80, "Hg" }, { 81, "Tl" }, { 82, "Pb" }, { 83, "Bi" }, { 84, "Po" }, { 85, "At" }, { 86, "Rn" }, { 87, "Fr" },
44  { 88, "Ra" }, { 89, "Ac" }, { 90, "Th" }, { 91, "Pa" }, { 92, "U" }, { 93, "Np" }, { 94, "Pu" }, { 95, "Am" }, { 96, "Cm" },
45  { 97, "Bk" }, { 98, "Cf" }, { 99, "Es" }, { 100, "Fm" }, { 101, "Md" }, { 102, "No" }, { 103, "Lr" }, { 104, "Rf" }, { 105, "Db" },
46  { 106, "Sg" }, { 107, "Bh" }, { 108, "Hs" }, { 109, "Mt" } };
47 
48 static int MCGIDI_miscNameToZAm_getLevel( statusMessageReporting *smr, const char *name, const char *p );
50  int length, double *data, char const *fromUnits[2], char const *toUnits[2] );
51 /*
52 ************************************************************
53 */
55 
56  return( sizeof( ZSymbols ) / sizeof( struct ZSymbol ) );
57 }
58 /*
59 ************************************************************
60 */
61 const char *MCGIDI_misc_ZToSymbol( int iZ ) {
62 
63  if( ( iZ < 0 ) || ( iZ >= MCGIDI_misc_NumberOfZSymbols( ) ) ) return( NULL );
64  return( ZSymbols[iZ].symbol );
65 }
66 /*
67 ************************************************************
68 */
69 int MCGIDI_misc_symbolToZ( const char *Z ) {
70 
71  int i, n = MCGIDI_misc_NumberOfZSymbols( );
72 
73  for( i = 0; i < n; i++ ) {
74  if( strcmp( Z, ZSymbols[i].symbol ) == 0 ) return( ZSymbols[i].Z );
75  }
76  return( -1 );
77 }
78 /*
79 ************************************************************
80 */
81 int MCGIDI_miscNameToZAm( statusMessageReporting *smr, const char *name, int *Z, int *A, int *m, int *level ) {
82 
83  const char *p;
84  char s[1024] = "", *q, *e; /* Note 1) routine will fail when parts of a particle name can be longer than 1024. */
85 
86  if( strlen( name ) >= ( sizeof( s ) - 1 ) ) {
87  smr_setReportError2( smr, smr_unknownID, 0, "particle name too long: '%s'", name );
88  return( 1 );
89  }
90 
91  *Z = *A = *m = *level = 0;
92  if( ( !strncmp( "FissionProduct", name, 14 ) ) || !strncmp( "99120", name, 5 ) ) {
93  *Z = 99;
94  *A = 120;
95  return( 0 );
96  }
97  if( strcmp( "gamma", name ) == 0 ) return( 0 );
98  if( strcmp( "n", name ) == 0 ) { *A = 1; return( 0 ); }
99 
100  for( p = name, q = s; ( *p != 0 ) && !isdigit( *p ) && ( *p != '_' ); p++, q++ ) *q = *p; /* '_' only for "natural". */
101  if( *p == 0 ) {
102  smr_setReportError2( smr, smr_unknownID, 0, "unsupported particle name = '%s'", name );
103  return( 1 );
104  }
105  *q = 0;
106  if( ( *Z = MCGIDI_misc_symbolToZ( s ) ) < 0 ) {
107  smr_setReportError2( smr, smr_unknownID, 1, "Particle %s's symbol = '%s' not found", name, s ); }
108  else { /* Getting here implies that *p is a digit. */
109  if( *p == '_' ) {
110  if( strncmp( p, "_natural", 8 ) == 0 ) {
111  p += 8;
112  if( *p ) *level = MCGIDI_miscNameToZAm_getLevel( smr, name, p ); }
113  else {
114  smr_setReportError2( smr, smr_unknownID, 0, "expecting 'natural': %s", name );
115  } }
116  else {
117  for( q = s; isdigit( *p ); p++, q++ ) *q = *p;
118  *q = 0;
119  if( strcmp( s, "natural" ) == 0 ) {
120  e = s;
121  while( *e ) e++; /* Loop checking, 11.06.2015, T. Koi*/ }
122  else {
123  *A = (int) strtol( s, &e, 10 );
124  }
125  if( *e != 0 ) {
126  smr_setReportError2( smr, smr_unknownID, 1, "Failed to convert A to integer in particle name %s", name ); }
127  else { /* Getting here implies that *p == '_' or 0. */
128  if( *p ) *level = MCGIDI_miscNameToZAm_getLevel( smr, name, p );
129  }
130  }
131  }
132 
133  return( !smr_isOk( smr ) );
134 }
135 /*
136 ************************************************************
137 */
138 static int MCGIDI_miscNameToZAm_getLevel( statusMessageReporting *smr, const char *name, const char *p ) {
139 
140  int level = 0;
141  char *e;
142 
143  if( *p == '_' ) {
144  p++;
145  switch( *p ) {
146  case 'e' :
147  p++;
148  level = (int) strtol( p, &e, 10 );
149  if( *e != 0 ) smr_setReportError2( smr, smr_unknownID, 1, "Failed to convert level to integer in particle name %s", name );
150  break;
151  case 'c' :
153  break;
154  case 's' :
155  level = MCGIDI_particleLevel_sum;
156  break;
157  default :
158  smr_setReportError2( smr, smr_unknownID, 0, "invalid 'natural': %s", name );
159  } }
160  else {
161  smr_setReportError2( smr, smr_unknownID, 0, "invalid level specifier: %s", name );
162  }
163  return( level );
164 }
165 /*
166 ************************************************************
167 */
168 char const *MCGIDI_misc_pointerToTOMAttributeIfAllOk( statusMessageReporting *smr, const char *path, int required,
169  xDataTOM_attributionList *attributes, const char *name, const char *file, int line ) {
170 
171  char const *value;
172 
173  if( !smr_isOk( smr ) ) return( NULL );
174  if( ( value = xDataTOMAL_getAttributesValue( attributes, name ) ) == NULL ) {
175  if( required ) {
176  smr_setReportError( smr, NULL, file, line, __func__, smr_unknownID, 1, "element does not have attribute named %s for file = %d", name, path );
177  }
178  }
179  return( value );
180 }
181 /*
182 ************************************************************
183 */
184 char const *MCGIDI_misc_pointerToAttributeIfAllOk( statusMessageReporting *smr, xDataXML_element *element, const char *path, int required,
185  xDataTOM_attributionList *attributes, const char *name, const char *file, int line ) {
186 
187  char const *value;
188 
189  if( !smr_isOk( smr ) ) return( NULL );
190  if( ( value = xDataTOMAL_getAttributesValue( attributes, name ) ) == NULL ) {
191  if( required ) {
192  if( element != NULL ) {
193  MCGIDI_misc_setMessageError_Element( smr, NULL, element, file, line, 1, "element does not have attribute named %s", name ); }
194  else {
195  smr_setReportError( smr, NULL, file, line, __func__, smr_unknownID, 1, "element does not have attribute named %s for file = %d", name, path );
196  }
197  }
198  }
199  return( value );
200 }
201 /*
202 ************************************************************
203 */
204 int MCGIDI_misc_setMessageError_Element( statusMessageReporting *smr, void *userInterface, xDataXML_element *element, const char *file, int line, int code,
205  const char *fmt, ... ) {
206 
207  int status = 0;
208  va_list args;
209  char *msg;
210 
211  va_start( args, fmt );
212  msg = smr_vallocateFormatMessage( fmt, &args );
213  va_end( args );
214  if( msg == NULL ) {
215  status = 1;
216  va_start( args, fmt );
217  smr_vsetReportError( smr, userInterface, file, line, __func__, smr_unknownID, code, fmt, &args );
218  va_end( args ); }
219  else {
220  status = smr_setReportError( smr, userInterface, file, line, __func__, smr_unknownID, code, "%s for element %s", msg, element->name );
221  smr_freeMemory( (void **) &msg );
222  }
223  return( status );
224 }
225 /*
226 ************************************************************
227 */
229 /*
230 * Returns -2 is d < first point of ds, -1 if > last point of ds and the lower index of ds otherwise.
231 */
232  xDataTOM_Int imin = 0, imid, imax = n - 1;
233 
234  if( d < ds[0] ) return( -2 );
235  if( d > ds[n-1] ) return( -1 );
236  while( 1 ) { // Loop checking, 11.06.2015, T. Koi
237  imid = ( imin + imax ) >> 1;
238  if( imid == imin ) break;
239  if( d < ds[imid] ) {
240  imax = imid; }
241  else {
242  imin = imid;
243  }
244  }
245  return( imin );
246 }
247 /*
248 ************************************************************
249 */
250 char *MCGIDI_misc_getAbsPath( statusMessageReporting *smr, const char *fileName ) {
251 /*
252 * User must free returned string.
253 */
254  int n = (int) strlen( fileName ) + 1, nCwd = 0;
255  char *absPath, cwd[4 * 1024] = "", *p, *needle;
256 
257  if( fileName[0] != '/' ) {
258  //if( getcwd( cwd, sizeof( cwd ) + 1 ) == NULL ) {
259  //TK modified above line for compiler(gcc.4.8) warning message
260  if( getcwd( cwd, sizeof( cwd ) ) == NULL ) {
261  smr_setReportError2p( smr, smr_unknownID, -1, "hardwired cwd too small" );
262  return( NULL );
263  }
264  nCwd = (int) strlen( cwd );
265  n += nCwd + 1; /* cwd + '/'. */
266  }
267  if( ( absPath = (char *) smr_malloc2( smr, n, 0, "absPath" ) ) == NULL ) return( NULL );
268  if( fileName[0] != '/' ) {
269  strcpy( absPath, cwd );
270  strcat( absPath, "/" );
271  strcat( absPath, fileName ); }
272  else {
273  strcpy( absPath, fileName );
274  }
275 
276  while( 1 ) { /* Remove all ./ from path. */ // Loop checking, 11.06.2015, T. Koi
277  if( ( needle = strstr( absPath, "/./" ) ) == NULL ) break;
278  p = needle;
279  for( needle += 2; *needle; p++, needle++ ) *p = *needle;
280  *p = 0;
281  }
282  while( 1 ) { /* Remove all ../ from path. */ // Loop checking, 11.06.2015, T. Koi
283  if( ( needle = strstr( absPath, "/../" ) ) == NULL ) break;
284  p = needle - 1;
285  while( ( p > absPath ) && ( *p != '/' ) ) p--; // Loop checking, 11.06.2015, T. Koi
286  if( *p != '/' ) break; /* This should not happen if path is legit, I think, and I do not know what to do so will leave it. */
287  if( p == absPath ) break; /* Ditto. */
288  for( needle += 3; *needle; p++, needle++ ) *p = *needle;
289  *p = 0;
290  }
291  return( absPath );
292 }
293 /*
294 ************************************************************
295 */
297 
298  int i;
299  xDataXML_attribute *attribute;
300 
301  xDataTOMAL_initial( smr, TOM );
302  for( i = 0; ; i++ ) {
303  if( ( attribute = xDataXML_attributeByIndex( XML, i ) ) == NULL ) break;
304  if( xDataTOMAL_addAttribute( smr, TOM, attribute->name, attribute->value ) != 0 ) goto err;
305  }
306  return( 0 );
307 
308 err:
309  xDataTOMAL_release( TOM );
310  return( 1 );
311 }
312 /*
313 ************************************************************
314 */
316 
317  char const *frameString;
319 
320  if( ( frameString = xDataTOM_getAttributesValueInElement( frameElement, MCGIDI_token_productFrame ) ) != NULL ) {
321  if( ( frame = xDataTOM_axis_stringToFrame( smr, frameString ) ) == xDataTOM_frame_invalid ) {
322  smr_setReportError2( smr, smr_unknownID, 1, "Invalid frame = '%s'", frameString );
323  }
324  }
325  return( frame );
326 }
327 /*
328 ************************************************************
329 */
330 int MCGIDI_misc_PQUStringToDouble( statusMessageReporting *smr, char const *str, char const *unit, double conversion, double *value ) {
331 /*
332 * Currently, white spaces are not allowed after the unit.
333 *
334 * Examples of allowed strings are: "2.39e6 eV", " 2.39e6eV" and " 2.39e6 eV".
335 */
336  char const *s = str;
337  char *e;
338 
339 
340  while( isspace( *s ) ) s++; // Loop checking, 11.06.2015, T. Koi
341  *value = strtod( s, &e ) * conversion;
342  if( e == s ) {
343  smr_setReportError2( smr, smr_unknownID, 1, "no number at start of string = <%s>", str );
344  return( 1 );
345  }
346  while( isspace( *e ) ) e++; // Loop checking, 11.06.2015, T. Koi
347  if( strcmp( e, unit ) != 0 ) {
348  smr_setReportError2( smr, smr_unknownID, 1, "unit = '%s' not '%s' in '%s'", e, unit, str );
349  return( 1 );
350  }
351  return( 0 );
352 }
353 /*
354 ************************************************************
355 */
356 int MCGIDI_misc_PQUStringToDoubleInUnitOf( statusMessageReporting *smr, char const *str, char const *toUnit, double *value ) {
357 /*
358 * Currently, white spaces are not allowed after the unit.
359 *
360 * Examples of allowed strings are: "2.39e6 eV", " 2.39e6eV" and " 2.39e6 eV".
361 */
362  char const *s1 = str;
363  char *e1;
364  double factor;
365 
366  while( isspace( *s1 ) ) s1++; // Loop checking, 11.06.2015, T. Koi
367  *value = strtod( s1, &e1 );
368  if( e1 == s1 ) {
369  smr_setReportError2( smr, smr_unknownID, 1, "no number at start of string = <%s>", str );
370  return( 1 );
371  }
372  while( isspace( *e1 ) ) e1++; // Loop checking, 11.06.2015, T. Koi
373 
374  factor = MCGIDI_misc_getUnitConversionFactor( smr, e1, toUnit );
375  *value *= factor;
376  return( !smr_isOk( smr ) );
377 }
378 /*
379 ************************************************************
380 */
381 double MCGIDI_misc_getUnitConversionFactor( statusMessageReporting *smr, char const *fromUnit, char const *toUnit ) {
382 /*
383 * This is a kludge until units are better supported.
384 */
385  if( strcmp( fromUnit, toUnit ) == 0 ) return( 1.0 );
386 
387  if( strcmp( fromUnit, "eV" ) == 0 ) {
388  if( strcmp( toUnit, "MeV" ) == 0 ) return( 1e-6 ); }
389  else if( strcmp( fromUnit, "MeV" ) == 0 ) {
390  if( strcmp( toUnit, "eV" ) == 0 ) return( 1e+6 ); }
391  else if( strcmp( fromUnit, "1/eV" ) == 0 ) {
392  if( strcmp( toUnit, "1/MeV" ) == 0 ) return( 1e+6 ); }
393  else if( strcmp( fromUnit, "1/MeV" ) == 0 ) {
394  if( strcmp( toUnit, "1/eV" ) == 0 ) return( 1e-6 ); }
395  else if( strcmp( fromUnit, "K" ) == 0 ) {
396  if( strcmp( toUnit, "MeV/k" ) == 0 ) return( 8.617343183775137e-11 );
397  }
398 
399  smr_setReportError2( smr, smr_unknownID, 1, "Cannot convert unit '%s' to unit '%s'", fromUnit, toUnit );
400  return( 1.0 );
401 }
402 /*
403 ************************************************************
404 */
406  ptwXY_interpolation interpolation, char const *toUnits[2] ) {
407 
408  int length;
409  double *data;
410  char const *fromUnits[2];
411 
412  fromUnits[0] = xDataTOM_subAxes_getUnit( smr, &(XYs->subAxes), 0 );
413  if( !smr_isOk( smr ) ) return( NULL );
414  fromUnits[1] = xDataTOM_subAxes_getUnit( smr, &(XYs->subAxes), 1 );
415  if( !smr_isOk( smr ) ) return( NULL );
416 
417  length = xDataTOM_XYs_getData( XYs, &data );
418 
419  return( MCGIDI_misc_Data2ptwXYPointsInUnitsOf( smr, interpolation, length, data, fromUnits, toUnits ) );
420 }
421 /*
422 ************************************************************
423 */
425 
426  int length;
427  double *data;
428  xDataTOM_axes *axes = &(linear->xDataInfo.axes);
429  char const *fromUnits[2];
430  ptwXY_interpolation interpolation;
431 
432  if( axes->numberOfAxes != 2 ) {
433  smr_setReportError2( smr, smr_unknownID, 1, "axes must have 2 axis, it has %d", axes->numberOfAxes );
434  return( NULL );
435  }
436 
437  if( MCGIDI_fromTOM_interpolation( smr, linear, 0, &interpolation ) != 0 ) return( NULL );
438  fromUnits[0] = axes->axis[0].unit;
439  fromUnits[1] = axes->axis[1].unit;
440 
441  length = xDataTOM_XYs_getDataFromXDataInfo( (xDataTOM_xDataInfo *) &(linear->xDataInfo), &data );
442  return( MCGIDI_misc_Data2ptwXYPointsInUnitsOf( smr, interpolation, length, data, fromUnits, toUnits ) );
443 }
444 /*
445 ************************************************************
446 */
448  int length, double *data, char const *fromUnits[2], char const *toUnits[2] ) {
449 
450  double xFactor, yFactor;
451  ptwXYPoints *ptwXY = NULL;
452  nfu_status status;
453 
454  xFactor = MCGIDI_misc_getUnitConversionFactor( smr, fromUnits[0], toUnits[0] );
455  if( !smr_isOk( smr ) ) goto err;
456  yFactor = MCGIDI_misc_getUnitConversionFactor( smr, fromUnits[1], toUnits[1] );
457  if( !smr_isOk( smr ) ) goto err;
458 
459 
460  ptwXY = ptwXY_create( interpolation, NULL, 2., 1e-3, length, 10, length, data, &status, 0 );
461  if( status != nfu_Okay ) {
462  smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_create err = %d: %s\n", status, nfu_statusMessage( status ) );
463  goto err;
464  }
465 
466  if( ( xFactor != 1. ) || ( yFactor != 1. ) ) {
467  if( ( status = ptwXY_scaleOffsetXAndY( ptwXY, xFactor, 0., yFactor, 0. ) ) != nfu_Okay ) {
468  smr_setReportError2( smr, smr_unknownID, 1, "ptwXY_scaleOffsetXAndY err = %d: %s\n", status, nfu_statusMessage( status ) );
469  goto err;
470  }
471  }
472 
473  return( ptwXY );
474 
475 err:
476  if( ptwXY != NULL ) ptwXY_free( ptwXY );
477  return( NULL );
478 }
479 /*
480 ************************************************************
481 */
482 void MCGIDI_misc_updateTransportabilitiesMap( transportabilitiesMap *transportabilities, int PoPID, enum MCGIDI_transportability transportability ) {
483 
484  transportabilitiesMap::iterator iter = transportabilities->find( PoPID );
485 
486  if( iter != transportabilities->end( ) ) {
487  switch ( iter->second ) {
489  break;
491  switch( transportability ) {
494  transportability = MCGIDI_transportability_none;
495  break;
497  break;
499  transportability = MCGIDI_transportability_partial;
500  break;
501  }
502  break;
504  transportability = MCGIDI_transportability_partial;
505  break;
507  switch( transportability ) {
510  transportability = MCGIDI_transportability_partial;
511  break;
514  break;
515  }
516  break;
517  }
518  }
519  (*transportabilities)[PoPID] = transportability;
520 }
521 /*
522 ************************************************************
523 */
524 void MCGIDI_misc_updateTransportabilitiesMap2( transportabilitiesMap *transportabilities, int PoPID, int transportable ) {
525 
527 }
528 
529 #if defined __cplusplus
530 }
531 #endif
532