Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MCGIDI_sampling.cc File Reference
#include <stdlib.h>
#include <cmath>
#include "MCGIDI.h"
Include dependency graph for MCGIDI_sampling.cc:

Go to the source code of this file.

Functions

int MCGIDI_sampling_pdfsOfXGivenW_initialize (statusMessageReporting *, MCGIDI_pdfsOfXGivenW *dists)
 
int MCGIDI_sampling_pdfsOfXGivenW_release (statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *dists)
 
int MCGIDI_sampling_pdfsOfX_release (statusMessageReporting *, MCGIDI_pdfOfX *dist)
 
int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW (MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue)
 
int MCGIDI_sampling_sampleX_from_pdfOfX (MCGIDI_pdfOfX *dist, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue)
 
int MCGIDI_sampling_doubleDistribution (statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *pdfOfWGivenV, MCGIDI_pdfsOfXGivenW *pdfOfXGivenVAndW, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
 
int MCGIDI_sampling_interpolationValues (statusMessageReporting *smr, ptwXY_interpolation interpolation, double *ws, double y1, double y2, double *y)
 
double MCGIDI_sampling_ptwXY_getValueAtX (ptwXYPoints *ptwXY, double x1)
 

Function Documentation

int MCGIDI_sampling_doubleDistribution ( statusMessageReporting smr,
MCGIDI_pdfsOfXGivenW pdfOfWGivenV,
MCGIDI_pdfsOfXGivenW pdfOfXGivenVAndW,
MCGIDI_quantitiesLookupModes modes,
MCGIDI_decaySamplingInfo decaySamplingInfo 
)

Definition at line 141 of file MCGIDI_sampling.cc.

142  {
143 
144  int iV;
145  double e_in = modes.getProjectileEnergy( );
146  double randomW = decaySamplingInfo->rng( decaySamplingInfo->rngState ), randomX = decaySamplingInfo->rng( decaySamplingInfo->rngState );
147  MCGIDI_pdfsOfXGivenW_sampled sampledX, sampledW;
148  ptwXY_interpolation interpolationWY = pdfOfWGivenV->interpolationWY;
149 
150  sampledX.smr = smr;
151  sampledW.smr = smr;
152  sampledW.interpolationXY = pdfOfWGivenV->interpolationXY;
153  iV = MCGIDI_misc_binarySearch( pdfOfWGivenV->numberOfWs, pdfOfWGivenV->Ws, e_in );
154  if( iV < 0 ) {
155  interpolationWY = ptwXY_interpolationFlat;
156  if( iV == -2 ) {
157  iV = 0; }
158  else {
159  iV = pdfOfWGivenV->numberOfWs - 1;
160  }
161  e_in = pdfOfWGivenV->Ws[iV];
162  }
163 
164  MCGIDI_sampling_sampleX_from_pdfOfX( &(pdfOfWGivenV->dist[iV]), &sampledW, randomW );
165  sampledX.w = sampledW.x;
166  MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(pdfOfXGivenVAndW[iV]), &sampledX, randomX );
167 
168  if( interpolationWY != ptwXY_interpolationFlat ) {
169  double x = sampledX.x, w = sampledW.x, Vs[3] = { e_in, pdfOfWGivenV->Ws[iV], pdfOfWGivenV->Ws[iV+1] };
170 
171  MCGIDI_sampling_sampleX_from_pdfOfX( &(pdfOfWGivenV->dist[iV+1]), &sampledW, randomW );
172  sampledX.w = sampledW.x;
173  MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(pdfOfXGivenVAndW[iV+1]), &sampledX, randomX );
174 
175  MCGIDI_sampling_interpolationValues( smr, interpolationWY, Vs, w, sampledW.x, &sampledW.x );
176  MCGIDI_sampling_interpolationValues( smr, interpolationWY, Vs, x, sampledX.x, &sampledX.x );
177  }
178 
179  decaySamplingInfo->mu = sampledW.x;
180  decaySamplingInfo->Ep = sampledX.x;
181 
182  return( 0 );
183 }
ptwXY_interpolation interpolationXY
Definition: MCGIDI.h:306
int MCGIDI_sampling_sampleX_from_pdfOfX(MCGIDI_pdfOfX *dist, MCGIDI_pdfsOfXGivenW_sampled *sampled, double r)
statusMessageReporting * smr
Definition: MCGIDI.h:312
int MCGIDI_sampling_interpolationValues(statusMessageReporting *smr, ptwXY_interpolation interpolation, double *ws, double y1, double y2, double *y)
ptwXY_interpolation interpolationWY
Definition: MCGIDI.h:306
xDataTOM_Int MCGIDI_misc_binarySearch(xDataTOM_Int n, double *ds, double d)
Definition: MCGIDI_misc.cc:228
enum ptwXY_interpolation_e ptwXY_interpolation
int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW(MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double r)
double(* rng)(void *)
Definition: MCGIDI.h:258
ptwXY_interpolation interpolationXY
Definition: MCGIDI.h:313
double getProjectileEnergy(void) const
Definition: MCGIDI.h:97
MCGIDI_pdfOfX * dist
Definition: MCGIDI.h:308

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_sampling_interpolationValues ( statusMessageReporting smr,
ptwXY_interpolation  interpolation,
double *  ws,
double  y1,
double  y2,
double *  y 
)

