00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include "Randomize.hh"
00038
00039 #include <stdio.h>
00040 #include <stdlib.h>
00041 #include <string.h>
00042 #include <cmath>
00043 #include <tpia_target.h>
00044 #include <tpia_misc.h>
00045
00046 #include <string>
00047
00048 #if defined __cplusplus
00049 namespace GIDI {
00050 using namespace GIDI;
00051 #endif
00052
00053 struct ZSymbol {
00054 int Z;
00055
00056 std::string symbol;
00057 };
00058
00059 static struct ZSymbol ZSymbols[] = { { 0, "n" }, { 1, "H" }, { 2, "He" }, { 3, "Li" }, { 4, "Be" }, { 5, "B" }, { 6, "C" },
00060 { 7, "N" }, { 8, "O" }, { 9, "F" }, { 10, "Ne" }, { 11, "Na" }, { 12, "Mg" }, { 13, "Al" }, { 14, "Si" }, { 15, "P" },
00061 { 16, "S" }, { 17, "Cl" }, { 18, "Ar" }, { 19, "K" }, { 20, "Ca" }, { 21, "Sc" }, { 22, "Ti" }, { 23, "V" }, { 24, "Cr" },
00062 { 25, "Mn" }, { 26, "Fe" }, { 27, "Co" }, { 28, "Ni" }, { 29, "Cu" }, { 30, "Zn" }, { 31, "Ga" }, { 32, "Ge" }, { 33, "As" },
00063 { 34, "Se" }, { 35, "Br" }, { 36, "Kr" }, { 37, "Rb" }, { 38, "Sr" }, { 39, "Y" }, { 40, "Zr" }, { 41, "Nb" }, { 42, "Mo" },
00064 { 43, "Tc" }, { 44, "Ru" }, { 45, "Rh" }, { 46, "Pd" }, { 47, "Ag" }, { 48, "Cd" }, { 49, "In" }, { 50, "Sn" }, { 51, "Sb" },
00065 { 52, "Te" }, { 53, "I" }, { 54, "Xe" }, { 55, "Cs" }, { 56, "Ba" }, { 57, "La" }, { 58, "Ce" }, { 59, "Pr" }, { 60, "Nd" },
00066 { 61, "Pm" }, { 62, "Sm" }, { 63, "Eu" }, { 64, "Gd" }, { 65, "Tb" }, { 66, "Dy" }, { 67, "Ho" }, { 68, "Er" }, { 69, "Tm" },
00067 { 70, "Yb" }, { 71, "Lu" }, { 72, "Hf" }, { 73, "Ta" }, { 74, "W" }, { 75, "Re" }, { 76, "Os" }, { 77, "Ir" }, { 78, "Pt" },
00068 { 79, "Au" }, { 80, "Hg" }, { 81, "Tl" }, { 82, "Pb" }, { 83, "Bi" }, { 84, "Po" }, { 85, "At" }, { 86, "Rn" }, { 87, "Fr" },
00069 { 88, "Ra" }, { 89, "Ac" }, { 90, "Th" }, { 91, "Pa" }, { 92, "U" }, { 93, "Np" }, { 94, "Pu" }, { 95, "Am" }, { 96, "Cm" },
00070 { 97, "Bk" }, { 98, "Cf" }, { 99, "Es" }, { 100, "Fm" }, { 101, "Md" }, { 102, "No" }, { 103, "Lr" }, { 104, "Rf" }, { 105, "Db" },
00071 { 106, "Sg" }, { 107, "Bh" }, { 108, "Hs" }, { 109, "Mt" } };
00072
00073
00074
00075
00076 int tpia_misc_NumberOfZSymbols( void ) {
00077
00078 return( sizeof( ZSymbols ) / sizeof( struct ZSymbol ) );
00079 }
00080
00081
00082
00083 const char *tpia_misc_ZToSymbol( int iZ ) {
00084
00085 if( ( iZ < 0 ) || ( iZ >= tpia_misc_NumberOfZSymbols( ) ) ) return( NULL );
00086
00087 return( ZSymbols[iZ].symbol.c_str() );
00088 }
00089
00090
00091
00092 int tpia_misc_symbolToZ( const char *Z ) {
00093
00094 int i, n = tpia_misc_NumberOfZSymbols( );
00095
00096 for( i = 0; i < n; i++ ) {
00097
00098 if( !strcmp( Z, ZSymbols[i].symbol.c_str() ) ) return( ZSymbols[i].Z );
00099 }
00100 return( -1 );
00101 }
00102
00103
00104
00105 int tpia_miscNameToZAm( statusMessageReporting *smr, const char *name, int *Z, int *A, int *m ) {
00106
00107 int i, n;
00108 const char *p;
00109 char s[1024] = "", *q, *e;
00110
00111 n = sizeof( s ) - 1;
00112
00113 *Z = *A = *m = 0;
00114 if( !strncmp( "FissionProduct", name, 14 ) ) {
00115 *Z = 99;
00116 *A = 120;
00117 return( 0 );
00118 }
00119 if( !strcmp( "gamma", name ) ) return( 0 );
00120 for( p = name, q = s, i = 0; ( *p != '_' ) && ( i != n ); p++, q++, i++ ) *q = *p;
00121 if( i == n ) {
00122 smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Failed to find first '_' in particle name %s", name ); }
00123 else {
00124 *q = 0;
00125 if( ( *Z = tpia_misc_symbolToZ( s ) ) < 0 ) {
00126 smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Particle %s's symbol = '%s' not found", name, s ); }
00127 else {
00128 for( p++, q = s; ( *p != '_' ) && ( *p != 0 ) && ( i != n ); p++, q++, i++ ) *q = *p;
00129 if( i == n ) {
00130 smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Failed to find second '_' in particle name %s", name ); }
00131 else {
00132 *q = 0;
00133 if( strcmp( s, "natural" ) == 0 ) {
00134 e = s;
00135 while( *e ) e++; }
00136 else {
00137 *A = (int) strtol( s, &e, 10 );
00138 }
00139 if( *e != 0 ) {
00140 smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Failed to convert A to integer in particle name %s", name ); }
00141 else {
00142 *m = 0;
00143 if( *p == '_' ) {
00144 p++;
00145 if( *p != 'm' ) {
00146 smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Particle name %s missing meta-stable label 'm'", name ); }
00147 else {
00148 p++;
00149 *m = (int) strtol( p, &e, 10 );
00150 if( *e != 0 ) smr_setMessageError( smr, NULL, __FILE__, __LINE__, 1, "Failed to convert m to integer in particle name %s", name );
00151 }
00152 }
00153 }
00154 }
00155 }
00156 }
00157
00158 return( !smr_isOk( smr ) );
00159 }
00160
00161
00162
00163 char *tpia_misc_pointerToAttributeIfAllOk( statusMessageReporting *smr, xData_element *element, const char *path, int required,
00164 xData_attributionList *attributes, const char *name, const char *file, int line ) {
00165
00166 char *value;
00167
00168 if( !smr_isOk( smr ) ) return( NULL );
00169 if( ( value = xData_getAttributesValue( attributes, name ) ) == NULL ) {
00170 if( required ) {
00171 if( element != NULL ) {
00172 tpia_misc_setMessageError_Element( smr, NULL, element, file, line, 1, "element does not have attribute named %s", name ); }
00173 else {
00174 smr_setMessageError( smr, NULL, file, line, 1, "element does not have attribute named %s for file = %d", name, path );
00175 }
00176 }
00177 }
00178 return( value );
00179 }
00180
00181
00182
00183 int tpia_misc_setMessageError_Element( statusMessageReporting *smr, void *userInterface, xData_element *element, const char *file, int line, int code,
00184 const char *fmt, ... ) {
00185
00186 int status = 0;
00187 va_list args;
00188 char *msg;
00189
00190 va_start( args, fmt );
00191 msg = smr_vallocateFormatMessage( fmt, &args );
00192 va_end( args );
00193 if( msg == NULL ) {
00194 status = 1;
00195 va_start( args, fmt );
00196 smr_vsetMessageError( smr, userInterface, file, line, code, fmt, &args );
00197 va_end( args ); }
00198 else {
00199 status = smr_setMessageError( smr, userInterface, file, line, code, "%s for element %s at line %d column %d", msg, element->fullName,
00200 (int) element->docInfo.line, (int) element->docInfo.column );
00201 free( msg );
00202 }
00203 return( status );
00204 }
00205
00206
00207
00208 xData_Int tpia_misc_binarySearch( xData_Int n, double *ds, double d ) {
00209
00210 xData_Int imin = 0, imid, imax = n - 1;
00211
00212 if( d < ds[0] ) return( -2 );
00213 if( d > ds[n-1] ) return( -1 );
00214 while( 1 ) {
00215 imid = ( imin + imax ) >> 1;
00216 if( imid == imin ) break;
00217 if( d < ds[imid] ) {
00218 imax = imid; }
00219 else {
00220 imin = imid;
00221 }
00222 }
00223 return( imin );
00224 }
00225
00226
00227
00228 double *tpia_misc_get2dx_y_data( statusMessageReporting *smr, xData_element *element, xData_Int *length ) {
00229
00230 xData_element *xDataElement;
00231 double *d = NULL;
00232 xData_Int length_;
00233
00234 xData_addToAccessed( smr, element, 1 );
00235
00236 if( ( xDataElement = xData_getOneElementByTagName( smr, element, (char*)"xData", 1 ) ) != NULL ) {
00237 xData_addToAccessed( smr, xDataElement, 1 );
00238 if( xData_is_2d_xy( smr, &(xDataElement->xDataTypeInfo), 1 ) ) {
00239 d = xData_2d_xy_allocateCopyData( smr, xDataElement, &length_ );
00240 *length = (xData_Int) length_;
00241 }
00242 }
00243 return( d );
00244 }
00245
00246
00247
00248 double *tpia_misc_get2dxindex_y_data( statusMessageReporting *smr, xData_element *element, xData_Int *start, xData_Int *end, double *xValues ) {
00249
00250 xData_element *xDataElement;
00251 double *d = NULL;
00252
00253 xData_addToAccessed( smr, element, 1 );
00254
00255 if( ( xDataElement = xData_getOneElementByTagName( smr, element, (char*) "xData", 1 ) ) != NULL ) {
00256 xData_addToAccessed( smr, xDataElement, 1 );
00257 if( xData_is_2d_xindex_y( smr, &(xDataElement->xDataTypeInfo), 1 ) ) {
00258 if( start != NULL ) *start = xDataElement->xDataTypeInfo.start;
00259 if( end != NULL ) *end = xDataElement->xDataTypeInfo.end;
00260 d = xData_2d_xindex_y_toFilledYs( smr, xDataElement, xValues );
00261 }
00262 }
00263 return( d );
00264 }
00265
00266
00267
00268 double *tpia_misc_get2d_xShared_yHistogram_data( statusMessageReporting *smr, xData_element *element, xData_Int *start, xData_Int *end, xData_Int *length ) {
00269
00270 xData_Int i;
00271 xData_element *xDataElement;
00272 double *d = NULL;
00273
00274 xData_addToAccessed( smr, element, 1 );
00275
00276 if( ( xDataElement = xData_getOneElementByTagName( smr, element, (char*) "xData", 1 ) ) != NULL ) {
00277 xData_addToAccessed( smr, xDataElement, 1 );
00278 if( ( d = xData_2d_xShared_yHistogram_copyData( smr, xDataElement, &i ) ) != NULL ) {
00279 if( start != NULL ) *start = xDataElement->xDataTypeInfo.start;
00280 if( end != NULL ) *end = xDataElement->xDataTypeInfo.end;
00281 if( length != NULL ) *length = xDataElement->xDataTypeInfo.length;
00282 }
00283 }
00284 return( d );
00285 }
00286
00287
00288
00289 int tpia_misc_get2d_xShared_yHistogram_data_Grouped( statusMessageReporting *smr, xData_element *element, tpia_1dData *group ) {
00290
00291 if( ( group->data = tpia_misc_get2d_xShared_yHistogram_data( smr, element, &(group->start), &(group->end), &(group->length) ) ) == NULL ) return( 1 );
00292 return( 0 );
00293 }
00294
00295
00296
00297
00298 double tpia_misc_getPointwiseCrossSectionAtE( statusMessageReporting *, tpia_1dData *crossSection, double *energyGrid, xData_Int index, double e_in ) {
00299
00300 double xsec = 0.0, e1, e2;
00301
00302 if( ( index >= crossSection->start ) && ( index < crossSection->end ) ) {
00303 e1 = energyGrid[index];
00304 e2 = energyGrid[index + 1];
00305 index -= crossSection->start;
00306 if( e1 == e2 ) {
00307 xsec = 0.5 * ( crossSection->data[index] + crossSection->data[index + 1] ); }
00308 else {
00309 xsec = ( crossSection->data[index] * ( e2 - e_in ) + crossSection->data[index + 1] * ( e_in - e1 ) ) / ( e2 - e1 );
00310 }
00311 }
00312 return( xsec );
00313 }
00314
00315
00316
00317 tpia_EqualProbableBinSpectrum *tpia_misc_getEqualProbableBin( statusMessageReporting *smr, xData_element *parent, xData_Int *n, xData_Int *nBins ) {
00318
00319 xData_element *element;
00320
00321 xData_addToAccessed( smr, parent, 1 );
00322
00323 if( ( element = xData_getOneElementByTagName( smr, parent, (char*) "equalProbableBins", 0 ) ) == NULL ) return( NULL );
00324 if( xData_convertAttributeTo_xData_Int( smr, element, "nBins", nBins ) != 0 ) {
00325 tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "missing or invalid nBins attribute" );
00326 return( NULL );
00327 }
00328 return( tpia_misc_getEqualProbableBins( smr, element, "energy", *nBins, n ) );
00329 }
00330
00331
00332
00333 tpia_EqualProbableBinSpectrum *tpia_misc_getEqualProbableBins( statusMessageReporting *smr, xData_element *parent, const char *name, xData_Int nBins,
00334 xData_Int *n ) {
00335
00336
00337 int i, j;
00338 xData_Int index, size;
00339 xData_elementList *list;
00340 xData_element *element, *xData;
00341 double *d;
00342 tpia_EqualProbableBinSpectrum *epbs = NULL, *epb;
00343
00344 xData_addToAccessed( smr, parent, 1 );
00345 list = xData_getElementsByTagNameAndSort( smr, parent, name, NULL, NULL );
00346 if( list->n == 0 ) {
00347 tpia_misc_setMessageError_Element( smr, NULL, parent, __FILE__, __LINE__, 1, "bins does not contain any %s elements", name ); }
00348 else {
00349 *n = list->n;
00350 size = list->n * ( sizeof( tpia_EqualProbableBinSpectrum ) + ( nBins + 1 ) * sizeof( double ) );
00351
00352 if( ( epbs = (tpia_EqualProbableBinSpectrum*) xData_malloc2( smr, size, 0, "energies" ) ) != NULL ) {
00353 d = (double *) &(epbs[list->n]);
00354 for( i = 0, epb = epbs; i < list->n; i++, epb++ ) {
00355 element = list->items[i].element;
00356 xData_addToAccessed( smr, element, 1 );
00357 if( xData_convertAttributeTo_xData_Int( smr, element, "index", &index ) != 0 ) {
00358 tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "missing or invalid index attribute" );
00359
00360 epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
00361 break;
00362 }
00363 if( index != i ) {
00364 tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "index = %lld is not incremental", index );
00365
00366 epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
00367 break;
00368 }
00369 if( ( j = xData_convertAttributeToDouble( smr, element, "value", &(epb->value) ) ) != 0 ) {
00370 if( j == 1 ) {
00371 tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "element does not have value attribute" ); }
00372 else {
00373 tpia_misc_setMessageError_Element( smr, NULL, element, __FILE__, __LINE__, 1, "failed to convert value attribute to double" );
00374 }
00375
00376 epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
00377 break;
00378 }
00379 if( ( xData = xData_getElements_xDataElement( smr, element ) ) == NULL ) {
00380
00381 epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
00382 break;
00383 }
00384 xData_addToAccessed( smr, xData, 1 );
00385 epb->index = index;
00386 epb->nBins = nBins;
00387 epb->bins = d;
00388 if( xData_1d_x_copyData( smr, xData, ( nBins + 1 ) * sizeof( double ), d ) != 0 ) {
00389
00390 epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
00391 break;
00392 }
00393 d += nBins + 1;
00394 }
00395 }
00396 }
00397 xData_freeElementList( smr, list );
00398 return( epbs );
00399 }
00400
00401
00402
00403 double tpia_misc_drng( double (*rng)( void * ), void *rngState ) {
00404
00405 double r;
00406
00407 if( rng != NULL ) {
00408 r = rng( rngState ); }
00409 else {
00410
00411 r = CLHEP::HepRandom::getTheEngine()->flat();
00412
00413 }
00414 return( r );
00415 }
00416
00417
00418
00419
00420 int tpia_misc_sampleEqualProbableBin( statusMessageReporting *, tpia_decaySamplingInfo *decaySamplingInfo, double e_in, int nBins,
00421 tpia_EqualProbableBinSpectra *binned, double *value_ ) {
00422
00423 int i, j, index1, index2, method = 0;
00424 double fE = 1., r, value1, value2, value3, P12, P23, value = -2.;
00425
00426 if( e_in <= binned->energies[0].value ) {
00427 index1 = 0;
00428 index2 = 0; }
00429 else if( e_in >= binned->energies[binned->numberOfEs-1].value ) {
00430 index1 = binned->numberOfEs - 1;
00431 index2 = binned->numberOfEs - 1; }
00432 else {
00433 for( i = 0; i < binned->numberOfEs - 1; i++ ) {
00434 if( e_in <= binned->energies[i].value ) break;
00435 }
00436 i--;
00437 index1 = i;
00438 index2 = i;
00439 if( e_in != binned->energies[i].value ) {
00440 index2++;
00441 fE = ( e_in - binned->energies[i].value ) / ( binned->energies[i+1].value - binned->energies[i].value );
00442 }
00443 }
00444 r = tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState );
00445 j = (int) (r * nBins);
00446 if( j >= nBins ) j = nBins - 1;
00447 if( j < 0 ) j = 0;
00448 r = tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState );
00449 if( tpia_samplingMethods_isLinear( decaySamplingInfo->samplingMethods->angular_equalProbableBinMethod ) ) {
00450 method = 1;
00451 if( ( ( j == 0 ) && ( r <= 0.5 ) ) || ( j == ( nBins - 1 ) && r > 0.5 ) ) method = 0;
00452 }
00453 if( method == 0 ) {
00454 value1 = ( 1. - fE ) * binned->energies[index1].bins[j] + fE * binned->energies[index2].bins[j];
00455 value2 = ( 1. - fE ) * binned->energies[index1].bins[j+1] + fE * binned->energies[index2].bins[j+1];
00456 value = ( 1. - r ) * value1 + r * value2; }
00457 else {
00458 if( r <= 0.5 ) j--;
00459 value1 = ( 1. - fE ) * binned->energies[index1].bins[j] + fE * binned->energies[index2].bins[j];
00460 value2 = ( 1. - fE ) * binned->energies[index1].bins[j+1] + fE * binned->energies[index2].bins[j+1];
00461 value3 = ( 1. - fE ) * binned->energies[index1].bins[j+2] + fE * binned->energies[index2].bins[j+2];
00462 P12 = 1. / ( value2 - value1 );
00463 P23 = 1. / ( value3 - value2 );
00464 r = tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState );
00465 if( 0.25 * ( 1.0 + 2.0 * ( value2 - value1 ) / ( value3 - value1 ) ) > r ) {
00466 P23 = 2. / ( value3 - value1 );
00467 value3 = value2; }
00468 else {
00469 P12 = 2. / ( value3 - value1 );
00470 value1 = value2;
00471 }
00472 r = tpia_misc_drng( decaySamplingInfo->rng, decaySamplingInfo->rngState );
00473 if( P23 != P12 ) r = ( -P12 + std::sqrt( P12 * P12 * ( 1. - r ) + r * P23 * P23 ) ) / ( P23 - P12 );
00474 value = 0.5 * ( value1 + value2 + r * ( value3 - value1 ) );
00475 }
00476 *value_ = value;
00477 return( 0 );
00478 }
00479
00480 #if defined __cplusplus
00481 }
00482 #endif