| Geant4 Cross Reference |
1 /* 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.
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( s
34 double energyInFactor, double energyOu
35 static double MCGIDI_KalbachMann_S_a_or_b( dou
36 /*
37 **********************************************
38 */
39 MCGIDI_KalbachMann *MCGIDI_KalbachMann_new( st
40 ptwXY_interpolation interpolationXY )
41
42 MCGIDI_KalbachMann *KalbachMann;
43
44 if( ( KalbachMann = (MCGIDI_KalbachMann *)
45 if( MCGIDI_KalbachMann_initialize( smr, Ka
46 return( KalbachMann );
47 }
48 /*
49 **********************************************
50 */
51 int MCGIDI_KalbachMann_initialize( statusMessa
52
53 memset( KalbachMann, 0, sizeof( MCGIDI_Kal
54 KalbachMann->dists.interpolationWY = inter
55 KalbachMann->dists.interpolationXY = inter
56 return( 0 );
57 }
58 /*
59 **********************************************
60 */
61 MCGIDI_KalbachMann *MCGIDI_KalbachMann_free( s
62
63 MCGIDI_KalbachMann_release( smr, KalbachMa
64 smr_freeMemory( (void **) &KalbachMann );
65 return( NULL );
66 }
67 /*
68 **********************************************
69 */
70 int MCGIDI_KalbachMann_release( statusMessageR
71
72 int i;
73 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMan
74
75 for( i = 0; i < dists->numberOfWs; i++ ) {
76 smr_freeMemory( (void **) &(KalbachMan
77 smr_freeMemory( (void **) &(dists->dis
78 }
79 smr_freeMemory( (void **) &(KalbachMann->r
80 smr_freeMemory( (void **) &(dists->Ws) );
81 smr_freeMemory( (void **) &(dists->dist) )
82
83 MCGIDI_KalbachMann_initialize( smr, Kalbac
84 return( 0 );
85 }
86 /*
87 **********************************************
88 */
89 int MCGIDI_KalbachMann_parseFromTOM( statusMes
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, inter
98 char const *energyFromUnit, *energyToUnit
99
100 MCGIDI_POP *productPOP = distribution->pro
101 double productZ = productPOP->Z, productA
102 MCGIDI_target_heated *targetHeated = MCGID
103 MCGIDI_POP *projectilePOP = MCGIDI_target_
104 double projectileZ = projectilePOP->Z, pro
105 MCGIDI_POP *targetPOP = MCGIDI_target_heat
106 double targetZ = targetPOP->Z, targetA = t
107 double Ia = 0., Ib = 0., Ma = -1, mb = -1;
108
109 if( ( targetA == 0 ) && ( targetZ == 6 ) )
110 targetN = 6;
111 targetA = 12;
112 }
113 if( ( KalbachMannElement = xDataTOME_getOn
114
115 if( MCGIDI_fromTOM_interpolation( smr, Kal
116 if( MCGIDI_fromTOM_interpolation( smr, Kal
117
118 xDataInfo = &(KalbachMannElement->xDataInf
119 KalbachMannXData = (xDataTOM_KalbachMann *
120 if( KalbachMannXData->type == xDataTOM_Kal
121
122 energyFromUnit = xDataTOM_axes_getUnit( sm
123 if( !smr_isOk( smr ) ) goto err;
124 energyInFactor = MCGIDI_misc_getUnitConver
125 if( !smr_isOk( smr ) ) goto err;
126
127 energyFromUnit = xDataTOM_axes_getUnit( sm
128 if( !smr_isOk( smr ) ) goto err;
129 energyOutFactor = MCGIDI_misc_getUnitConve
130 if( !smr_isOk( smr ) ) goto err;
131
132 if( ( KalbachMann = distribution->KalbachM
133
134 /*
135 double productMass MCGIDI_product_getMass_
136 */
137 KalbachMann->energyToMeVFactor = MCGIDI_mi
138 KalbachMann->massFactor = (double) product
139 KalbachMann->massFactor /= projectileN + p
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_o
182 KalbachMann->Sb = MCGIDI_KalbachMann_S_a_o
183 targetZ + projectileZ, targetN + proje
184
185 KalbachMann->dists.numberOfWs = 0;
186 if( ( KalbachMann->dists.Ws = (double *) s
187 if( ( KalbachMann->dists.dist = (MCGIDI_pd
188 if( ( KalbachMann->ras = (MCGIDI_KalbachMa
189
190 for( index = 0; index < KalbachMannXData->
191 if( MCGIDI_KalbachMann_parseFromTOM2(
192 energyInFactor, energyOutFactor, K
193 }
194
195 if( ( KalbachMann->frame = MCGIDI_misc_get
196 distribution->type = MCGIDI_distributionTy
197
198 return( 0 );
199
200 err:
201 if( KalbachMann != NULL ) MCGIDI_KalbachMa
202 return( 1 );
203 }
204 /*
205 **********************************************
206 */
207 static int MCGIDI_KalbachMann_parseFromTOM2( s
208 double energyInFactor, double energyOu
209
210 int i, j, n = coefficientsXData->length /
211 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMan
212 MCGIDI_pdfOfX *dist = &(dists->dist[index]
213 double norm, *p, *rs = NULL, *as_ = NULL,
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
221 pdf = &(Xs[n]);
222 cdf = &(pdf[n]);
223
224 if( ( rs = (double *) smr_malloc2( smr, (
225 if( dataPerEout == 4 ) as_ = &(rs[n]);
226
227 ptwFunc = "ptwXY_new";
228 if( ( pdfXY = ptwXY_new( KalbachMann->dist
229
230 ptwFunc = "ptwXY_setXYPairAtIndex";
231 for( i = 0, p = coefficientsXData->coeffic
232 if( ( status = ptwXY_setValueAtX( pdfX
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
239 Xs[j] = energyOutFactor * point->x;
240 pdf[j] = point->y / energyOutFactor;
241 }
242
243 ptwFunc = "ptwXY_runningIntegral";
244 if( ( cdfX = ptwXY_runningIntegral( pdfXY,
245 norm = ptwX_getPointAtIndex_Unsafely( cdfX
246 if( std::fabs( 1. - norm ) > 0.99 ) {
247 smr_setReportError2( smr, smr_unknownI
248 goto err;
249 }
250 for( j = 0; j < n; j++ ) cdf[j] = ptwX_get
251 for( j = 0; j < n; j++ ) pdf[j] /= norm;
252
253 dists->numberOfWs++;
254 dists->Ws[index] = energyInFactor * coeffi
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
268
269 err:
270 if( Xs != NULL ) smr_freeMemory( (void **)
271 if( rs != NULL ) smr_freeMemory( (void **)
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( dou
280
281 double A_AB = Z_AB + N_AB, A_C = Z_C + N_C
282 double invA_AB_third = 1.0 / G4Pow::GetIns
283 double NZA_AB = ( N_AB - Z_AB ) * ( N_AB -
284
285 S = 15.68 * ( A_C - A_AB )
286 - 18.56 * ( A_C * invA_C_third - A_AB
287 - 0.717 * ( Z_C * Z_C * invA_C_third -
288 - I_ab;
289 return( S );
290 }
291 /*
292 **********************************************
293 */
294 int MCGIDI_KalbachMann_sampleEp( statusMessage
295 MCGIDI_decaySamplingInfo *decaySamplin
296
297 double Epl, Epu, Ep, r, r2, rl, ru, a, a2,
298 MCGIDI_pdfsOfXGivenW_sampled sampled;
299 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMan
300 ptwXY_interpolation interpolationWY;
301
302 sampled.smr = smr;
303 sampled.w = modes.getProjectileEnergy( );
304 MCGIDI_sampling_sampleX_from_pdfsOfXGivenW
305
306 interpolationWY = sampled.interpolationWY;
307 if( sampled.iW < 0 ) {
308 interpolationWY = ptwXY_interpolationF
309 if( sampled.iW == -2 ) { /* ???
310 sampled.iW = 0; }
311 else if( sampled.iW == -1 ) {
312 sampled.iW = dists->numberOfWs - 1
313 }
314 }
315
316 Ep = sampled.x;
317 if( sampled.interpolationXY == ptwXY_inter
318 r = KalbachMann->ras[sampled.iW].rs[sa
319 else {
320 Epl = dists->dist[sampled.iW].Xs[sampl
321 Epu = dists->dist[sampled.iW].Xs[sampl
322 rl = KalbachMann->ras[sampled.iW].rs[s
323 ru = KalbachMann->ras[sampled.iW].rs[s
324 r = ( ru - rl ) / ( Epu - Epl ) * ( Ep
325 }
326 if( interpolationWY == ptwXY_interpolation
327 if( sampled.interpolationXY == ptwXY_i
328 r2 = KalbachMann->ras[sampled.iW+1
329 else {
330 Epl = dists->dist[sampled.iW+1].Xs
331 Epu = dists->dist[sampled.iW+1].Xs
332 rl = KalbachMann->ras[sampled.iW+1
333 ru = KalbachMann->ras[sampled.iW+1
334 r2 = ( ru - rl ) / ( Epu - Epl ) *
335 }
336 r = sampled.frac * r + ( 1. - sampled.
337 }
338
339 if( KalbachMann->ras[0].as == NULL ) {
340 double X1, X3_2;
341 double eb = KalbachMann->massFactor *
342
343 X1 = eb; /* Not valid for
344 X3_2 = eb * eb; /* Not valid for
345 a = X1 * ( C1 + C2 * X1 * X1 ) + C2 *
346 else {
347 if( sampled.interpolationXY == ptwXY_i
348 a = KalbachMann->ras[sampled.iW].a
349 else {
350 Epl = dists->dist[sampled.iW].Xs[s
351 Epu = dists->dist[sampled.iW].Xs[s
352 al = KalbachMann->ras[sampled.iW].
353 au = KalbachMann->ras[sampled.iW].
354 a = ( au - al ) / ( Epu - Epl ) *
355 }
356 a2 = 0.;
357 if( interpolationWY == ptwXY_interpola
358 if( sampled.interpolationXY == ptw
359 a2 = KalbachMann->ras[sampled.
360 else {
361 Epl = dists->dist[sampled.iW+1
362 Epu = dists->dist[sampled.iW+1
363 al = KalbachMann->ras[sampled.
364 au = KalbachMann->ras[sampled.
365 a2 = ( au - al ) / ( Epu - Epl
366 }
367 }
368 a = sampled.frac * a + ( 1. - sampled.
369 }
370
371 /* In the following: Cosh[ a mu ] + r
372 if( decaySamplingInfo->rng( decaySamplingI
373 double T = ( 2. * decaySamplingInfo->r
374
375 mu = G4Log( T + std::sqrt( T * T + 1.
376 else {
377 double rng1 = decaySamplingInfo->rng(
378
379 mu = G4Log( rng1 * exp_a + ( 1. - rng1
380 }
381 if( mu < -1 ) {
382 mu = -1;}
383 else if( mu > 1 ) {
384 mu = 1;
385 }
386
387 decaySamplingInfo->frame = KalbachMann->fr
388 decaySamplingInfo->Ep = Ep;
389 decaySamplingInfo->mu = mu;
390 return( !smr_isOk( smr ) );
391 }
392
393 #if defined __cplusplus
394 }
395 #endif
396