Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/MCGIDI_energy.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/lend/src/MCGIDI_energy.cc (Version 11.3.0) and /processes/hadronic/models/lend/src/MCGIDI_energy.cc (Version 4.0.p2)


  1 /*                                                  1 
  2 # <<BEGIN-copyright>>                             
  3 # <<END-copyright>>                               
  4 */                                                
  5 #include <string.h>                               
  6 #define _USE_MATH_DEFINES                         
  7 #include <cmath>                                  
  8                                                   
  9                                                   
 10 #include "MCGIDI_fromTOM.h"                       
 11 #include "MCGIDI_misc.h"                          
 12 #include "MCGIDI_private.h"                       
 13 #include <nf_specialFunctions.h>                  
 14                                                   
 15 #if defined __cplusplus                           
 16 #include "G4Exp.hh"                               
 17 #include "G4Log.hh"                               
 18 #include "G4Pow.hh"                               
 19 namespace GIDI {                                  
 20 using namespace GIDI;                             
 21 #endif                                            
 22                                                   
 23 static int MCGIDI_energy_parseWeightFromTOM( s    
 24 static int MCGIDI_energy_parseWeightedFunction    
 25 static int MCGIDI_energy_parseGeneralEvaporati    
 26 static int MCGIDI_energy_parseEvaporationFromT    
 27 static int MCGIDI_energy_parseWattFromTOM( sta    
 28 static int MCGIDI_energy_parseSimpleMaxwellian    
 29 static int MCGIDI_energy_parseMadlandNixFromTO    
 30 static nfu_status MCGIDI_energy_parseMadlandNi    
 31 static double MCGIDI_energy_parseMadlandNixFro    
 32 static int MCGIDI_energy_parseNBodyPhaseSpaceF    
 33         MCGIDI_distribution *distribution );      
 34                                                   
 35 static int MCGIDI_energy_sampleSimpleMaxwellia    
 36 static int MCGIDI_energy_sampleEvaporation( st    
 37 static int MCGIDI_energy_sampleWatt( statusMes    
 38 static int MCGIDI_energy_sampleWeightedFunctio    
 39     MCGIDI_quantitiesLookupModes &modes, MCGID    
 40 static nfu_status MCGIDI_energy_NBodyPhaseSpac    
 41 /*                                                
 42 **********************************************    
 43 */                                                
 44 MCGIDI_energy *MCGIDI_energy_new( statusMessag    
 45                                                   
 46     MCGIDI_energy *energy;                        
 47                                                   
 48     if( ( energy = (MCGIDI_energy *) smr_mallo    
 49     if( MCGIDI_energy_initialize( smr, energy     
 50     return( energy );                             
 51 }                                                 
 52 /*                                                
 53 **********************************************    
 54 */                                                
 55 int MCGIDI_energy_initialize( statusMessageRep    
 56                                                   
 57     memset( energy, 0, sizeof( MCGIDI_energy )    
 58     return( 0 );                                  
 59 }                                                 
 60 /*                                                
 61 **********************************************    
 62 */                                                
 63 MCGIDI_energy *MCGIDI_energy_free( statusMessa    
 64                                                   
 65     MCGIDI_energy_release( smr, energy );         
 66     smr_freeMemory( (void **) &energy );          
 67     return( NULL );                               
 68 }                                                 
 69 /*                                                
 70 **********************************************    
 71 */                                                
 72 int MCGIDI_energy_release( statusMessageReport    
 73                                                   
 74     int i;                                        
 75                                                   
 76     MCGIDI_sampling_pdfsOfXGivenW_release( smr    
 77     if( energy->theta ) energy->theta = ptwXY_    
 78     if( energy->Watt_a ) energy->Watt_a = ptwX    
 79     if( energy->Watt_b ) energy->Watt_b = ptwX    
 80     if( ( energy->type == MCGIDI_energyType_ge    
 81         MCGIDI_sampling_pdfsOfX_release( smr,     
 82     else if( energy->type == MCGIDI_energyType    
 83         for( i = 0; i < energy->weightedFuncti    
 84             ptwXY_free( energy->weightedFuncti    
 85             MCGIDI_energy_free( smr, energy->w    
 86         }                                         
 87     }                                             
 88                                                   
 89     MCGIDI_energy_initialize( smr, energy );      
 90     return( 0 );                                  
 91 }                                                 
 92 /*                                                
 93 **********************************************    
 94 */                                                
 95 int MCGIDI_energy_parseFromTOM( statusMessageR    
 96         enum MCGIDI_energyType energyType, dou    
 97                                                   
 98     MCGIDI_energy *energy = NULL;                 
 99     xDataTOM_element *energyElement, *linearEl    
100     char const *nativeData;                       
101     double projectileMass_MeV, targetMass_MeV;    
102                                                   
103     if( ( energy = MCGIDI_energy_new( smr ) )     
104                                                   
105     projectileMass_MeV = MCGIDI_product_getPro    
106     targetMass_MeV = MCGIDI_product_getTargetM    
107     energy->e_inCOMFactor = targetMass_MeV / (    
108                                                   
109     if( ( energyType == MCGIDI_energyType_prim    
110         energy->type = energyType;                
111         energy->gammaEnergy_MeV = gammaEnergy_    
112         energy->frame = xDataTOM_frame_lab;       
113         if( energyType == MCGIDI_energyType_pr    
114     else {                                        
115         if( ( energyElement = xDataTOME_getOne    
116         if( ( nativeData = xDataTOM_getAttribu    
117         if( ( linearElement = xDataTOME_getOne    
118             linearElement = xDataTOME_getOneEl    
119         if( linearElement == NULL ) {             
120             if( ( functional = xDataTOME_getOn    
121                 if( MCGIDI_energy_parseGeneral    
122             else if( ( functional = xDataTOME_    
123                 if( MCGIDI_energy_parseSimpleM    
124             else if( ( functional = xDataTOME_    
125                 if( MCGIDI_energy_parseEvapora    
126             else if( ( functional = xDataTOME_    
127                 if( MCGIDI_energy_parseWattFro    
128             else if( ( functional = xDataTOME_    
129                 if( MCGIDI_energy_parseMadland    
130             else if( ( functional = xDataTOME_    
131                 if( MCGIDI_energy_parseNBodyPh    
132             else if( ( functional = xDataTOME_    
133                 if( MCGIDI_energy_parseWeighte    
134             else {                                
135                 smr_setReportError2( smr, smr_    
136                 goto err;                         
137             }                                     
138             frameElement = functional; }          
139         else {                                    
140             char const *toUnits[3] = { "MeV",     
141                                                   
142             frameElement = linearElement;         
143             if( MCGIDI_fromTOM_pdfsOfXGivenW(     
144             energy->type = MCGIDI_energyType_l    
145         }                                         
146         if( ( energy->frame = MCGIDI_misc_getP    
147     }                                             
148     distribution->energy = energy;                
149                                                   
150     return( 0 );                                  
151                                                   
152 err:                                              
153     if( energy != NULL ) MCGIDI_energy_free( s    
154     return( 1 );                                  
155 }                                                 
156 /*                                                
157 **********************************************    
158 */                                                
159 static int MCGIDI_energy_parseWeightedFunction    
160                                                   
161     int i;                                        
162     xDataTOM_element *child;                      
163                                                   
164     for( i = 0, child = xDataTOME_getFirstElem    
165         if( strcmp( child->name, "weighted" )     
166         if( MCGIDI_energy_parseWeightFromTOM(     
167         energy->weightedFunctionals.numberOfWe    
168     }                                             
169     energy->type = MCGIDI_energyType_weightedF    
170     return( 0 );                                  
171                                                   
172 err:                                              
173     return( 1 );                                  
174 }                                                 
175 /*                                                
176 **********************************************    
177 */                                                
178 static int MCGIDI_energy_parseWeightFromTOM( s    
179                                                   
180     xDataTOM_element *child;                      
181     MCGIDI_energy *energy = NULL;                 
182     ptwXYPoints *weight = NULL;                   
183     char const *toUnits[2] = { "MeV", "" };       
184                                                   
185     if( ( energy = MCGIDI_energy_new( smr ) )     
186     for( child = xDataTOME_getFirstElement( el    
187         if( strcmp( child->name, "weight" ) ==    
188             if( ( weight = MCGIDI_misc_dataFro    
189         else if( strcmp( child->name, "evapora    
190             if( MCGIDI_energy_parseEvaporation    
191         else {                                    
192             smr_setReportError2( smr, smr_unkn    
193             goto err;                             
194         }                                         
195     }                                             
196     weightedFunctional->weight = weight;          
197     weightedFunctional->energy = energy;          
198     return( 0 );                                  
199                                                   
200 err:                                              
201     if( weight != NULL ) ptwXY_free( weight );    
202     if( energy != NULL ) MCGIDI_energy_free( s    
203     return( 1 );                                  
204 }                                                 
205 /*                                                
206 **********************************************    
207 */                                                
208 static int MCGIDI_energy_parseGeneralEvaporati    
209                                                   
210     double norm;                                  
211     xDataTOM_element *thetaTOM, *gTOM;            
212     ptwXYPoints *theta = NULL, *g = NULL;         
213     char const *toUnits[2] = { "MeV", "MeV" };    
214                                                   
215     if( ( thetaTOM = xDataTOME_getOneElementBy    
216     if( ( theta = MCGIDI_misc_dataFromElement2    
217                                                   
218     if( ( gTOM = xDataTOME_getOneElementByName    
219     toUnits[0] = "";                              
220     toUnits[1] = "";                              
221     if( ( g = MCGIDI_misc_dataFromElement2ptwX    
222     if( MCGIDI_fromTOM_pdfOfX( smr, g, &(energ    
223     energy->gInterpolation = ptwXY_getInterpol    
224     g = ptwXY_free( g );                          
225     if( std::fabs( 1. - norm ) > 0.001 ) print    
226                                                   
227     energy->type = MCGIDI_energyType_generalEv    
228     energy->theta = theta;                        
229     return( 0 );                                  
230                                                   
231 err:                                              
232     if( theta != NULL ) ptwXY_free( theta );      
233     if( g != NULL ) ptwXY_free( g );              
234     return( 1 );                                  
235 }                                                 
236 /*                                                
237 **********************************************    
238 */                                                
239 static int MCGIDI_energy_parseSimpleMaxwellian    
240                                                   
241     char const *U, *toUnits[2] = { "MeV", "MeV    
242     xDataTOM_element *thetaTOM;                   
243                                                   
244     if( ( U = xDataTOM_getAttributesValueInEle    
245         smr_setReportError2( smr, smr_unknownI    
246         goto err;                                 
247     }                                             
248     if( MCGIDI_misc_PQUStringToDoubleInUnitOf(    
249     if( ( thetaTOM = xDataTOME_getOneElementBy    
250     if( ( energy->theta = MCGIDI_misc_dataFrom    
251     energy->type = MCGIDI_energyType_simpleMax    
252     return( 0 );                                  
253                                                   
254 err:                                              
255     return( 1 );                                  
256 }                                                 
257 /*                                                
258 **********************************************    
259 */                                                
260 static int MCGIDI_energy_parseEvaporationFromT    
261                                                   
262     char const *U, *toUnits[2] = { "MeV", "MeV    
263     xDataTOM_element *thetaTOM;                   
264                                                   
265     if( ( U = xDataTOM_getAttributesValueInEle    
266         smr_setReportError2( smr, smr_unknownI    
267         goto err;                                 
268     }                                             
269     if( MCGIDI_misc_PQUStringToDoubleInUnitOf(    
270     if( ( thetaTOM = xDataTOME_getOneElementBy    
271     if( ( energy->theta = MCGIDI_misc_dataFrom    
272     energy->type = MCGIDI_energyType_evaporati    
273     return( 0 );                                  
274                                                   
275 err:                                              
276     return( 1 );                                  
277 }                                                 
278 /*                                                
279 **********************************************    
280 */                                                
281 static int MCGIDI_energy_parseWattFromTOM( sta    
282                                                   
283     char const *U, *toUnits[2] = { "MeV", "MeV    
284     xDataTOM_element *aOrBTOM;                    
285                                                   
286     if( ( U = xDataTOM_getAttributesValueInEle    
287         smr_setReportError2( smr, smr_unknownI    
288         goto err;                                 
289     }                                             
290     if( MCGIDI_misc_PQUStringToDoubleInUnitOf(    
291                                                   
292     if( ( aOrBTOM = xDataTOME_getOneElementByN    
293     if( ( energy->Watt_a = MCGIDI_misc_dataFro    
294                                                   
295     toUnits[1] = "1/MeV";                         
296     if( ( aOrBTOM = xDataTOME_getOneElementByN    
297     if( ( energy->Watt_b = MCGIDI_misc_dataFro    
298                                                   
299     energy->type = MCGIDI_energyType_Watt;        
300     return( 0 );                                  
301                                                   
302 err:                                              
303     return( 1 );                                  
304 }                                                 
305 /*                                                
306 **********************************************    
307 */                                                
308 static int MCGIDI_energy_parseMadlandNixFromTO    
309                                                   
310     int iE, length, nXs, i1, n;                   
311     double E=0., T_M=0., EFL=0., EFH=0., argLi    
312            xs[] = { 1e-5, 1e-3, 1e-1, 1e1, 1e3    
313     ptwXYPoints *ptwXY_TM = NULL, *pdfXY = NUL    
314     ptwXYPoint *point;                            
315     ptwXPoints *cdfX = NULL;                      
316     nfu_status status = nfu_Okay;                 
317     xDataTOM_element *TM_TOM;                     
318     xDataTOM_XYs *XYs;                            
319     MCGIDI_pdfsOfXGivenW *dists = &(energy->di    
320     MCGIDI_pdfOfX *dist;                          
321     char const *EF, *TMUnits[2] = { "MeV", "Me    
322                                                   
323     nXs = sizeof( xs ) / sizeof( xs[0] );         
324                                                   
325     if( ( EF = xDataTOM_getAttributesValueInEl    
326         smr_setReportError2( smr, smr_unknownI    
327         goto err;                                 
328     }                                             
329     if( MCGIDI_misc_PQUStringToDoubleInUnitOf(    
330     argList[0] = EFL;                             
331                                                   
332     if( ( EF = xDataTOM_getAttributesValueInEl    
333         smr_setReportError2( smr, smr_unknownI    
334         goto err;                                 
335     }                                             
336     if( MCGIDI_misc_PQUStringToDoubleInUnitOf(    
337     argList[1] = EFH;                             
338                                                   
339     if( ( TM_TOM = xDataTOME_getOneElementByNa    
340     if( ( XYs = (xDataTOM_XYs *) xDataTOME_get    
341     if( ( ptwXY_TM = MCGIDI_misc_dataFromXYs2p    
342                                                   
343     length = (int) ptwXY_length( ptwXY_TM );      
344     dists->interpolationWY = ptwXY_interpolati    
345     dists->interpolationXY = ptwXY_interpolati    
346     if( ( dists->Ws = (double *) smr_malloc2(     
347     if( ( dists->dist = (MCGIDI_pdfOfX *) smr_    
348                                                   
349     for( iE = 0; iE < length; iE++ ) {            
350         ptwXY_getXYPairAtIndex( ptwXY_TM, iE,     
351         argList[2] = T_M;                         
352         dist = &(dists->dist[iE]);                
353         dists->Ws[iE] = E;                        
354                                                   
355         if( ( pdfXY = ptwXY_createFromFunction    
356             (void *) argList, 1e-3, 0, 12, &st    
357         if( ( status = ptwXY_normalize( pdfXY     
358             smr_setReportError2( smr, smr_unkn    
359             goto err;                             
360         }                                         
361                                                   
362         if( ptwXY_simpleCoalescePoints( pdfXY     
363         dist->numberOfXs = n = (int) ptwXY_len    
364                                                   
365         if( ( dist->Xs = (double *) smr_malloc    
366         dists->numberOfWs++;                      
367         dist->pdf = &(dist->Xs[n]);               
368         dist->cdf = &(dist->pdf[n]);              
369                                                   
370         for( i1 = 0; i1 < n; i1++ ) {             
371             point = ptwXY_getPointAtIndex_Unsa    
372             dist->Xs[i1] = point->x;              
373             dist->pdf[i1] = point->y;             
374         }                                         
375                                                   
376         if( ( cdfX = ptwXY_runningIntegral( pd    
377             smr_setReportError2( smr, smr_unkn    
378             goto err;                             
379         }                                         
380                                                   
381         norm = ptwX_getPointAtIndex_Unsafely(     
382         for( i1 = 0; i1 < n; i1++ ) dist->cdf[    
383         for( i1 = 0; i1 < n; i1++ ) dist->pdf[    
384         pdfXY = ptwXY_free( pdfXY );              
385         cdfX = ptwX_free( cdfX );                 
386     }                                             
387                                                   
388     energy->type = MCGIDI_energyType_MadlandNi    
389                                                   
390     ptwXY_free( ptwXY_TM );                       
391     return( 0 );                                  
392                                                   
393 err:                                              
394     if( ptwXY_TM != NULL ) ptwXY_free( ptwXY_T    
395     if( pdfXY != NULL ) ptwXY_free( pdfXY );      
396     if( cdfX != NULL ) cdfX = ptwX_free( cdfX     
397                                                   
398     return( 1 );                                  
399 }                                                 
400 /*                                                
401 **********************************************    
402 */                                                
403 static nfu_status MCGIDI_energy_parseMadlandNi    
404                                                   
405     double *parameters = (double *) argList, E    
406     nfu_status status = nfu_Okay;                 
407                                                   
408     EFL = parameters[0];                          
409     EFH = parameters[1];                          
410     T_M = parameters[2];                          
411     *y = MCGIDI_energy_parseMadlandNixFromTOM_    
412     if( status == nfu_Okay ) *y += MCGIDI_ener    
413     *y *= 0.5;                                    
414     return( status );                             
415 }                                                 
416 /*                                                
417 **********************************************    
418 */                                                
419 static double MCGIDI_energy_parseMadlandNixFro    
420                                                   
421     double u1, u2, E1, E2 = 0., gamma1 = 0., g    
422                                                   
423     u1 = std::sqrt( Ep ) - std::sqrt( E_F );      
424     u1 *= u1 / T_M;                               
425     u2 = std::sqrt( Ep ) + std::sqrt( E_F );      
426     u2 *= u2 / T_M;                               
427     E1 = 0;                      /* u1^3/2 * E    
428     if( u1 != 0 ) E1 = nf_exponentialIntegral(    
429     if( *status == nfu_Okay ) E2 = nf_exponent    
430     if( *status != nfu_Okay ) return( 0. );       
431     if( u1 > 2. ) {                               
432         signG = -1;                               
433         gamma1 = nf_incompleteGammaFunctionCom    
434         if( *status == nfu_Okay ) gamma2 = nf_    
435     else {                                        
436         gamma1 = nf_incompleteGammaFunction( 1    
437         if( *status == nfu_Okay ) gamma2 = nf_    
438     }                                             
439     if( *status != nfu_Okay ) return( 0. );       
440     return( ( u2 * std::sqrt( u2 ) * E2 - u1 *    
441 }                                                 
442 /*                                                
443 **********************************************    
444 */                                                
445 static int MCGIDI_energy_parseNBodyPhaseSpaceF    
446         MCGIDI_distribution *distribution ) {     
447                                                   
448     int argList[1];                               
449     double xs[2] = { 0.0, 1.0 }, productMass_M    
450     ptwXYPoints *pdf = NULL;                      
451     nfu_status status;                            
452     char const *mass;                             
453                                                   
454     if( xDataTOME_convertAttributeToInteger( N    
455     if( ( mass = xDataTOM_getAttributesValueIn    
456         smr_setReportError2( smr, smr_unknownI    
457         goto err;                                 
458     }                                             
459     if( MCGIDI_misc_PQUStringToDouble( smr, ma    
460     argList[0] = energy->NBodyPhaseSpace.numbe    
461     if( ( pdf = ptwXY_createFromFunction( 2, x    
462         smr_setReportError2( smr, smr_unknownI    
463             status, nfu_statusMessage( status     
464         goto err;                                 
465     }                                             
466     if( MCGIDI_fromTOM_pdfOfX( smr, pdf, &(ene    
467     productMass_MeV = MCGIDI_product_getMass_M    
468     if( !smr_isOk( smr ) ) goto err;              
469     energy->NBodyPhaseSpace.massFactor = ( 1.     
470     energy->NBodyPhaseSpace.Q_MeV = MCGIDI_out    
471     if( !smr_isOk( smr ) ) goto err;              
472                                                   
473     ptwXY_free( pdf );                            
474     energy->type = MCGIDI_energyType_NBodyPhas    
475                                                   
476     return( 0 );                                  
477                                                   
478 err:                                              
479     if( pdf != NULL ) ptwXY_free( pdf );          
480     return( 1 );                                  
481 }                                                 
482 /*                                                
483 **********************************************    
484 */                                                
485 static nfu_status MCGIDI_energy_NBodyPhaseSpac    
486                                                   
487     int numberOfProducts = *((int *) argList);    
488     double e = 0.5 * ( 3 * numberOfProducts -     
489                                                   
490     *y = std::sqrt( x ) * G4Pow::GetInstance()    
491     return( nfu_Okay );                           
492 }                                                 
493 /*                                                
494 **********************************************    
495 */                                                
496 int MCGIDI_energy_sampleEnergy( statusMessageR    
497         MCGIDI_decaySamplingInfo *decaySamplin    
498 /*                                                
499 *   This function must be called before angula    
500 */                                                
501     double theta, randomEp, Watt_a, Watt_b, e_    
502     MCGIDI_pdfsOfXGivenW_sampled sampled;         
503                                                   
504     decaySamplingInfo->frame = energy->frame;     
505     switch( energy->type ) {                      
506     case MCGIDI_energyType_primaryGamma :         
507         decaySamplingInfo->Ep = energy->gammaE    
508         break;                                    
509     case MCGIDI_energyType_discreteGamma :        
510         decaySamplingInfo->Ep = energy->gammaE    
511         break;                                    
512     case MCGIDI_energyType_linear :               
513         randomEp = decaySamplingInfo->rng( dec    
514         sampled.smr = smr;                        
515         sampled.w = e_in;                         
516         MCGIDI_sampling_sampleX_from_pdfsOfXGi    
517         decaySamplingInfo->Ep = sampled.x;        
518         break;                                    
519     case MCGIDI_energyType_generalEvaporation     
520         sampled.interpolationXY = energy->gInt    
521         MCGIDI_sampling_sampleX_from_pdfOfX( &    
522         theta = MCGIDI_sampling_ptwXY_getValue    
523         decaySamplingInfo->Ep = theta * sample    
524         break;                                    
525     case MCGIDI_energyType_simpleMaxwellianFis    
526         theta = MCGIDI_sampling_ptwXY_getValue    
527         MCGIDI_energy_sampleSimpleMaxwellianFi    
528         decaySamplingInfo->Ep *= theta;           
529         break;                                    
530     case MCGIDI_energyType_evaporation :          
531         theta = MCGIDI_sampling_ptwXY_getValue    
532         MCGIDI_energy_sampleEvaporation( smr,     
533         decaySamplingInfo->Ep *= theta;           
534         break;                                    
535     case MCGIDI_energyType_Watt :                 
536         Watt_a = MCGIDI_sampling_ptwXY_getValu    
537         Watt_b = MCGIDI_sampling_ptwXY_getValu    
538         MCGIDI_energy_sampleWatt( smr, e_in -     
539         break;                                    
540     case MCGIDI_energyType_MadlandNix :           
541         MCGIDI_sampling_sampleX_from_pdfsOfXGi    
542         decaySamplingInfo->Ep = sampled.x;        
543         break;                                    
544     case MCGIDI_energyType_NBodyPhaseSpace :      
545         MCGIDI_sampling_sampleX_from_pdfOfX( &    
546         decaySamplingInfo->Ep = ( energy->e_in    
547         break;                                    
548     case MCGIDI_energyType_weightedFunctional     
549         MCGIDI_energy_sampleWeightedFunctional    
550         break;                                    
551     default :                                     
552         smr_setReportError2( smr, smr_unknownI    
553     }                                             
554                                                   
555     return( !smr_isOk( smr ) );                   
556 }                                                 
557 /*                                                
558 **********************************************    
559 */                                                
560 static int MCGIDI_energy_sampleSimpleMaxwellia    
561                                                   
562     int i1;                                       
563     double a = e_in_U_theta, b, c, x, norm_a,     
564                                                   
565     sqrt_x = std::sqrt( a );                      
566     norm_a = sqrt_pi_2 * erf( sqrt_x ) - sqrt_    
567     b = norm_a * decaySamplingInfo->rng( decay    
568     for( i1 = 0; i1 < 16; i1++ ) {                
569         x = 0.5 * ( xMin + xMax );                
570         sqrt_x = std::sqrt( x );                  
571         c = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x    
572         if( b < c ) {                             
573             xMax = x; }                           
574         else {                                    
575             xMin = x;                             
576         }                                         
577     }                                             
578         /* To order e, the correct x is x + e     
579     decaySamplingInfo->Ep = x;                    
580                                                   
581     return( 0 );                                  
582 }                                                 
583 /*                                                
584 **********************************************    
585 */                                                
586 static int MCGIDI_energy_sampleEvaporation( st    
587                                                   
588     int i1;                                       
589     double a = e_in_U_theta, b, c, x, norm_a,     
590                                                   
591     norm_a = 1 - ( 1 + a ) * G4Exp( -a );         
592     b = 1. - norm_a * decaySamplingInfo->rng(     
593     for( i1 = 0; i1 < 16; i1++ ) {                
594         x = 0.5 * ( xMin + xMax );                
595         c = ( 1 + x ) * G4Exp( -x );              
596         if( b > c ) {                             
597             xMax = x; }                           
598         else {                                    
599             xMin = x;                             
600         }                                         
601     }                                             
602         /* To order e, the correct x is x + e     
603     decaySamplingInfo->Ep = x;                    
604                                                   
605     return( 0 );                                  
606 }                                                 
607 /*                                                
608 **********************************************    
609 */                                                
610 static int MCGIDI_energy_sampleWatt( statusMes    
611 /*                                                
612 *   From MCAPM via Sample Watt Spectrum as in     
613 */                                                
614     double WattMin = 0., WattMax = e_in_U, x,     
615                                                   
616     x = 1. + ( Watt_b / ( 8. * Watt_a ) );        
617     y = ( x + std::sqrt( x * x - 1. ) ) / Watt    
618     z = Watt_a * y - 1.;                          
619    G4int icounter=0;                              
620     G4int icounter_max=1024;                      
621     do                                            
622     {                                             
623       icounter++;                                 
624       if ( icounter > icounter_max ) {            
625          G4cout << "Loop-counter exceeded the     
626          break;                                   
627       }                                           
628         rand1 = -G4Log( decaySamplingInfo->rng    
629         rand2 = -G4Log( decaySamplingInfo->rng    
630         energyOut = y * rand1;                    
631     }                                             
632     while( ( ( rand2 - z * ( rand1 + 1. ) ) *     
633     decaySamplingInfo->Ep = energyOut;            
634                                                   
635     return( 0 );                                  
636 }                                                 
637 /*                                                
638 **********************************************    
639 */                                                
640 static int MCGIDI_energy_sampleWeightedFunctio    
641         MCGIDI_quantitiesLookupModes &modes, M    
642 /*                                                
643 c   This routine assumes that the weights sum     
644 */                                                
645     int iW;                                       
646     double rW = decaySamplingInfo->rng( decayS    
647     MCGIDI_energyWeightedFunctional *weightedF    
648                                                   
649     for( iW = 0; iW < energy->weightedFunction    
650         weightedFunctional = &(energy->weighte    
651         weight = MCGIDI_sampling_ptwXY_getValu    
652         cumulativeW += weight;                    
653         if( cumulativeW >= rW ) break;            
654     }                                             
655     return( MCGIDI_energy_sampleEnergy( smr, w    
656 }                                                 
657                                                   
658 #if defined __cplusplus                           
659 }                                                 
660 #endif                                            
661                                                   
662