Definition at line 187 of file MCGIDI_sampling.cc.

187  {
188 
189  double frac;
190 
191  if( interpolation == ptwXY_interpolationLinLin ) {
192  frac = ( ws[2] - ws[0] ) / ( ws[2] - ws[1] );
193  *y = frac * y1 + ( 1 - frac ) * y2; }
194  else if( interpolation == ptwXY_interpolationLogLin ) {
195  frac = G4Log( ws[2] / ws[0] ) / G4Log( ws[2] / ws[1] );
196  *y = frac * y1 + ( 1 - frac ) * y2; }
197  else if( interpolation == ptwXY_interpolationLinLog ) {
198  frac = ( ws[2] - ws[0] ) / ( ws[2] - ws[1] );
199  *y = y1 * G4Pow::GetInstance()->powA( y2 / y1, frac ); }
200  else if( interpolation == ptwXY_interpolationLogLog ) {
201  frac = G4Log( ws[2] / ws[0] ) / G4Log( ws[2] / ws[1] );
202  *y = y2 * G4Pow::GetInstance()->powA( y2 / y1, frac ); }
203  else { // This should never happen.
204  smr_setReportError2( smr, smr_unknownID, 1, "bad interpolation = %d\n", interpolation );
205  return( 1 );
206  }
207  return( 0 );
208 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
#define smr_setReportError2(smr, libraryID, code, fmt,...)
#define smr_unknownID
G4double G4Log(G4double x)
Definition: G4Log.hh:230

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_sampling_pdfsOfX_release ( statusMessageReporting ,
MCGIDI_pdfOfX dist 
)

Definition at line 41 of file MCGIDI_sampling.cc.

41  {
42 
43  smr_freeMemory( (void **) &(dist->Xs) );
44  return( 0 );
45 }
double * Xs
Definition: MCGIDI.h:299
void * smr_freeMemory(void **p)

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_sampling_pdfsOfXGivenW_initialize ( statusMessageReporting ,
MCGIDI_pdfsOfXGivenW dists 
)

Definition at line 20 of file MCGIDI_sampling.cc.

20  {
21 
22  memset( dists, 0, sizeof( MCGIDI_pdfsOfXGivenW ) );
23  return( 0 );
24 }

Here is the caller graph for this function:

int MCGIDI_sampling_pdfsOfXGivenW_release ( statusMessageReporting smr,
MCGIDI_pdfsOfXGivenW dists 
)

Definition at line 28 of file MCGIDI_sampling.cc.

28  {
29 
30  int i;
31 
32  for( i = 0; i < dists->numberOfWs; i++ ) MCGIDI_sampling_pdfsOfX_release( smr, &(dists->dist[i]) );
33  smr_freeMemory( (void **) &(dists->Ws) );
34  smr_freeMemory( (void **) &(dists->dist) );
36  return( 0 );
37 }
int MCGIDI_sampling_pdfsOfXGivenW_initialize(statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *dists)
void * smr_freeMemory(void **p)
int MCGIDI_sampling_pdfsOfX_release(statusMessageReporting *smr, MCGIDI_pdfOfX *dist)
MCGIDI_pdfOfX * dist
Definition: MCGIDI.h:308

Here is the call graph for this function:

Here is the caller graph for this function:

double MCGIDI_sampling_ptwXY_getValueAtX ( ptwXYPoints ptwXY,
double  x1 
)

Definition at line 212 of file MCGIDI_sampling.cc.

212  {
213 
214  double y1;
215 
216  if( ptwXY_getValueAtX( ptwXY, x1, &y1 ) == nfu_XOutsideDomain ) {
217  if( x1 < ptwXY_getXMin( ptwXY ) ) {
218  ptwXY_getValueAtX( ptwXY, ptwXY_getXMin( ptwXY ), &y1 ); }
219  else {
220  ptwXY_getValueAtX( ptwXY, ptwXY_getXMax( ptwXY ), &y1 );
221  }
222  }
223  return( y1 );
224 }
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
Definition: ptwXY_core.cc:844
double ptwXY_getXMin(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1206
double ptwXY_getXMax(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:1239

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_sampling_sampleX_from_pdfOfX ( MCGIDI_pdfOfX dist,
MCGIDI_pdfsOfXGivenW_sampled sampled,
double  rngValue 
)

Definition at line 98 of file MCGIDI_sampling.cc.

98  {
99 
100  int iX;
101  double d1, d2, frac;
102 
103  iX = sampled->iX1 = MCGIDI_misc_binarySearch( dist->numberOfXs, dist->cdf, rngValue );
104 
105  if( iX < 0 ) { /* This should never happen. */
106  smr_setReportError2( sampled->smr, smr_unknownID, 1, "bad iX = %d\n", iX );
107  sampled->x = dist->Xs[0];
108  return( 1 );
109  }
110  if( sampled->interpolationXY == ptwXY_interpolationFlat ) {
111  frac = ( dist->cdf[iX+1] - rngValue ) / ( dist->cdf[iX+1] - dist->cdf[iX] );
112  sampled->x = frac * dist->Xs[iX] + ( 1 - frac ) * dist->Xs[iX+1]; }
113  else {
114  double s1 = dist->pdf[iX+1] - dist->pdf[iX];
115 
116  if( s1 == 0. ) {
117  if( dist->pdf[iX] == 0 ) {
118  sampled->x = dist->Xs[iX];
119  if( iX == 0 ) sampled->x = dist->Xs[1]; }
120  else {
121  frac = ( dist->cdf[iX+1] - rngValue ) / ( dist->cdf[iX+1] - dist->cdf[iX] );
122  sampled->x = frac * dist->Xs[iX] + ( 1 - frac ) * dist->Xs[iX+1];
123  } }
124  else {
125  s1 = s1 / ( dist->Xs[iX+1] - dist->Xs[iX] );
126  d1 = rngValue - dist->cdf[iX];
127  d2 = dist->cdf[iX+1] - rngValue;
128  if( d2 > d1 ) { /* Closer to iX. */
129  sampled->x = dist->Xs[iX] + ( std::sqrt( dist->pdf[iX] * dist->pdf[iX] + 2. * s1 * d1 ) - dist->pdf[iX] ) / s1; }
130  else { /* Closer to iX + 1. */
131  sampled->x = dist->Xs[iX+1] - ( dist->pdf[iX+1] - std::sqrt( dist->pdf[iX+1] * dist->pdf[iX+1] - 2. * s1 * d2 ) ) / s1;
132  }
133  }
134  }
135 
136  return( 0 );
137 }
double * cdf
Definition: MCGIDI.h:301
static const G4double d2
double * Xs
Definition: MCGIDI.h:299
statusMessageReporting * smr
Definition: MCGIDI.h:312
#define smr_setReportError2(smr, libraryID, code, fmt,...)
#define smr_unknownID
static const G4double d1
xDataTOM_Int MCGIDI_misc_binarySearch(xDataTOM_Int n, double *ds, double d)
Definition: MCGIDI_misc.cc:228
int numberOfXs
Definition: MCGIDI.h:298
double * pdf
Definition: MCGIDI.h:300
ptwXY_interpolation interpolationXY
Definition: MCGIDI.h:313

Here is the call graph for this function:

Here is the caller graph for this function:

int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW ( MCGIDI_pdfsOfXGivenW dists,
MCGIDI_pdfsOfXGivenW_sampled sampled,
double  rngValue 
)

Definition at line 49 of file MCGIDI_sampling.cc.

49  {
50 
51  int iW, iX1;
52 
53  sampled->interpolationWY = dists->interpolationWY;
54  sampled->interpolationXY = dists->interpolationXY;
55  iW = sampled->iW = MCGIDI_misc_binarySearch( dists->numberOfWs, dists->Ws, sampled->w );
56  sampled->frac = 1;
57 
58  if( iW == -2 ) { /* w < first value of Ws. */
59  return( MCGIDI_sampling_sampleX_from_pdfOfX( dists->dist, sampled, rngValue ) ); }
60  else if( iW == -1 ) { /* w > last value of Ws. */
61  return( MCGIDI_sampling_sampleX_from_pdfOfX( &(dists->dist[dists->numberOfWs-1]), sampled, rngValue ) ); }
62  else {
63  if( MCGIDI_sampling_sampleX_from_pdfOfX( &(dists->dist[iW]), sampled, rngValue ) ) return( 1 );
64  if( dists->interpolationWY != ptwXY_interpolationFlat ) { // ptwXY_interpolationOther was not allowed at startup.
65  double xSampled = sampled->x, frac = 1.;
66 
67  iX1 = sampled->iX1;
68  if( MCGIDI_sampling_sampleX_from_pdfOfX( &(dists->dist[iW+1]), sampled, rngValue ) ) return( 1 );
69 
71  frac = ( dists->Ws[iW+1] - sampled->w ) / ( dists->Ws[iW+1] - dists->Ws[iW] );
72  sampled->x = frac * xSampled + ( 1 - frac ) * sampled->x; }
73  else if( dists->interpolationWY == ptwXY_interpolationLogLin ) {
74  frac = G4Log( dists->Ws[iW+1] / sampled->w ) / G4Log( dists->Ws[iW+1] / dists->Ws[iW] );
75  sampled->x = frac * xSampled + ( 1 - frac ) * sampled->x; }
76  else if( dists->interpolationWY == ptwXY_interpolationLinLog ) {
77  frac = ( dists->Ws[iW+1] - sampled->w ) / ( dists->Ws[iW+1] - dists->Ws[iW] );
78  sampled->x = xSampled * G4Pow::GetInstance()->powA( sampled->x / xSampled, frac ); }
79  else if( dists->interpolationWY == ptwXY_interpolationLogLog ) {
80  frac = G4Log( dists->Ws[iW+1] / sampled->w ) / G4Log( dists->Ws[iW+1] / dists->Ws[iW] );
81  sampled->x = xSampled * G4Pow::GetInstance()->powA( sampled->x / xSampled, frac ); }
82  else { // This should never happen.
83  smr_setReportError2( sampled->smr, smr_unknownID, 1, "bad interpolation = %d\n", dists->interpolationWY );
84  return( 1 );
85  }
86 
87  sampled->iX2 = sampled->iX1;
88  sampled->iX1 = iX1;
89  sampled->frac = frac;
90  }
91  }
92 
93  return( 0 );
94 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
ptwXY_interpolation interpolationXY
Definition: MCGIDI.h:306
ptwXY_interpolation interpolationWY
Definition: MCGIDI.h:313
int MCGIDI_sampling_sampleX_from_pdfOfX(MCGIDI_pdfOfX *dist, MCGIDI_pdfsOfXGivenW_sampled *sampled, double r)
statusMessageReporting * smr
Definition: MCGIDI.h:312
#define smr_setReportError2(smr, libraryID, code, fmt,...)
ptwXY_interpolation interpolationWY
Definition: MCGIDI.h:306
#define smr_unknownID
xDataTOM_Int MCGIDI_misc_binarySearch(xDataTOM_Int n, double *ds, double d)
Definition: MCGIDI_misc.cc:228
G4double G4Log(G4double x)
Definition: G4Log.hh:230
ptwXY_interpolation interpolationXY
Definition: MCGIDI.h:313
MCGIDI_pdfOfX * dist
Definition: MCGIDI.h:308

Here is the call graph for this function:

Here is the caller graph for this function: