| Geant4 Cross Reference |
1 /* 1 /*
2 # <<BEGIN-copyright>> 2 # <<BEGIN-copyright>>
3 # <<END-copyright>> 3 # <<END-copyright>>
4 */ 4 */
5 #include <string.h> 5 #include <string.h>
6 #include <cmath> 6 #include <cmath>
7 7
8 #if defined __cplusplus 8 #if defined __cplusplus
9 #include "G4Exp.hh" 9 #include "G4Exp.hh"
10 #include "G4Log.hh" 10 #include "G4Log.hh"
11 #include "G4Pow.hh" 11 #include "G4Pow.hh"
12 namespace GIDI { 12 namespace GIDI {
13 using namespace GIDI; 13 using namespace GIDI;
14 #endif 14 #endif
15 const double C1 = 0.04, C2 = 1.8e-6/*, C3 = 6. 15 const double C1 = 0.04, C2 = 1.8e-6/*, C3 = 6.7e-7*/;
16 /* 16 /*
17 const double Et1 = 130., Et3 = 41.; 17 const double Et1 = 130., Et3 = 41.;
18 */ 18 */
19 #if defined __cplusplus 19 #if defined __cplusplus
20 } 20 }
21 #endif 21 #endif
22 22
23 #include "MCGIDI_fromTOM.h" 23 #include "MCGIDI_fromTOM.h"
24 #include "MCGIDI_misc.h" 24 #include "MCGIDI_misc.h"
25 #include "MCGIDI_private.h" 25 #include "MCGIDI_private.h"
26 26
27 #if defined __cplusplus 27 #if defined __cplusplus
28 namespace GIDI { 28 namespace GIDI {
29 using namespace GIDI; 29 using namespace GIDI;
30 #endif 30 #endif
31 31
32 32
33 static int MCGIDI_KalbachMann_parseFromTOM2( s 33 static int MCGIDI_KalbachMann_parseFromTOM2( statusMessageReporting *smr, int dataPerEout, int index, xDataTOM_KalbachMannCoefficients *coefficientsXData,
34 double energyInFactor, double energyOu 34 double energyInFactor, double energyOutFactor, MCGIDI_KalbachMann *KalbachMann );
35 static double MCGIDI_KalbachMann_S_a_or_b( dou 35 static double MCGIDI_KalbachMann_S_a_or_b( double Z_AB, double N_AB, double Z_C, double N_C, double I_ab );
36 /* 36 /*
37 ********************************************** 37 ************************************************************
38 */ 38 */
39 MCGIDI_KalbachMann *MCGIDI_KalbachMann_new( st 39 MCGIDI_KalbachMann *MCGIDI_KalbachMann_new( statusMessageReporting *smr, ptwXY_interpolation interpolationWY,
40 ptwXY_interpolation interpolationXY ) 40 ptwXY_interpolation interpolationXY ) {
41 41
42 MCGIDI_KalbachMann *KalbachMann; 42 MCGIDI_KalbachMann *KalbachMann;
43 43
44 if( ( KalbachMann = (MCGIDI_KalbachMann *) 44 if( ( KalbachMann = (MCGIDI_KalbachMann *) smr_malloc2( smr, sizeof( MCGIDI_KalbachMann ), 0, "KalbachMann" ) ) == NULL ) return( NULL );
45 if( MCGIDI_KalbachMann_initialize( smr, Ka 45 if( MCGIDI_KalbachMann_initialize( smr, KalbachMann, interpolationWY, interpolationXY ) ) KalbachMann = MCGIDI_KalbachMann_free( smr, KalbachMann );
46 return( KalbachMann ); 46 return( KalbachMann );
47 } 47 }
48 /* 48 /*
49 ********************************************** 49 ************************************************************
50 */ 50 */
51 int MCGIDI_KalbachMann_initialize( statusMessa 51 int MCGIDI_KalbachMann_initialize( statusMessageReporting * /*smr*/, MCGIDI_KalbachMann *KalbachMann, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY ) {
52 52
53 memset( KalbachMann, 0, sizeof( MCGIDI_Kal 53 memset( KalbachMann, 0, sizeof( MCGIDI_KalbachMann ) );
54 KalbachMann->dists.interpolationWY = inter 54 KalbachMann->dists.interpolationWY = interpolationWY;
55 KalbachMann->dists.interpolationXY = inter 55 KalbachMann->dists.interpolationXY = interpolationXY;
56 return( 0 ); 56 return( 0 );
57 } 57 }
58 /* 58 /*
59 ********************************************** 59 ************************************************************
60 */ 60 */
61 MCGIDI_KalbachMann *MCGIDI_KalbachMann_free( s 61 MCGIDI_KalbachMann *MCGIDI_KalbachMann_free( statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann ) {
62 62
63 MCGIDI_KalbachMann_release( smr, KalbachMa 63 MCGIDI_KalbachMann_release( smr, KalbachMann );
64 smr_freeMemory( (void **) &KalbachMann ); 64 smr_freeMemory( (void **) &KalbachMann );
65 return( NULL ); 65 return( NULL );
66 } 66 }
67 /* 67 /*
68 ********************************************** 68 ************************************************************
69 */ 69 */
70 int MCGIDI_KalbachMann_release( statusMessageR 70 int MCGIDI_KalbachMann_release( statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann ) {
71 71
72 int i; 72 int i;
73 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMan 73 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMann->dists);
74 74
75 for( i = 0; i < dists->numberOfWs; i++ ) { 75 for( i = 0; i < dists->numberOfWs; i++ ) {
76 smr_freeMemory( (void **) &(KalbachMan 76 smr_freeMemory( (void **) &(KalbachMann->ras[i].rs) );
77 smr_freeMemory( (void **) &(dists->dis 77 smr_freeMemory( (void **) &(dists->dist[i].Xs) );
78 } 78 }
79 smr_freeMemory( (void **) &(KalbachMann->r 79 smr_freeMemory( (void **) &(KalbachMann->ras) );
80 smr_freeMemory( (void **) &(dists->Ws) ); 80 smr_freeMemory( (void **) &(dists->Ws) );
81 smr_freeMemory( (void **) &(dists->dist) ) 81 smr_freeMemory( (void **) &(dists->dist) );
82 82
83 MCGIDI_KalbachMann_initialize( smr, Kalbac 83 MCGIDI_KalbachMann_initialize( smr, KalbachMann, ptwXY_interpolationLinLin, ptwXY_interpolationLinLin );
84 return( 0 ); 84 return( 0 );
85 } 85 }
86 /* 86 /*
87 ********************************************** 87 ************************************************************
88 */ 88 */
89 int MCGIDI_KalbachMann_parseFromTOM( statusMes 89 int MCGIDI_KalbachMann_parseFromTOM( statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution ) {
90 90
91 MCGIDI_KalbachMann *KalbachMann = NULL; 91 MCGIDI_KalbachMann *KalbachMann = NULL;
92 xDataTOM_element *KalbachMannElement; 92 xDataTOM_element *KalbachMannElement;
93 int index, dataPerEout = 3; 93 int index, dataPerEout = 3;
94 double energyInFactor, energyOutFactor; 94 double energyInFactor, energyOutFactor;
95 xDataTOM_xDataInfo *xDataInfo; 95 xDataTOM_xDataInfo *xDataInfo;
96 xDataTOM_KalbachMann *KalbachMannXData; 96 xDataTOM_KalbachMann *KalbachMannXData;
97 ptwXY_interpolation interpolationXY, inter 97 ptwXY_interpolation interpolationXY, interpolationWY;
98 char const *energyFromUnit, *energyToUnit 98 char const *energyFromUnit, *energyToUnit = "MeV";
99 99
100 MCGIDI_POP *productPOP = distribution->pro 100 MCGIDI_POP *productPOP = distribution->product->pop;
101 double productZ = productPOP->Z, productA 101 double productZ = productPOP->Z, productA = productPOP->A, productN = productA - productZ;
102 MCGIDI_target_heated *targetHeated = MCGID 102 MCGIDI_target_heated *targetHeated = MCGIDI_product_getTargetHeated( smr, distribution->product );
103 MCGIDI_POP *projectilePOP = MCGIDI_target_ 103 MCGIDI_POP *projectilePOP = MCGIDI_target_heated_getPOPForProjectile( smr, targetHeated );
104 double projectileZ = projectilePOP->Z, pro 104 double projectileZ = projectilePOP->Z, projectileA = projectilePOP->A, projectileN = projectileA - projectileZ;
105 MCGIDI_POP *targetPOP = MCGIDI_target_heat 105 MCGIDI_POP *targetPOP = MCGIDI_target_heated_getPOPForTarget( smr, targetHeated );
106 double targetZ = targetPOP->Z, targetA = t 106 double targetZ = targetPOP->Z, targetA = targetPOP->A, targetN = targetA - targetZ;
107 double Ia = 0., Ib = 0., Ma = -1, mb = -1; 107 double Ia = 0., Ib = 0., Ma = -1, mb = -1;
108 108
109 if( ( targetA == 0 ) && ( targetZ == 6 ) ) 109 if( ( targetA == 0 ) && ( targetZ == 6 ) ) { /* Special case for C_000 evaluation. */
110 targetN = 6; 110 targetN = 6;
111 targetA = 12; 111 targetA = 12;
112 } 112 }
113 if( ( KalbachMannElement = xDataTOME_getOn 113 if( ( KalbachMannElement = xDataTOME_getOneElementByName( smr, element, "KalbachMann", 1 ) ) == NULL ) goto err;
114 114
115 if( MCGIDI_fromTOM_interpolation( smr, Kal 115 if( MCGIDI_fromTOM_interpolation( smr, KalbachMannElement, 0, &interpolationWY ) ) goto err;
116 if( MCGIDI_fromTOM_interpolation( smr, Kal 116 if( MCGIDI_fromTOM_interpolation( smr, KalbachMannElement, 1, &interpolationXY ) ) goto err;
117 117
118 xDataInfo = &(KalbachMannElement->xDataInf 118 xDataInfo = &(KalbachMannElement->xDataInfo);
119 KalbachMannXData = (xDataTOM_KalbachMann * 119 KalbachMannXData = (xDataTOM_KalbachMann *) xDataInfo->data;
120 if( KalbachMannXData->type == xDataTOM_Kal 120 if( KalbachMannXData->type == xDataTOM_KalbachMannType_fra ) dataPerEout = 4;
121 121
122 energyFromUnit = xDataTOM_axes_getUnit( sm 122 energyFromUnit = xDataTOM_axes_getUnit( smr, &(xDataInfo->axes), 0 );
123 if( !smr_isOk( smr ) ) goto err; 123 if( !smr_isOk( smr ) ) goto err;
124 energyInFactor = MCGIDI_misc_getUnitConver 124 energyInFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyFromUnit, energyToUnit );
125 if( !smr_isOk( smr ) ) goto err; 125 if( !smr_isOk( smr ) ) goto err;
126 126
127 energyFromUnit = xDataTOM_axes_getUnit( sm 127 energyFromUnit = xDataTOM_axes_getUnit( smr, &(xDataInfo->axes), 1 );
128 if( !smr_isOk( smr ) ) goto err; 128 if( !smr_isOk( smr ) ) goto err;
129 energyOutFactor = MCGIDI_misc_getUnitConve 129 energyOutFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyFromUnit, energyToUnit );
130 if( !smr_isOk( smr ) ) goto err; 130 if( !smr_isOk( smr ) ) goto err;
131 131
132 if( ( KalbachMann = distribution->KalbachM 132 if( ( KalbachMann = distribution->KalbachMann = MCGIDI_KalbachMann_new( smr, interpolationWY, interpolationXY ) ) == NULL ) goto err;
133 133
134 /* 134 /*
135 double productMass MCGIDI_product_getMass_ 135 double productMass MCGIDI_product_getMass_MeV( smr, distribution->product ), residualMass;
136 */ 136 */
137 KalbachMann->energyToMeVFactor = MCGIDI_mi 137 KalbachMann->energyToMeVFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyToUnit, "MeV" );
138 KalbachMann->massFactor = (double) product 138 KalbachMann->massFactor = (double) productZ + productN; /* This is not correct as masses are needed not Z and N. */
139 KalbachMann->massFactor /= projectileN + p 139 KalbachMann->massFactor /= projectileN + projectileZ + targetZ + targetN - productZ + productN;
140 KalbachMann->massFactor += 1.; 140 KalbachMann->massFactor += 1.;
141 141
142 if( projectileZ == 0 ) { 142 if( projectileZ == 0 ) {
143 if( projectileN == 1 ) Ma = 1; } 143 if( projectileN == 1 ) Ma = 1; }
144 else if( projectileZ == 1 ) { 144 else if( projectileZ == 1 ) {
145 if( projectileN == 1 ) { 145 if( projectileN == 1 ) {
146 Ma = 1; } 146 Ma = 1; }
147 else if( projectileN == 2 ) { 147 else if( projectileN == 2 ) {
148 Ia = 2.22; 148 Ia = 2.22;
149 Ma = 1; } } 149 Ma = 1; } }
150 else if( projectileZ == 2 ) { 150 else if( projectileZ == 2 ) {
151 if( projectileN == 2 ) { 151 if( projectileN == 2 ) {
152 Ia = 28.3; 152 Ia = 28.3;
153 Ma = 0; 153 Ma = 0;
154 } 154 }
155 } 155 }
156 156
157 if( productZ == 0 ) { 157 if( productZ == 0 ) {
158 if( productN == 1 ) mb = 0.5; } 158 if( productN == 1 ) mb = 0.5; }
159 else if( productZ == 1 ) { 159 else if( productZ == 1 ) {
160 if( productN == 1 ) { 160 if( productN == 1 ) {
161 mb = 1; } 161 mb = 1; }
162 else if( productN == 2 ) { 162 else if( productN == 2 ) {
163 Ia = 2.22; 163 Ia = 2.22;
164 mb = 1; } 164 mb = 1; }
165 else if( productN == 3 ) { 165 else if( productN == 3 ) {
166 Ib = 8.48; 166 Ib = 8.48;
167 mb = 1; } } 167 mb = 1; } }
168 else if( productZ == 2 ) { 168 else if( productZ == 2 ) {
169 if( productN == 1 ) { 169 if( productN == 1 ) {
170 Ib = 7.72; 170 Ib = 7.72;
171 mb = 1; } 171 mb = 1; }
172 else if( productN == 2 ) { 172 else if( productN == 2 ) {
173 Ib = 28.3; 173 Ib = 28.3;
174 mb = 2; 174 mb = 2;
175 } 175 }
176 } 176 }
177 177
178 KalbachMann->Ma = Ma; 178 KalbachMann->Ma = Ma;
179 KalbachMann->mb = mb; 179 KalbachMann->mb = mb;
180 180
181 KalbachMann->Sa = MCGIDI_KalbachMann_S_a_o 181 KalbachMann->Sa = MCGIDI_KalbachMann_S_a_or_b( targetZ, targetN, targetZ + projectileZ, targetN + projectileN, Ia );
182 KalbachMann->Sb = MCGIDI_KalbachMann_S_a_o 182 KalbachMann->Sb = MCGIDI_KalbachMann_S_a_or_b( projectileZ + targetZ - productZ, projectileN + targetN - productN,
183 targetZ + projectileZ, targetN + proje 183 targetZ + projectileZ, targetN + projectileN, Ib );
184 184
185 KalbachMann->dists.numberOfWs = 0; 185 KalbachMann->dists.numberOfWs = 0;
186 if( ( KalbachMann->dists.Ws = (double *) s 186 if( ( KalbachMann->dists.Ws = (double *) smr_malloc2( smr, KalbachMannXData->numberOfEnergies * sizeof( double ), 0, "KalbachMann->dists->Ws" ) ) == NULL ) goto err;
187 if( ( KalbachMann->dists.dist = (MCGIDI_pd 187 if( ( KalbachMann->dists.dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, KalbachMannXData->numberOfEnergies * sizeof( MCGIDI_pdfOfX ), 0, "KalbachMann->dists->dist" ) ) == NULL ) goto err;
188 if( ( KalbachMann->ras = (MCGIDI_KalbachMa 188 if( ( KalbachMann->ras = (MCGIDI_KalbachMann_ras *) smr_malloc2( smr, KalbachMannXData->numberOfEnergies * sizeof( MCGIDI_KalbachMann_ras ), 0, "KalbachMann->ras" ) ) == NULL ) goto err;
189 189
190 for( index = 0; index < KalbachMannXData-> 190 for( index = 0; index < KalbachMannXData->numberOfEnergies; index++ ) {
191 if( MCGIDI_KalbachMann_parseFromTOM2( 191 if( MCGIDI_KalbachMann_parseFromTOM2( smr, dataPerEout, index, &(KalbachMannXData->coefficients[index]),
192 energyInFactor, energyOutFactor, K 192 energyInFactor, energyOutFactor, KalbachMann ) ) goto err;
193 } 193 }
194 194
195 if( ( KalbachMann->frame = MCGIDI_misc_get 195 if( ( KalbachMann->frame = MCGIDI_misc_getProductFrame( smr, KalbachMannElement ) ) == xDataTOM_frame_invalid ) goto err;
196 distribution->type = MCGIDI_distributionTy 196 distribution->type = MCGIDI_distributionType_KalbachMann_e;
197 197
198 return( 0 ); 198 return( 0 );
199 199
200 err: 200 err:
201 if( KalbachMann != NULL ) MCGIDI_KalbachMa 201 if( KalbachMann != NULL ) MCGIDI_KalbachMann_free( smr, KalbachMann );
202 return( 1 ); 202 return( 1 );
203 } 203 }
204 /* 204 /*
205 ********************************************** 205 ************************************************************
206 */ 206 */
207 static int MCGIDI_KalbachMann_parseFromTOM2( s 207 static int MCGIDI_KalbachMann_parseFromTOM2( statusMessageReporting *smr, int dataPerEout, int index, xDataTOM_KalbachMannCoefficients *coefficientsXData,
208 double energyInFactor, double energyOu 208 double energyInFactor, double energyOutFactor, MCGIDI_KalbachMann *KalbachMann ) {
209 209
210 int i, j, n = coefficientsXData->length / 210 int i, j, n = coefficientsXData->length / dataPerEout;
211 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMan 211 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMann->dists);
212 MCGIDI_pdfOfX *dist = &(dists->dist[index] 212 MCGIDI_pdfOfX *dist = &(dists->dist[index]);
213 double norm, *p, *rs = NULL, *as_ = NULL, 213 double norm, *p, *rs = NULL, *as_ = NULL, *Xs = NULL, *pdf, *cdf;
214 nfu_status status; 214 nfu_status status;
215 ptwXYPoints *pdfXY = NULL; 215 ptwXYPoints *pdfXY = NULL;
216 ptwXYPoint *point; 216 ptwXYPoint *point;
217 ptwXPoints *cdfX = NULL; 217 ptwXPoints *cdfX = NULL;
218 char const *ptwFunc = ""; 218 char const *ptwFunc = "";
219 219
220 if( ( Xs = (double *) smr_malloc2( smr, 3 220 if( ( Xs = (double *) smr_malloc2( smr, 3 * n * sizeof( double ), 0, "Xs" ) ) == NULL ) goto err;
221 pdf = &(Xs[n]); 221 pdf = &(Xs[n]);
222 cdf = &(pdf[n]); 222 cdf = &(pdf[n]);
223 223
224 if( ( rs = (double *) smr_malloc2( smr, ( 224 if( ( rs = (double *) smr_malloc2( smr, ( dataPerEout - 2 ) * n * sizeof( double ), 0, "rs" ) ) == NULL ) goto err;
225 if( dataPerEout == 4 ) as_ = &(rs[n]); 225 if( dataPerEout == 4 ) as_ = &(rs[n]);
226 226
227 ptwFunc = "ptwXY_new"; 227 ptwFunc = "ptwXY_new";
228 if( ( pdfXY = ptwXY_new( KalbachMann->dist 228 if( ( pdfXY = ptwXY_new( KalbachMann->dists.interpolationXY, NULL, 2., 1e-3, n, 10, &status, 0 ) ) == NULL ) goto errXY;
229 229
230 ptwFunc = "ptwXY_setXYPairAtIndex"; 230 ptwFunc = "ptwXY_setXYPairAtIndex";
231 for( i = 0, p = coefficientsXData->coeffic 231 for( i = 0, p = coefficientsXData->coefficients; i < n; i++, p += dataPerEout ) {
232 if( ( status = ptwXY_setValueAtX( pdfX 232 if( ( status = ptwXY_setValueAtX( pdfXY, p[0], p[1] ) ) != nfu_Okay ) goto errXY;
233 rs[i] = p[2]; 233 rs[i] = p[2];
234 if( dataPerEout == 4 ) as_[i] = p[3]; 234 if( dataPerEout == 4 ) as_[i] = p[3];
235 } 235 }
236 236
237 for( j = 0; j < n; j++ ) { 237 for( j = 0; j < n; j++ ) {
238 point = ptwXY_getPointAtIndex_Unsafely 238 point = ptwXY_getPointAtIndex_Unsafely( pdfXY, j );
239 Xs[j] = energyOutFactor * point->x; 239 Xs[j] = energyOutFactor * point->x;
240 pdf[j] = point->y / energyOutFactor; 240 pdf[j] = point->y / energyOutFactor;
241 } 241 }
242 242
243 ptwFunc = "ptwXY_runningIntegral"; 243 ptwFunc = "ptwXY_runningIntegral";
244 if( ( cdfX = ptwXY_runningIntegral( pdfXY, 244 if( ( cdfX = ptwXY_runningIntegral( pdfXY, &status ) ) == NULL ) goto errXY;
245 norm = ptwX_getPointAtIndex_Unsafely( cdfX 245 norm = ptwX_getPointAtIndex_Unsafely( cdfX, n - 1 );
246 if( std::fabs( 1. - norm ) > 0.99 ) { 246 if( std::fabs( 1. - norm ) > 0.99 ) {
247 smr_setReportError2( smr, smr_unknownI 247 smr_setReportError2( smr, smr_unknownID, 1, "bad norm = %e for angular.linear data", norm );
248 goto err; 248 goto err;
249 } 249 }
250 for( j = 0; j < n; j++ ) cdf[j] = ptwX_get 250 for( j = 0; j < n; j++ ) cdf[j] = ptwX_getPointAtIndex_Unsafely( cdfX, j ) / norm;
251 for( j = 0; j < n; j++ ) pdf[j] /= norm; 251 for( j = 0; j < n; j++ ) pdf[j] /= norm;
252 252
253 dists->numberOfWs++; 253 dists->numberOfWs++;
254 dists->Ws[index] = energyInFactor * coeffi 254 dists->Ws[index] = energyInFactor * coefficientsXData->value;
255 dist->numberOfXs = n; 255 dist->numberOfXs = n;
256 dist->Xs = Xs; 256 dist->Xs = Xs;
257 dist->pdf = pdf; 257 dist->pdf = pdf;
258 dist->cdf = cdf; 258 dist->cdf = cdf;
259 KalbachMann->ras[index].rs = rs; 259 KalbachMann->ras[index].rs = rs;
260 KalbachMann->ras[index].as = as_; 260 KalbachMann->ras[index].as = as_;
261 261
262 pdfXY = ptwXY_free( pdfXY ); 262 pdfXY = ptwXY_free( pdfXY );
263 cdfX = ptwX_free( cdfX ); 263 cdfX = ptwX_free( cdfX );
264 return( 0 ); 264 return( 0 );
265 265
266 errXY: 266 errXY:
267 smr_setReportError2( smr, smr_unknownID, 1 267 smr_setReportError2( smr, smr_unknownID, 1, "%s error = %d: %s\n", ptwFunc, status, nfu_statusMessage( status ) );
268 268
269 err: 269 err:
270 if( Xs != NULL ) smr_freeMemory( (void **) 270 if( Xs != NULL ) smr_freeMemory( (void **) &Xs);
271 if( rs != NULL ) smr_freeMemory( (void **) 271 if( rs != NULL ) smr_freeMemory( (void **) &rs);
272 if( pdfXY != NULL ) ptwXY_free( pdfXY ); 272 if( pdfXY != NULL ) ptwXY_free( pdfXY );
273 if( cdfX != NULL ) cdfX = ptwX_free( cdfX 273 if( cdfX != NULL ) cdfX = ptwX_free( cdfX );
274 return( 1 ); 274 return( 1 );
275 } 275 }
276 /* 276 /*
277 ********************************************** 277 ************************************************************
278 */ 278 */
279 static double MCGIDI_KalbachMann_S_a_or_b( dou 279 static double MCGIDI_KalbachMann_S_a_or_b( double Z_AB, double N_AB, double Z_C, double N_C, double I_ab ) {
280 280
281 double A_AB = Z_AB + N_AB, A_C = Z_C + N_C 281 double A_AB = Z_AB + N_AB, A_C = Z_C + N_C;
282 double invA_AB_third = 1.0 / G4Pow::GetIns 282 double invA_AB_third = 1.0 / G4Pow::GetInstance()->A13( A_AB ), invA_C_third = 1.0 / G4Pow::GetInstance()->A13 ( A_C );
283 double NZA_AB = ( N_AB - Z_AB ) * ( N_AB - 283 double NZA_AB = ( N_AB - Z_AB ) * ( N_AB - Z_AB ) / A_AB, NZA_C = ( N_C - Z_C ) * ( N_C - Z_C ) / A_C, S;
284 284
285 S = 15.68 * ( A_C - A_AB ) 285 S = 15.68 * ( A_C - A_AB ) - 28.07 * ( NZA_C - NZA_AB )
286 - 18.56 * ( A_C * invA_C_third - A_AB 286 - 18.56 * ( A_C * invA_C_third - A_AB * invA_AB_third ) + 33.22 * ( NZA_C * invA_C_third - NZA_AB * invA_AB_third )
287 - 0.717 * ( Z_C * Z_C * invA_C_third - 287 - 0.717 * ( Z_C * Z_C * invA_C_third - Z_AB * Z_AB * invA_AB_third ) + 1.211 * ( Z_C * Z_C / A_C - Z_AB * Z_AB / A_AB )
288 - I_ab; 288 - I_ab;
289 return( S ); 289 return( S );
290 } 290 }
291 /* 291 /*
292 ********************************************** 292 ************************************************************
293 */ 293 */
294 int MCGIDI_KalbachMann_sampleEp( statusMessage 294 int MCGIDI_KalbachMann_sampleEp( statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, MCGIDI_quantitiesLookupModes &modes,
295 MCGIDI_decaySamplingInfo *decaySamplin 295 MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
296 296
297 double Epl, Epu, Ep, r, r2, rl, ru, a, a2, 297 double Epl, Epu, Ep, r, r2, rl, ru, a, a2, al, au, mu, randomEp = decaySamplingInfo->rng( decaySamplingInfo->rngState );
298 MCGIDI_pdfsOfXGivenW_sampled sampled; 298 MCGIDI_pdfsOfXGivenW_sampled sampled;
299 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMan 299 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMann->dists);
300 ptwXY_interpolation interpolationWY; 300 ptwXY_interpolation interpolationWY;
301 301
302 sampled.smr = smr; 302 sampled.smr = smr;
303 sampled.w = modes.getProjectileEnergy( ); 303 sampled.w = modes.getProjectileEnergy( );
304 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW 304 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( dists, &sampled, randomEp );
305 305
306 interpolationWY = sampled.interpolationWY; 306 interpolationWY = sampled.interpolationWY;
307 if( sampled.iW < 0 ) { 307 if( sampled.iW < 0 ) {
308 interpolationWY = ptwXY_interpolationF 308 interpolationWY = ptwXY_interpolationFlat;
309 if( sampled.iW == -2 ) { /* ??? 309 if( sampled.iW == -2 ) { /* ???????????? This should probably report a warning. */
310 sampled.iW = 0; } 310 sampled.iW = 0; }
311 else if( sampled.iW == -1 ) { 311 else if( sampled.iW == -1 ) {
312 sampled.iW = dists->numberOfWs - 1 312 sampled.iW = dists->numberOfWs - 1;
313 } 313 }
314 } 314 }
315 315
316 Ep = sampled.x; 316 Ep = sampled.x; /* Sampled Ep. */
317 if( sampled.interpolationXY == ptwXY_inter 317 if( sampled.interpolationXY == ptwXY_interpolationFlat ) { /* Now sample r. */
318 r = KalbachMann->ras[sampled.iW].rs[sa 318 r = KalbachMann->ras[sampled.iW].rs[sampled.iX1]; }
319 else { 319 else {
320 Epl = dists->dist[sampled.iW].Xs[sampl 320 Epl = dists->dist[sampled.iW].Xs[sampled.iX1];
321 Epu = dists->dist[sampled.iW].Xs[sampl 321 Epu = dists->dist[sampled.iW].Xs[sampled.iX1+1];
322 rl = KalbachMann->ras[sampled.iW].rs[s 322 rl = KalbachMann->ras[sampled.iW].rs[sampled.iX1];
323 ru = KalbachMann->ras[sampled.iW].rs[s 323 ru = KalbachMann->ras[sampled.iW].rs[sampled.iX1+1];
324 r = ( ru - rl ) / ( Epu - Epl ) * ( Ep 324 r = ( ru - rl ) / ( Epu - Epl ) * ( Ep - Epl ) + rl;
325 } 325 }
326 if( interpolationWY == ptwXY_interpolation 326 if( interpolationWY == ptwXY_interpolationLinLin ) {
327 if( sampled.interpolationXY == ptwXY_i 327 if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
328 r2 = KalbachMann->ras[sampled.iW+1 328 r2 = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2]; }
329 else { 329 else {
330 Epl = dists->dist[sampled.iW+1].Xs 330 Epl = dists->dist[sampled.iW+1].Xs[sampled.iX2];
331 Epu = dists->dist[sampled.iW+1].Xs 331 Epu = dists->dist[sampled.iW+1].Xs[sampled.iX2+1];
332 rl = KalbachMann->ras[sampled.iW+1 332 rl = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2];
333 ru = KalbachMann->ras[sampled.iW+1 333 ru = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2+1];
334 r2 = ( ru - rl ) / ( Epu - Epl ) * 334 r2 = ( ru - rl ) / ( Epu - Epl ) * ( Ep - Epl ) + rl;
335 } 335 }
336 r = sampled.frac * r + ( 1. - sampled. 336 r = sampled.frac * r + ( 1. - sampled.frac ) * r2;
337 } 337 }
338 338
339 if( KalbachMann->ras[0].as == NULL ) { 339 if( KalbachMann->ras[0].as == NULL ) { /* Now determine a. */
340 double X1, X3_2; 340 double X1, X3_2;
341 double eb = KalbachMann->massFactor * 341 double eb = KalbachMann->massFactor * KalbachMann->energyToMeVFactor * Ep + KalbachMann->Sb;
342 342
343 X1 = eb; /* Not valid for 343 X1 = eb; /* Not valid for ea > Et1. */
344 X3_2 = eb * eb; /* Not valid for 344 X3_2 = eb * eb; /* Not valid for ea > Et3. */
345 a = X1 * ( C1 + C2 * X1 * X1 ) + C2 * 345 a = X1 * ( C1 + C2 * X1 * X1 ) + C2 * KalbachMann->Ma * KalbachMann->mb * X3_2 * X3_2; }
346 else { 346 else {
347 if( sampled.interpolationXY == ptwXY_i 347 if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
348 a = KalbachMann->ras[sampled.iW].a 348 a = KalbachMann->ras[sampled.iW].as[sampled.iX1]; }
349 else { 349 else {
350 Epl = dists->dist[sampled.iW].Xs[s 350 Epl = dists->dist[sampled.iW].Xs[sampled.iX1];
351 Epu = dists->dist[sampled.iW].Xs[s 351 Epu = dists->dist[sampled.iW].Xs[sampled.iX1+1];
352 al = KalbachMann->ras[sampled.iW]. 352 al = KalbachMann->ras[sampled.iW].as[sampled.iX1];
353 au = KalbachMann->ras[sampled.iW]. 353 au = KalbachMann->ras[sampled.iW].as[sampled.iX1+1];
354 a = ( au - al ) / ( Epu - Epl ) * 354 a = ( au - al ) / ( Epu - Epl ) * ( Ep - Epl ) + al;
355 } 355 }
356 a2 = 0.; 356 a2 = 0.;
357 if( interpolationWY == ptwXY_interpola 357 if( interpolationWY == ptwXY_interpolationLinLin ) {
358 if( sampled.interpolationXY == ptw 358 if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
359 a2 = KalbachMann->ras[sampled. 359 a2 = KalbachMann->ras[sampled.iW+1].as[sampled.iX2]; }
360 else { 360 else {
361 Epl = dists->dist[sampled.iW+1 361 Epl = dists->dist[sampled.iW+1].Xs[sampled.iX2];
362 Epu = dists->dist[sampled.iW+1 362 Epu = dists->dist[sampled.iW+1].Xs[sampled.iX2+1];
363 al = KalbachMann->ras[sampled. 363 al = KalbachMann->ras[sampled.iW+1].as[sampled.iX2];
364 au = KalbachMann->ras[sampled. 364 au = KalbachMann->ras[sampled.iW+1].as[sampled.iX2+1];
365 a2 = ( au - al ) / ( Epu - Epl 365 a2 = ( au - al ) / ( Epu - Epl ) * ( Ep - Epl ) + al;
366 } 366 }
367 } 367 }
368 a = sampled.frac * a + ( 1. - sampled. 368 a = sampled.frac * a + ( 1. - sampled.frac ) * a2;
369 } 369 }
370 370
371 /* In the following: Cosh[ a mu ] + r 371 /* In the following: Cosh[ a mu ] + r Sinh[ a mu ] = ( 1 - r ) Cosh[ a mu ] + r ( Cosh[ a mu ] + Sinh[ a mu ] ). */
372 if( decaySamplingInfo->rng( decaySamplingI 372 if( decaySamplingInfo->rng( decaySamplingInfo->rngState ) >= r ) { /* Sample the '( 1 - r ) Cosh[ a mu ]' term. */
373 double T = ( 2. * decaySamplingInfo->r 373 double T = ( 2. * decaySamplingInfo->rng( decaySamplingInfo->rngState ) - 1. ) * std::sinh( a );
374 374
375 mu = G4Log( T + std::sqrt( T * T + 1. 375 mu = G4Log( T + std::sqrt( T * T + 1. ) ) / a; }
376 else { 376 else { /* Sample the 'r ( Cosh[ a mu ] + Sinh[ a mu ]' term. */
377 double rng1 = decaySamplingInfo->rng( 377 double rng1 = decaySamplingInfo->rng( decaySamplingInfo->rngState ), exp_a = G4Exp( a );
378 378
379 mu = G4Log( rng1 * exp_a + ( 1. - rng1 379 mu = G4Log( rng1 * exp_a + ( 1. - rng1 ) / exp_a ) / a;
380 } 380 }
381 if( mu < -1 ) { 381 if( mu < -1 ) {
382 mu = -1;} 382 mu = -1;}
383 else if( mu > 1 ) { 383 else if( mu > 1 ) {
384 mu = 1; 384 mu = 1;
385 } 385 }
386 386
387 decaySamplingInfo->frame = KalbachMann->fr 387 decaySamplingInfo->frame = KalbachMann->frame;
388 decaySamplingInfo->Ep = Ep; 388 decaySamplingInfo->Ep = Ep;
389 decaySamplingInfo->mu = mu; 389 decaySamplingInfo->mu = mu;
390 return( !smr_isOk( smr ) ); 390 return( !smr_isOk( smr ) );
391 } 391 }
392 392
393 #if defined __cplusplus 393 #if defined __cplusplus
394 } 394 }
395 #endif 395 #endif
396 396