| Geant4 Cross Reference |
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>
18 #include <xDataTOM_importXML_private.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 );
49 static ptwXYPoints *MCGIDI_misc_Data2ptwXYPointsInUnitsOf( statusMessageReporting *smr, ptwXY_interpolation interpolation,
50 int length, double *data, char const *fromUnits[2], char const *toUnits[2] );
51 /*
52 ************************************************************
53 */
54 int MCGIDI_misc_NumberOfZSymbols( void ) {
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' :
152 level = MCGIDI_particleLevel_continuum;
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 */
228 xDataTOM_Int MCGIDI_misc_binarySearch( xDataTOM_Int n, double *ds, double d ) {
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 */
296 int MCGIDI_misc_copyXMLAttributesToTOM( statusMessageReporting *smr, xDataTOM_attributionList *TOM, xDataXML_attributionList *XML ) {
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 */
315 enum xDataTOM_frame MCGIDI_misc_getProductFrame( statusMessageReporting *smr, xDataTOM_element *frameElement ) {
316
317 char const *frameString;
318 enum xDataTOM_frame frame = xDataTOM_frame_invalid;
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 */
405 ptwXYPoints *MCGIDI_misc_dataFromXYs2ptwXYPointsInUnitsOf( statusMessageReporting *smr, xDataTOM_XYs *XYs,
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 */
424 ptwXYPoints *MCGIDI_misc_dataFromElement2ptwXYPointsInUnitsOf( statusMessageReporting *smr, xDataTOM_element *linear, char const *toUnits[2] ) {
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 */
447 static ptwXYPoints *MCGIDI_misc_Data2ptwXYPointsInUnitsOf( statusMessageReporting *smr, ptwXY_interpolation interpolation,
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 ) {
488 case MCGIDI_transportability_unknown :
489 break;
490 case MCGIDI_transportability_none :
491 switch( transportability ) {
492 case MCGIDI_transportability_unknown :
493 case MCGIDI_transportability_none :
494 transportability = MCGIDI_transportability_none;
495 break;
496 case MCGIDI_transportability_partial :
497 break;
498 case MCGIDI_transportability_full :
499 transportability = MCGIDI_transportability_partial;
500 break;
501 }
502 break;
503 case MCGIDI_transportability_partial :
504 transportability = MCGIDI_transportability_partial;
505 break;
506 case MCGIDI_transportability_full :
507 switch( transportability ) {
508 case MCGIDI_transportability_none :
509 case MCGIDI_transportability_partial :
510 transportability = MCGIDI_transportability_partial;
511 break;
512 case MCGIDI_transportability_unknown :
513 case MCGIDI_transportability_full :
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
526 MCGIDI_misc_updateTransportabilitiesMap( transportabilities, PoPID, ( transportable ? MCGIDI_transportability_full : MCGIDI_transportability_none ) );
527 }
528
529 #if defined __cplusplus
530 }
531 #endif
532
533