| Geant4 Cross Reference |
1 /*
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 #include <stdlib.h>
6 #include <cmath>
7
8 #include "MCGIDI.h"
9
10 #if defined __cplusplus
11 #include "G4Log.hh"
12 #include "G4Pow.hh"
13 namespace GIDI {
14 using namespace GIDI;
15 #endif
16
17 /*
18 ************************************************************
19 */
20 int MCGIDI_sampling_pdfsOfXGivenW_initialize( statusMessageReporting * /*smr*/, MCGIDI_pdfsOfXGivenW *dists ) {
21
22 memset( dists, 0, sizeof( MCGIDI_pdfsOfXGivenW ) );
23 return( 0 );
24 }
25 /*
26 ************************************************************
27 */
28 int MCGIDI_sampling_pdfsOfXGivenW_release( statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *dists ) {
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) );
35 MCGIDI_sampling_pdfsOfXGivenW_initialize( smr, dists );
36 return( 0 );
37 }
38 /*
39 ************************************************************
40 */
41 int MCGIDI_sampling_pdfsOfX_release( statusMessageReporting * /*smr*/, MCGIDI_pdfOfX *dist ) {
42
43 smr_freeMemory( (void **) &(dist->Xs) );
44 return( 0 );
45 }
46 /*
47 ************************************************************
48 */
49 int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue ) {
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
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] );
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 }
95 /*
96 ************************************************************
97 */
98 int MCGIDI_sampling_sampleX_from_pdfOfX( MCGIDI_pdfOfX *dist, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue ) {
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 }
138 /*
139 ************************************************************
140 */
141 int MCGIDI_sampling_doubleDistribution( statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *pdfOfWGivenV, MCGIDI_pdfsOfXGivenW *pdfOfXGivenVAndW,
142 MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
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 }
184 /*
185 ************************************************************
186 */
187 int MCGIDI_sampling_interpolationValues( statusMessageReporting *smr, ptwXY_interpolation interpolation, double *ws, double y1, double y2, double *y ) {
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 }
209 /*
210 ************************************************************
211 */
212 double MCGIDI_sampling_ptwXY_getValueAtX( ptwXYPoints *ptwXY, double x1 ) {
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 }
225
226 #if defined __cplusplus
227 }
228 #endif
229
230