| Geant4 Cross Reference |
1 /* 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(
21
22 memset( dists, 0, sizeof( MCGIDI_pdfsOfXGi
23 return( 0 );
24 }
25 /*
26 **********************************************
27 */
28 int MCGIDI_sampling_pdfsOfXGivenW_release( sta
29
30 int i;
31
32 for( i = 0; i < dists->numberOfWs; i++ ) M
33 smr_freeMemory( (void **) &(dists->Ws) );
34 smr_freeMemory( (void **) &(dists->dist) )
35 MCGIDI_sampling_pdfsOfXGivenW_initialize(
36 return( 0 );
37 }
38 /*
39 **********************************************
40 */
41 int MCGIDI_sampling_pdfsOfX_release( statusMes
42
43 smr_freeMemory( (void **) &(dist->Xs) );
44 return( 0 );
45 }
46 /*
47 **********************************************
48 */
49 int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW
50
51 int iW, iX1;
52
53 sampled->interpolationWY = dists->interpol
54 sampled->interpolationXY = dists->interpol
55 iW = sampled->iW = MCGIDI_misc_binarySearc
56 sampled->frac = 1;
57
58 if( iW == -2 ) { /* w < first v
59 return( MCGIDI_sampling_sampleX_from_p
60 else if( iW == -1 ) { /* w > last va
61 return( MCGIDI_sampling_sampleX_from_p
62 else {
63 if( MCGIDI_sampling_sampleX_from_pdfOf
64 if( dists->interpolationWY != ptwXY_in
65 double xSampled = sampled->x, frac
66
67 iX1 = sampled->iX1;
68 if( MCGIDI_sampling_sampleX_from_p
69
70 if( dists->interpolationWY == ptwX
71 frac = ( dists->Ws[iW+1] - sam
72 sampled->x = frac * xSampled +
73 else if( dists->interpolationWY ==
74 frac = G4Log( dists->Ws[iW+1]
75 sampled->x = frac * xSampled +
76 else if( dists->interpolationWY ==
77 frac = ( dists->Ws[iW+1] - sam
78 sampled->x = xSampled * G4Pow:
79 else if( dists->interpolationWY ==
80 frac = G4Log( dists->Ws[iW+1]
81 sampled->x = xSampled * G4Pow:
82 else { // This should never happe
83 smr_setReportError2( sampled->
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( MCGID
99
100 int iX;
101 double d1, d2, frac;
102
103 iX = sampled->iX1 = MCGIDI_misc_binarySear
104
105 if( iX < 0 ) { /* This sho
106 smr_setReportError2( sampled->smr, smr
107 sampled->x = dist->Xs[0];
108 return( 1 );
109 }
110 if( sampled->interpolationXY == ptwXY_inte
111 frac = ( dist->cdf[iX+1] - rngValue )
112 sampled->x = frac * dist->Xs[iX] + ( 1
113 else {
114 double s1 = dist->pdf[iX+1] - dist->pd
115
116 if( s1 == 0. ) {
117 if( dist->pdf[iX] == 0 ) {
118 sampled->x = dist->Xs[iX];
119 if( iX == 0 ) sampled->x = dis
120 else {
121 frac = ( dist->cdf[iX+1] - rng
122 sampled->x = frac * dist->Xs[i
123 } }
124 else {
125 s1 = s1 / ( dist->Xs[iX+1] - dist-
126 d1 = rngValue - dist->cdf[iX];
127 d2 = dist->cdf[iX+1] - rngValue;
128 if( d2 > d1 ) { /* Closer to i
129 sampled->x = dist->Xs[iX] + (
130 else { /* Closer to i
131 sampled->x = dist->Xs[iX+1] -
132 }
133 }
134 }
135
136 return( 0 );
137 }
138 /*
139 **********************************************
140 */
141 int MCGIDI_sampling_doubleDistribution( status
142 MCGIDI_quantitiesLookupModes &modes, M
143
144 int iV;
145 double e_in = modes.getProjectileEnergy( )
146 double randomW = decaySamplingInfo->rng( d
147 MCGIDI_pdfsOfXGivenW_sampled sampledX, sam
148 ptwXY_interpolation interpolationWY = pdfO
149
150 sampledX.smr = smr;
151 sampledW.smr = smr;
152 sampledW.interpolationXY = pdfOfWGivenV->i
153 iV = MCGIDI_misc_binarySearch( pdfOfWGiven
154 if( iV < 0 ) {
155 interpolationWY = ptwXY_interpolationF
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( &(pdf
165 sampledX.w = sampledW.x;
166 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW
167
168 if( interpolationWY != ptwXY_interpolation
169 double x = sampledX.x, w = sampledW.x,
170
171 MCGIDI_sampling_sampleX_from_pdfOfX( &
172 sampledX.w = sampledW.x;
173 MCGIDI_sampling_sampleX_from_pdfsOfXGi
174
175 MCGIDI_sampling_interpolationValues( s
176 MCGIDI_sampling_interpolationValues( s
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( statu
188
189 double frac;
190
191 if( interpolation == ptwXY_interpolationLi
192 frac = ( ws[2] - ws[0] ) / ( ws[2] - w
193 *y = frac * y1 + ( 1 - frac ) * y2; }
194 else if( interpolation == ptwXY_interpolat
195 frac = G4Log( ws[2] / ws[0] ) / G4Log(
196 *y = frac * y1 + ( 1 - frac ) * y2; }
197 else if( interpolation == ptwXY_interpolat
198 frac = ( ws[2] - ws[0] ) / ( ws[2] - w
199 *y = y1 * G4Pow::GetInstance()->powA(
200 else if( interpolation == ptwXY_interpolat
201 frac = G4Log( ws[2] / ws[0] ) / G4Log(
202 *y = y2 * G4Pow::GetInstance()->powA(
203 else { // This should never happen.
204 smr_setReportError2( smr, smr_unknownI
205 return( 1 );
206 }
207 return( 0 );
208 }
209 /*
210 **********************************************
211 */
212 double MCGIDI_sampling_ptwXY_getValueAtX( ptwX
213
214 double y1;
215
216 if( ptwXY_getValueAtX( ptwXY, x1, &y1 ) ==
217 if( x1 < ptwXY_getXMin( ptwXY ) ) {
218 ptwXY_getValueAtX( ptwXY, ptwXY_ge
219 else {
220 ptwXY_getValueAtX( ptwXY, ptwXY_ge
221 }
222 }
223 return( y1 );
224 }
225
226 #if defined __cplusplus
227 }
228 #endif
229
230