tpia_misc.cc

Go to the documentation of this file.
00001 /*
00002 # <<BEGIN-copyright>>
00003 # Copyright (c) 2010, Lawrence Livermore National Security, LLC. 
00004 # Produced at the Lawrence Livermore National Laboratory 
00005 # Written by Bret R. Beck, beck6@llnl.gov. 
00006 # CODE-461393
00007 # All rights reserved. 
00008 #  
00009 # This file is part of GIDI. For details, see nuclear.llnl.gov. 
00010 # Please also read the "Additional BSD Notice" at nuclear.llnl.gov. 
00011 # 
00012 # Redistribution and use in source and binary forms, with or without modification, 
00013 # are permitted provided that the following conditions are met: 
00014 #
00015 #      1) Redistributions of source code must retain the above copyright notice, 
00016 #         this list of conditions and the disclaimer below.
00017 #      2) Redistributions in binary form must reproduce the above copyright notice, 
00018 #         this list of conditions and the disclaimer (as noted below) in the 
00019 #          documentation and/or other materials provided with the distribution.
00020 #      3) Neither the name of the LLNS/LLNL nor the names of its contributors may be 
00021 #         used to endorse or promote products derived from this software without 
00022 #         specific prior written permission. 
00023 #
00024 # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 
00025 # EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 
00026 # OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT 
00027 # SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR 
00028 # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
00029 # CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 
00030 # OR SERVICES;  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 
00031 # AND ON  ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
00032 # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, 
00033 # EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
00034 # <<END-copyright>>
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     //char *symbol;
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     //return( ZSymbols[iZ].symbol );
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         //if( !strcmp( Z, ZSymbols[i].symbol ) ) return( ZSymbols[i].Z );
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;   /* Note 1) routine will fail when parts of a particle name can be longer than 1024. */
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 ) {              /* See note 1 above. */
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 {                  /* Getting here implies that *p == '_'. */
00128             for( p++, q = s; ( *p != '_' ) && ( *p != 0 ) && ( i != n ); p++, q++, i++ ) *q = *p;
00129             if( i == n ) {      /* See note 1 above. */
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 {          /* Getting here implies that *p == '_' or 0. */
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     //if( ( xDataElement = xData_getOneElementByTagName( smr, element, "xData", 1 ) ) != NULL ) {
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     //if( ( xDataElement = xData_getOneElementByTagName( smr, element, "xData", 1 ) ) != NULL ) {
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     //if( ( xDataElement = xData_getOneElementByTagName( smr, element, "xData", 1 ) ) != NULL ) {
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 //double tpia_misc_getPointwiseCrossSectionAtE( statusMessageReporting *smr, tpia_1dData *crossSection, double *energyGrid, xData_Int index, double e_in ) {
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 ) {                            /* Allow for future where step function may be allowed. */
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     //if( ( element = xData_getOneElementByTagName( smr, parent, "equalProbableBins", 0 ) ) == NULL ) return( NULL );
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         //if( ( epbs = xData_malloc2( smr, size, 0, "energies" ) ) != NULL ) {
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++ ) {    /* Loop to test nBins and index are proper. */
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                     //epbs = xData_free( smr, epbs );
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                     //epbs = xData_free( smr, epbs );
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                     //epbs = xData_free( smr, epbs );
00376                     epbs = (tpia_EqualProbableBinSpectrum*) xData_free( smr, epbs );
00377                     break;
00378                 }
00379                 if( ( xData = xData_getElements_xDataElement( smr, element ) ) == NULL ) {
00380                     //epbs = xData_free( smr, epbs );
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                     //epbs = xData_free( smr, epbs );
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         //r = drand48( );
00411         r = CLHEP::HepRandom::getTheEngine()->flat();
00412           
00413     }
00414     return( r );
00415 }
00416 /*
00417 ************************************************************
00418 */
00419 //int tpia_misc_sampleEqualProbableBin( statusMessageReporting *smr, tpia_decaySamplingInfo *decaySamplingInfo, double e_in, int nBins, 
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 );          // Do not change r until after Point1 below.
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 ) {                 /* Constant interpolaton. */
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 {                              /* Linear interpolation. */
00458         if( r <= 0.5 ) j--;             /* Point1: Above test insures that j is not 0 (nBins-1) if r <= 0.5 (> 0.5); */
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

Generated on Mon May 27 17:50:35 2013 for Geant4 by  doxygen 1.4.7