Geant4 Cross Reference |
1 /* 1 /* 2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>> 3 # <<END-copyright>> 3 # <<END-copyright>> 4 */ 4 */ 5 #include <stdlib.h> 5 #include <stdlib.h> 6 #include <cmath> 6 #include <cmath> 7 7 8 #include "MCGIDI.h" 8 #include "MCGIDI.h" 9 9 10 #if defined __cplusplus 10 #if defined __cplusplus 11 #include "G4Log.hh" 11 #include "G4Log.hh" 12 #include "G4Pow.hh" 12 #include "G4Pow.hh" 13 namespace GIDI { 13 namespace GIDI { 14 using namespace GIDI; 14 using namespace GIDI; 15 #endif 15 #endif 16 16 17 /* 17 /* 18 ********************************************** 18 ************************************************************ 19 */ 19 */ 20 int MCGIDI_sampling_pdfsOfXGivenW_initialize( 20 int MCGIDI_sampling_pdfsOfXGivenW_initialize( statusMessageReporting * /*smr*/, MCGIDI_pdfsOfXGivenW *dists ) { 21 21 22 memset( dists, 0, sizeof( MCGIDI_pdfsOfXGi 22 memset( dists, 0, sizeof( MCGIDI_pdfsOfXGivenW ) ); 23 return( 0 ); 23 return( 0 ); 24 } 24 } 25 /* 25 /* 26 ********************************************** 26 ************************************************************ 27 */ 27 */ 28 int MCGIDI_sampling_pdfsOfXGivenW_release( sta 28 int MCGIDI_sampling_pdfsOfXGivenW_release( statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *dists ) { 29 29 30 int i; 30 int i; 31 31 32 for( i = 0; i < dists->numberOfWs; i++ ) M 32 for( i = 0; i < dists->numberOfWs; i++ ) MCGIDI_sampling_pdfsOfX_release( smr, &(dists->dist[i]) ); 33 smr_freeMemory( (void **) &(dists->Ws) ); 33 smr_freeMemory( (void **) &(dists->Ws) ); 34 smr_freeMemory( (void **) &(dists->dist) ) 34 smr_freeMemory( (void **) &(dists->dist) ); 35 MCGIDI_sampling_pdfsOfXGivenW_initialize( 35 MCGIDI_sampling_pdfsOfXGivenW_initialize( smr, dists ); 36 return( 0 ); 36 return( 0 ); 37 } 37 } 38 /* 38 /* 39 ********************************************** 39 ************************************************************ 40 */ 40 */ 41 int MCGIDI_sampling_pdfsOfX_release( statusMes 41 int MCGIDI_sampling_pdfsOfX_release( statusMessageReporting * /*smr*/, MCGIDI_pdfOfX *dist ) { 42 42 43 smr_freeMemory( (void **) &(dist->Xs) ); 43 smr_freeMemory( (void **) &(dist->Xs) ); 44 return( 0 ); 44 return( 0 ); 45 } 45 } 46 /* 46 /* 47 ********************************************** 47 ************************************************************ 48 */ 48 */ 49 int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW 49 int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue ) { 50 50 51 int iW, iX1; 51 int iW, iX1; 52 52 53 sampled->interpolationWY = dists->interpol 53 sampled->interpolationWY = dists->interpolationWY; 54 sampled->interpolationXY = dists->interpol 54 sampled->interpolationXY = dists->interpolationXY; 55 iW = sampled->iW = MCGIDI_misc_binarySearc 55 iW = sampled->iW = MCGIDI_misc_binarySearch( dists->numberOfWs, dists->Ws, sampled->w ); 56 sampled->frac = 1; 56 sampled->frac = 1; 57 57 58 if( iW == -2 ) { /* w < first v 58 if( iW == -2 ) { /* w < first value of Ws. */ 59 return( MCGIDI_sampling_sampleX_from_p 59 return( MCGIDI_sampling_sampleX_from_pdfOfX( dists->dist, sampled, rngValue ) ); } 60 else if( iW == -1 ) { /* w > last va 60 else if( iW == -1 ) { /* w > last value of Ws. */ 61 return( MCGIDI_sampling_sampleX_from_p 61 return( MCGIDI_sampling_sampleX_from_pdfOfX( &(dists->dist[dists->numberOfWs-1]), sampled, rngValue ) ); } 62 else { 62 else { 63 if( MCGIDI_sampling_sampleX_from_pdfOf 63 if( MCGIDI_sampling_sampleX_from_pdfOfX( &(dists->dist[iW]), sampled, rngValue ) ) return( 1 ); 64 if( dists->interpolationWY != ptwXY_in 64 if( dists->interpolationWY != ptwXY_interpolationFlat ) { // ptwXY_interpolationOther was not allowed at startup. 65 double xSampled = sampled->x, frac 65 double xSampled = sampled->x, frac = 1.; 66 66 67 iX1 = sampled->iX1; 67 iX1 = sampled->iX1; 68 if( MCGIDI_sampling_sampleX_from_p 68 if( MCGIDI_sampling_sampleX_from_pdfOfX( &(dists->dist[iW+1]), sampled, rngValue ) ) return( 1 ); 69 69 70 if( dists->interpolationWY == ptwX 70 if( dists->interpolationWY == ptwXY_interpolationLinLin ) { 71 frac = ( dists->Ws[iW+1] - sam 71 frac = ( dists->Ws[iW+1] - sampled->w ) / ( dists->Ws[iW+1] - dists->Ws[iW] ); 72 sampled->x = frac * xSampled + 72 sampled->x = frac * xSampled + ( 1 - frac ) * sampled->x; } 73 else if( dists->interpolationWY == 73 else if( dists->interpolationWY == ptwXY_interpolationLogLin ) { 74 frac = G4Log( dists->Ws[iW+1] 74 frac = G4Log( dists->Ws[iW+1] / sampled->w ) / G4Log( dists->Ws[iW+1] / dists->Ws[iW] ); 75 sampled->x = frac * xSampled + 75 sampled->x = frac * xSampled + ( 1 - frac ) * sampled->x; } 76 else if( dists->interpolationWY == 76 else if( dists->interpolationWY == ptwXY_interpolationLinLog ) { 77 frac = ( dists->Ws[iW+1] - sam 77 frac = ( dists->Ws[iW+1] - sampled->w ) / ( dists->Ws[iW+1] - dists->Ws[iW] ); 78 sampled->x = xSampled * G4Pow: 78 sampled->x = xSampled * G4Pow::GetInstance()->powA( sampled->x / xSampled, frac ); } 79 else if( dists->interpolationWY == 79 else if( dists->interpolationWY == ptwXY_interpolationLogLog ) { 80 frac = G4Log( dists->Ws[iW+1] 80 frac = G4Log( dists->Ws[iW+1] / sampled->w ) / G4Log( dists->Ws[iW+1] / dists->Ws[iW] ); 81 sampled->x = xSampled * G4Pow: 81 sampled->x = xSampled * G4Pow::GetInstance()->powA( sampled->x / xSampled, frac ); } 82 else { // This should never happe 82 else { // This should never happen. 83 smr_setReportError2( sampled-> 83 smr_setReportError2( sampled->smr, smr_unknownID, 1, "bad interpolation = %d\n", dists->interpolationWY ); 84 return( 1 ); 84 return( 1 ); 85 } 85 } 86 86 87 sampled->iX2 = sampled->iX1; 87 sampled->iX2 = sampled->iX1; 88 sampled->iX1 = iX1; 88 sampled->iX1 = iX1; 89 sampled->frac = frac; 89 sampled->frac = frac; 90 } 90 } 91 } 91 } 92 92 93 return( 0 ); 93 return( 0 ); 94 } 94 } 95 /* 95 /* 96 ********************************************** 96 ************************************************************ 97 */ 97 */ 98 int MCGIDI_sampling_sampleX_from_pdfOfX( MCGID 98 int MCGIDI_sampling_sampleX_from_pdfOfX( MCGIDI_pdfOfX *dist, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue ) { 99 99 100 int iX; 100 int iX; 101 double d1, d2, frac; 101 double d1, d2, frac; 102 102 103 iX = sampled->iX1 = MCGIDI_misc_binarySear 103 iX = sampled->iX1 = MCGIDI_misc_binarySearch( dist->numberOfXs, dist->cdf, rngValue ); 104 104 105 if( iX < 0 ) { /* This sho 105 if( iX < 0 ) { /* This should never happen. */ 106 smr_setReportError2( sampled->smr, smr 106 smr_setReportError2( sampled->smr, smr_unknownID, 1, "bad iX = %d\n", iX ); 107 sampled->x = dist->Xs[0]; 107 sampled->x = dist->Xs[0]; 108 return( 1 ); 108 return( 1 ); 109 } 109 } 110 if( sampled->interpolationXY == ptwXY_inte 110 if( sampled->interpolationXY == ptwXY_interpolationFlat ) { 111 frac = ( dist->cdf[iX+1] - rngValue ) 111 frac = ( dist->cdf[iX+1] - rngValue ) / ( dist->cdf[iX+1] - dist->cdf[iX] ); 112 sampled->x = frac * dist->Xs[iX] + ( 1 112 sampled->x = frac * dist->Xs[iX] + ( 1 - frac ) * dist->Xs[iX+1]; } 113 else { 113 else { 114 double s1 = dist->pdf[iX+1] - dist->pd 114 double s1 = dist->pdf[iX+1] - dist->pdf[iX]; 115 115 116 if( s1 == 0. ) { 116 if( s1 == 0. ) { 117 if( dist->pdf[iX] == 0 ) { 117 if( dist->pdf[iX] == 0 ) { 118 sampled->x = dist->Xs[iX]; 118 sampled->x = dist->Xs[iX]; 119 if( iX == 0 ) sampled->x = dis 119 if( iX == 0 ) sampled->x = dist->Xs[1]; } 120 else { 120 else { 121 frac = ( dist->cdf[iX+1] - rng 121 frac = ( dist->cdf[iX+1] - rngValue ) / ( dist->cdf[iX+1] - dist->cdf[iX] ); 122 sampled->x = frac * dist->Xs[i 122 sampled->x = frac * dist->Xs[iX] + ( 1 - frac ) * dist->Xs[iX+1]; 123 } } 123 } } 124 else { 124 else { 125 s1 = s1 / ( dist->Xs[iX+1] - dist- 125 s1 = s1 / ( dist->Xs[iX+1] - dist->Xs[iX] ); 126 d1 = rngValue - dist->cdf[iX]; 126 d1 = rngValue - dist->cdf[iX]; 127 d2 = dist->cdf[iX+1] - rngValue; 127 d2 = dist->cdf[iX+1] - rngValue; 128 if( d2 > d1 ) { /* Closer to i 128 if( d2 > d1 ) { /* Closer to iX. */ 129 sampled->x = dist->Xs[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 i 130 else { /* Closer to iX + 1. */ 131 sampled->x = dist->Xs[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 } 132 } 133 } 133 } 134 } 134 } 135 135 136 return( 0 ); 136 return( 0 ); 137 } 137 } 138 /* 138 /* 139 ********************************************** 139 ************************************************************ 140 */ 140 */ 141 int MCGIDI_sampling_doubleDistribution( status 141 int MCGIDI_sampling_doubleDistribution( statusMessageReporting *smr, MCGIDI_pdfsOfXGivenW *pdfOfWGivenV, MCGIDI_pdfsOfXGivenW *pdfOfXGivenVAndW, 142 MCGIDI_quantitiesLookupModes &modes, M 142 MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo ) { 143 143 144 int iV; 144 int iV; 145 double e_in = modes.getProjectileEnergy( ) 145 double e_in = modes.getProjectileEnergy( ); 146 double randomW = decaySamplingInfo->rng( d 146 double randomW = decaySamplingInfo->rng( decaySamplingInfo->rngState ), randomX = decaySamplingInfo->rng( decaySamplingInfo->rngState ); 147 MCGIDI_pdfsOfXGivenW_sampled sampledX, sam 147 MCGIDI_pdfsOfXGivenW_sampled sampledX, sampledW; 148 ptwXY_interpolation interpolationWY = pdfO 148 ptwXY_interpolation interpolationWY = pdfOfWGivenV->interpolationWY; 149 149 150 sampledX.smr = smr; 150 sampledX.smr = smr; 151 sampledW.smr = smr; 151 sampledW.smr = smr; 152 sampledW.interpolationXY = pdfOfWGivenV->i 152 sampledW.interpolationXY = pdfOfWGivenV->interpolationXY; 153 iV = MCGIDI_misc_binarySearch( pdfOfWGiven 153 iV = MCGIDI_misc_binarySearch( pdfOfWGivenV->numberOfWs, pdfOfWGivenV->Ws, e_in ); 154 if( iV < 0 ) { 154 if( iV < 0 ) { 155 interpolationWY = ptwXY_interpolationF 155 interpolationWY = ptwXY_interpolationFlat; 156 if( iV == -2 ) { 156 if( iV == -2 ) { 157 iV = 0; } 157 iV = 0; } 158 else { 158 else { 159 iV = pdfOfWGivenV->numberOfWs - 1; 159 iV = pdfOfWGivenV->numberOfWs - 1; 160 } 160 } 161 e_in = pdfOfWGivenV->Ws[iV]; 161 e_in = pdfOfWGivenV->Ws[iV]; 162 } 162 } 163 163 164 MCGIDI_sampling_sampleX_from_pdfOfX( &(pdf 164 MCGIDI_sampling_sampleX_from_pdfOfX( &(pdfOfWGivenV->dist[iV]), &sampledW, randomW ); 165 sampledX.w = sampledW.x; 165 sampledX.w = sampledW.x; 166 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW 166 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(pdfOfXGivenVAndW[iV]), &sampledX, randomX ); 167 167 168 if( interpolationWY != ptwXY_interpolation 168 if( interpolationWY != ptwXY_interpolationFlat ) { 169 double x = sampledX.x, w = sampledW.x, 169 double x = sampledX.x, w = sampledW.x, Vs[3] = { e_in, pdfOfWGivenV->Ws[iV], pdfOfWGivenV->Ws[iV+1] }; 170 170 171 MCGIDI_sampling_sampleX_from_pdfOfX( & 171 MCGIDI_sampling_sampleX_from_pdfOfX( &(pdfOfWGivenV->dist[iV+1]), &sampledW, randomW ); 172 sampledX.w = sampledW.x; 172 sampledX.w = sampledW.x; 173 MCGIDI_sampling_sampleX_from_pdfsOfXGi 173 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( &(pdfOfXGivenVAndW[iV+1]), &sampledX, randomX ); 174 174 175 MCGIDI_sampling_interpolationValues( s 175 MCGIDI_sampling_interpolationValues( smr, interpolationWY, Vs, w, sampledW.x, &sampledW.x ); 176 MCGIDI_sampling_interpolationValues( s 176 MCGIDI_sampling_interpolationValues( smr, interpolationWY, Vs, x, sampledX.x, &sampledX.x ); 177 } 177 } 178 178 179 decaySamplingInfo->mu = sampledW.x; 179 decaySamplingInfo->mu = sampledW.x; 180 decaySamplingInfo->Ep = sampledX.x; 180 decaySamplingInfo->Ep = sampledX.x; 181 181 182 return( 0 ); 182 return( 0 ); 183 } 183 } 184 /* 184 /* 185 ********************************************** 185 ************************************************************ 186 */ 186 */ 187 int MCGIDI_sampling_interpolationValues( statu 187 int MCGIDI_sampling_interpolationValues( statusMessageReporting *smr, ptwXY_interpolation interpolation, double *ws, double y1, double y2, double *y ) { 188 188 189 double frac; 189 double frac; 190 190 191 if( interpolation == ptwXY_interpolationLi 191 if( interpolation == ptwXY_interpolationLinLin ) { 192 frac = ( ws[2] - ws[0] ) / ( ws[2] - w 192 frac = ( ws[2] - ws[0] ) / ( ws[2] - ws[1] ); 193 *y = frac * y1 + ( 1 - frac ) * y2; } 193 *y = frac * y1 + ( 1 - frac ) * y2; } 194 else if( interpolation == ptwXY_interpolat 194 else if( interpolation == ptwXY_interpolationLogLin ) { 195 frac = G4Log( ws[2] / ws[0] ) / G4Log( 195 frac = G4Log( ws[2] / ws[0] ) / G4Log( ws[2] / ws[1] ); 196 *y = frac * y1 + ( 1 - frac ) * y2; } 196 *y = frac * y1 + ( 1 - frac ) * y2; } 197 else if( interpolation == ptwXY_interpolat 197 else if( interpolation == ptwXY_interpolationLinLog ) { 198 frac = ( ws[2] - ws[0] ) / ( ws[2] - w 198 frac = ( ws[2] - ws[0] ) / ( ws[2] - ws[1] ); 199 *y = y1 * G4Pow::GetInstance()->powA( 199 *y = y1 * G4Pow::GetInstance()->powA( y2 / y1, frac ); } 200 else if( interpolation == ptwXY_interpolat 200 else if( interpolation == ptwXY_interpolationLogLog ) { 201 frac = G4Log( ws[2] / ws[0] ) / G4Log( 201 frac = G4Log( ws[2] / ws[0] ) / G4Log( ws[2] / ws[1] ); 202 *y = y2 * G4Pow::GetInstance()->powA( 202 *y = y2 * G4Pow::GetInstance()->powA( y2 / y1, frac ); } 203 else { // This should never happen. 203 else { // This should never happen. 204 smr_setReportError2( smr, smr_unknownI 204 smr_setReportError2( smr, smr_unknownID, 1, "bad interpolation = %d\n", interpolation ); 205 return( 1 ); 205 return( 1 ); 206 } 206 } 207 return( 0 ); 207 return( 0 ); 208 } 208 } 209 /* 209 /* 210 ********************************************** 210 ************************************************************ 211 */ 211 */ 212 double MCGIDI_sampling_ptwXY_getValueAtX( ptwX 212 double MCGIDI_sampling_ptwXY_getValueAtX( ptwXYPoints *ptwXY, double x1 ) { 213 213 214 double y1; 214 double y1; 215 215 216 if( ptwXY_getValueAtX( ptwXY, x1, &y1 ) == 216 if( ptwXY_getValueAtX( ptwXY, x1, &y1 ) == nfu_XOutsideDomain ) { 217 if( x1 < ptwXY_getXMin( ptwXY ) ) { 217 if( x1 < ptwXY_getXMin( ptwXY ) ) { 218 ptwXY_getValueAtX( ptwXY, ptwXY_ge 218 ptwXY_getValueAtX( ptwXY, ptwXY_getXMin( ptwXY ), &y1 ); } 219 else { 219 else { 220 ptwXY_getValueAtX( ptwXY, ptwXY_ge 220 ptwXY_getValueAtX( ptwXY, ptwXY_getXMax( ptwXY ), &y1 ); 221 } 221 } 222 } 222 } 223 return( y1 ); 223 return( y1 ); 224 } 224 } 225 225 226 #if defined __cplusplus 226 #if defined __cplusplus 227 } 227 } 228 #endif 228 #endif 229 229 230 230