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