10 #if defined __cplusplus
22 memset( dists, 0,
sizeof( MCGIDI_pdfsOfXGivenW ) );
53 sampled->interpolationWY = dists->interpolationWY;
54 sampled->interpolationXY = dists->interpolationXY;
64 if( dists->interpolationWY != ptwXY_interpolationFlat ) {
65 double xSampled = sampled->x, frac = 1.;
70 if( dists->interpolationWY == ptwXY_interpolationLinLin ) {
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] );
79 else if( dists->interpolationWY == ptwXY_interpolationLogLog ) {
80 frac =
G4Log( dists->Ws[iW+1] / sampled->w ) /
G4Log( dists->Ws[iW+1] / dists->Ws[iW] );
83 smr_setReportError2( sampled->smr, smr_unknownID, 1,
"bad interpolation = %d\n", dists->interpolationWY );
87 sampled->iX2 = sampled->iX1;
106 smr_setReportError2( sampled->smr, smr_unknownID, 1,
"bad iX = %d\n", iX );
107 sampled->x = dist->Xs[0];
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]; }
114 double s1 = dist->pdf[iX+1] - dist->pdf[iX];
117 if( dist->pdf[iX] == 0 ) {
118 sampled->x = dist->Xs[iX];
119 if( iX == 0 ) sampled->x = dist->Xs[1]; }
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];
125 s1 = s1 / ( dist->Xs[iX+1] - dist->Xs[iX] );
126 d1 = rngValue - dist->cdf[iX];
127 d2 = dist->cdf[iX+1] - rngValue;
129 sampled->x = dist->Xs[iX] + ( std::sqrt( dist->pdf[iX] * dist->pdf[iX] + 2. * s1 * d1 ) - dist->pdf[iX] ) / s1; }
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;
142 MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
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;
152 sampledW.interpolationXY = pdfOfWGivenV->interpolationXY;
155 interpolationWY = ptwXY_interpolationFlat;
159 iV = pdfOfWGivenV->numberOfWs - 1;
161 e_in = pdfOfWGivenV->Ws[iV];
165 sampledX.w = sampledW.x;
168 if( interpolationWY != ptwXY_interpolationFlat ) {
169 double x = sampledX.x,
w = sampledW.x, Vs[3] = { e_in, pdfOfWGivenV->Ws[iV], pdfOfWGivenV->Ws[iV+1] };
172 sampledX.w = sampledW.x;
179 decaySamplingInfo->mu = sampledW.x;
180 decaySamplingInfo->Ep = sampledX.x;
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] );
200 else if( interpolation == ptwXY_interpolationLogLog ) {
201 frac =
G4Log( ws[2] / ws[0] ) /
G4Log( ws[2] / ws[1] );
204 smr_setReportError2( smr, smr_unknownID, 1,
"bad interpolation = %d\n", interpolation );
226 #if defined __cplusplus
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
xDataTOM_Int MCGIDI_misc_binarySearch(xDataTOM_Int n, double *ds, double d)
const G4double w[NPOINTSGL]
int MCGIDI_sampling_pdfsOfX_release(statusMessageReporting *, MCGIDI_pdfOfX *dist)
int MCGIDI_sampling_pdfsOfXGivenW_release(statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *dists)
double ptwXY_getXMax(ptwXYPoints *ptwXY)
int MCGIDI_sampling_doubleDistribution(statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *pdfOfWGivenV, MCGIDI_pdfsOfXGivenW *pdfOfXGivenVAndW, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
void * smr_freeMemory(void **p)
int MCGIDI_sampling_interpolationValues(statusMessageReporting *smr, ptwXY_interpolation interpolation, double *ws, double y1, double y2, double *y)
G4double G4Log(G4double x)
int MCGIDI_sampling_pdfsOfXGivenW_initialize(statusMessageReporting *, MCGIDI_pdfsOfXGivenW *dists)
nfu_status ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x, double *y)
const G4double x[NPOINTSGL]
int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW(MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue)
double MCGIDI_sampling_ptwXY_getValueAtX(ptwXYPoints *ptwXY, double x1)
double ptwXY_getXMin(ptwXYPoints *ptwXY)
int MCGIDI_sampling_sampleX_from_pdfOfX(MCGIDI_pdfOfX *dist, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue)