| Geant4 Cross Reference |
1 /* 1
2 # <<BEGIN-copyright>>
3 # <<END-copyright>>
4 */
5 #include <string.h>
6 #define _USE_MATH_DEFINES
7 #include <cmath>
8
9
10 #include "MCGIDI_fromTOM.h"
11 #include "MCGIDI_misc.h"
12 #include "MCGIDI_private.h"
13 #include <nf_specialFunctions.h>
14
15 #if defined __cplusplus
16 #include "G4Exp.hh"
17 #include "G4Log.hh"
18 #include "G4Pow.hh"
19 namespace GIDI {
20 using namespace GIDI;
21 #endif
22
23 static int MCGIDI_energy_parseWeightFromTOM( s
24 static int MCGIDI_energy_parseWeightedFunction
25 static int MCGIDI_energy_parseGeneralEvaporati
26 static int MCGIDI_energy_parseEvaporationFromT
27 static int MCGIDI_energy_parseWattFromTOM( sta
28 static int MCGIDI_energy_parseSimpleMaxwellian
29 static int MCGIDI_energy_parseMadlandNixFromTO
30 static nfu_status MCGIDI_energy_parseMadlandNi
31 static double MCGIDI_energy_parseMadlandNixFro
32 static int MCGIDI_energy_parseNBodyPhaseSpaceF
33 MCGIDI_distribution *distribution );
34
35 static int MCGIDI_energy_sampleSimpleMaxwellia
36 static int MCGIDI_energy_sampleEvaporation( st
37 static int MCGIDI_energy_sampleWatt( statusMes
38 static int MCGIDI_energy_sampleWeightedFunctio
39 MCGIDI_quantitiesLookupModes &modes, MCGID
40 static nfu_status MCGIDI_energy_NBodyPhaseSpac
41 /*
42 **********************************************
43 */
44 MCGIDI_energy *MCGIDI_energy_new( statusMessag
45
46 MCGIDI_energy *energy;
47
48 if( ( energy = (MCGIDI_energy *) smr_mallo
49 if( MCGIDI_energy_initialize( smr, energy
50 return( energy );
51 }
52 /*
53 **********************************************
54 */
55 int MCGIDI_energy_initialize( statusMessageRep
56
57 memset( energy, 0, sizeof( MCGIDI_energy )
58 return( 0 );
59 }
60 /*
61 **********************************************
62 */
63 MCGIDI_energy *MCGIDI_energy_free( statusMessa
64
65 MCGIDI_energy_release( smr, energy );
66 smr_freeMemory( (void **) &energy );
67 return( NULL );
68 }
69 /*
70 **********************************************
71 */
72 int MCGIDI_energy_release( statusMessageReport
73
74 int i;
75
76 MCGIDI_sampling_pdfsOfXGivenW_release( smr
77 if( energy->theta ) energy->theta = ptwXY_
78 if( energy->Watt_a ) energy->Watt_a = ptwX
79 if( energy->Watt_b ) energy->Watt_b = ptwX
80 if( ( energy->type == MCGIDI_energyType_ge
81 MCGIDI_sampling_pdfsOfX_release( smr,
82 else if( energy->type == MCGIDI_energyType
83 for( i = 0; i < energy->weightedFuncti
84 ptwXY_free( energy->weightedFuncti
85 MCGIDI_energy_free( smr, energy->w
86 }
87 }
88
89 MCGIDI_energy_initialize( smr, energy );
90 return( 0 );
91 }
92 /*
93 **********************************************
94 */
95 int MCGIDI_energy_parseFromTOM( statusMessageR
96 enum MCGIDI_energyType energyType, dou
97
98 MCGIDI_energy *energy = NULL;
99 xDataTOM_element *energyElement, *linearEl
100 char const *nativeData;
101 double projectileMass_MeV, targetMass_MeV;
102
103 if( ( energy = MCGIDI_energy_new( smr ) )
104
105 projectileMass_MeV = MCGIDI_product_getPro
106 targetMass_MeV = MCGIDI_product_getTargetM
107 energy->e_inCOMFactor = targetMass_MeV / (
108
109 if( ( energyType == MCGIDI_energyType_prim
110 energy->type = energyType;
111 energy->gammaEnergy_MeV = gammaEnergy_
112 energy->frame = xDataTOM_frame_lab;
113 if( energyType == MCGIDI_energyType_pr
114 else {
115 if( ( energyElement = xDataTOME_getOne
116 if( ( nativeData = xDataTOM_getAttribu
117 if( ( linearElement = xDataTOME_getOne
118 linearElement = xDataTOME_getOneEl
119 if( linearElement == NULL ) {
120 if( ( functional = xDataTOME_getOn
121 if( MCGIDI_energy_parseGeneral
122 else if( ( functional = xDataTOME_
123 if( MCGIDI_energy_parseSimpleM
124 else if( ( functional = xDataTOME_
125 if( MCGIDI_energy_parseEvapora
126 else if( ( functional = xDataTOME_
127 if( MCGIDI_energy_parseWattFro
128 else if( ( functional = xDataTOME_
129 if( MCGIDI_energy_parseMadland
130 else if( ( functional = xDataTOME_
131 if( MCGIDI_energy_parseNBodyPh
132 else if( ( functional = xDataTOME_
133 if( MCGIDI_energy_parseWeighte
134 else {
135 smr_setReportError2( smr, smr_
136 goto err;
137 }
138 frameElement = functional; }
139 else {
140 char const *toUnits[3] = { "MeV",
141
142 frameElement = linearElement;
143 if( MCGIDI_fromTOM_pdfsOfXGivenW(
144 energy->type = MCGIDI_energyType_l
145 }
146 if( ( energy->frame = MCGIDI_misc_getP
147 }
148 distribution->energy = energy;
149
150 return( 0 );
151
152 err:
153 if( energy != NULL ) MCGIDI_energy_free( s
154 return( 1 );
155 }
156 /*
157 **********************************************
158 */
159 static int MCGIDI_energy_parseWeightedFunction
160
161 int i;
162 xDataTOM_element *child;
163
164 for( i = 0, child = xDataTOME_getFirstElem
165 if( strcmp( child->name, "weighted" )
166 if( MCGIDI_energy_parseWeightFromTOM(
167 energy->weightedFunctionals.numberOfWe
168 }
169 energy->type = MCGIDI_energyType_weightedF
170 return( 0 );
171
172 err:
173 return( 1 );
174 }
175 /*
176 **********************************************
177 */
178 static int MCGIDI_energy_parseWeightFromTOM( s
179
180 xDataTOM_element *child;
181 MCGIDI_energy *energy = NULL;
182 ptwXYPoints *weight = NULL;
183 char const *toUnits[2] = { "MeV", "" };
184
185 if( ( energy = MCGIDI_energy_new( smr ) )
186 for( child = xDataTOME_getFirstElement( el
187 if( strcmp( child->name, "weight" ) ==
188 if( ( weight = MCGIDI_misc_dataFro
189 else if( strcmp( child->name, "evapora
190 if( MCGIDI_energy_parseEvaporation
191 else {
192 smr_setReportError2( smr, smr_unkn
193 goto err;
194 }
195 }
196 weightedFunctional->weight = weight;
197 weightedFunctional->energy = energy;
198 return( 0 );
199
200 err:
201 if( weight != NULL ) ptwXY_free( weight );
202 if( energy != NULL ) MCGIDI_energy_free( s
203 return( 1 );
204 }
205 /*
206 **********************************************
207 */
208 static int MCGIDI_energy_parseGeneralEvaporati
209
210 double norm;
211 xDataTOM_element *thetaTOM, *gTOM;
212 ptwXYPoints *theta = NULL, *g = NULL;
213 char const *toUnits[2] = { "MeV", "MeV" };
214
215 if( ( thetaTOM = xDataTOME_getOneElementBy
216 if( ( theta = MCGIDI_misc_dataFromElement2
217
218 if( ( gTOM = xDataTOME_getOneElementByName
219 toUnits[0] = "";
220 toUnits[1] = "";
221 if( ( g = MCGIDI_misc_dataFromElement2ptwX
222 if( MCGIDI_fromTOM_pdfOfX( smr, g, &(energ
223 energy->gInterpolation = ptwXY_getInterpol
224 g = ptwXY_free( g );
225 if( std::fabs( 1. - norm ) > 0.001 ) print
226
227 energy->type = MCGIDI_energyType_generalEv
228 energy->theta = theta;
229 return( 0 );
230
231 err:
232 if( theta != NULL ) ptwXY_free( theta );
233 if( g != NULL ) ptwXY_free( g );
234 return( 1 );
235 }
236 /*
237 **********************************************
238 */
239 static int MCGIDI_energy_parseSimpleMaxwellian
240
241 char const *U, *toUnits[2] = { "MeV", "MeV
242 xDataTOM_element *thetaTOM;
243
244 if( ( U = xDataTOM_getAttributesValueInEle
245 smr_setReportError2( smr, smr_unknownI
246 goto err;
247 }
248 if( MCGIDI_misc_PQUStringToDoubleInUnitOf(
249 if( ( thetaTOM = xDataTOME_getOneElementBy
250 if( ( energy->theta = MCGIDI_misc_dataFrom
251 energy->type = MCGIDI_energyType_simpleMax
252 return( 0 );
253
254 err:
255 return( 1 );
256 }
257 /*
258 **********************************************
259 */
260 static int MCGIDI_energy_parseEvaporationFromT
261
262 char const *U, *toUnits[2] = { "MeV", "MeV
263 xDataTOM_element *thetaTOM;
264
265 if( ( U = xDataTOM_getAttributesValueInEle
266 smr_setReportError2( smr, smr_unknownI
267 goto err;
268 }
269 if( MCGIDI_misc_PQUStringToDoubleInUnitOf(
270 if( ( thetaTOM = xDataTOME_getOneElementBy
271 if( ( energy->theta = MCGIDI_misc_dataFrom
272 energy->type = MCGIDI_energyType_evaporati
273 return( 0 );
274
275 err:
276 return( 1 );
277 }
278 /*
279 **********************************************
280 */
281 static int MCGIDI_energy_parseWattFromTOM( sta
282
283 char const *U, *toUnits[2] = { "MeV", "MeV
284 xDataTOM_element *aOrBTOM;
285
286 if( ( U = xDataTOM_getAttributesValueInEle
287 smr_setReportError2( smr, smr_unknownI
288 goto err;
289 }
290 if( MCGIDI_misc_PQUStringToDoubleInUnitOf(
291
292 if( ( aOrBTOM = xDataTOME_getOneElementByN
293 if( ( energy->Watt_a = MCGIDI_misc_dataFro
294
295 toUnits[1] = "1/MeV";
296 if( ( aOrBTOM = xDataTOME_getOneElementByN
297 if( ( energy->Watt_b = MCGIDI_misc_dataFro
298
299 energy->type = MCGIDI_energyType_Watt;
300 return( 0 );
301
302 err:
303 return( 1 );
304 }
305 /*
306 **********************************************
307 */
308 static int MCGIDI_energy_parseMadlandNixFromTO
309
310 int iE, length, nXs, i1, n;
311 double E=0., T_M=0., EFL=0., EFH=0., argLi
312 xs[] = { 1e-5, 1e-3, 1e-1, 1e1, 1e3
313 ptwXYPoints *ptwXY_TM = NULL, *pdfXY = NUL
314 ptwXYPoint *point;
315 ptwXPoints *cdfX = NULL;
316 nfu_status status = nfu_Okay;
317 xDataTOM_element *TM_TOM;
318 xDataTOM_XYs *XYs;
319 MCGIDI_pdfsOfXGivenW *dists = &(energy->di
320 MCGIDI_pdfOfX *dist;
321 char const *EF, *TMUnits[2] = { "MeV", "Me
322
323 nXs = sizeof( xs ) / sizeof( xs[0] );
324
325 if( ( EF = xDataTOM_getAttributesValueInEl
326 smr_setReportError2( smr, smr_unknownI
327 goto err;
328 }
329 if( MCGIDI_misc_PQUStringToDoubleInUnitOf(
330 argList[0] = EFL;
331
332 if( ( EF = xDataTOM_getAttributesValueInEl
333 smr_setReportError2( smr, smr_unknownI
334 goto err;
335 }
336 if( MCGIDI_misc_PQUStringToDoubleInUnitOf(
337 argList[1] = EFH;
338
339 if( ( TM_TOM = xDataTOME_getOneElementByNa
340 if( ( XYs = (xDataTOM_XYs *) xDataTOME_get
341 if( ( ptwXY_TM = MCGIDI_misc_dataFromXYs2p
342
343 length = (int) ptwXY_length( ptwXY_TM );
344 dists->interpolationWY = ptwXY_interpolati
345 dists->interpolationXY = ptwXY_interpolati
346 if( ( dists->Ws = (double *) smr_malloc2(
347 if( ( dists->dist = (MCGIDI_pdfOfX *) smr_
348
349 for( iE = 0; iE < length; iE++ ) {
350 ptwXY_getXYPairAtIndex( ptwXY_TM, iE,
351 argList[2] = T_M;
352 dist = &(dists->dist[iE]);
353 dists->Ws[iE] = E;
354
355 if( ( pdfXY = ptwXY_createFromFunction
356 (void *) argList, 1e-3, 0, 12, &st
357 if( ( status = ptwXY_normalize( pdfXY
358 smr_setReportError2( smr, smr_unkn
359 goto err;
360 }
361
362 if( ptwXY_simpleCoalescePoints( pdfXY
363 dist->numberOfXs = n = (int) ptwXY_len
364
365 if( ( dist->Xs = (double *) smr_malloc
366 dists->numberOfWs++;
367 dist->pdf = &(dist->Xs[n]);
368 dist->cdf = &(dist->pdf[n]);
369
370 for( i1 = 0; i1 < n; i1++ ) {
371 point = ptwXY_getPointAtIndex_Unsa
372 dist->Xs[i1] = point->x;
373 dist->pdf[i1] = point->y;
374 }
375
376 if( ( cdfX = ptwXY_runningIntegral( pd
377 smr_setReportError2( smr, smr_unkn
378 goto err;
379 }
380
381 norm = ptwX_getPointAtIndex_Unsafely(
382 for( i1 = 0; i1 < n; i1++ ) dist->cdf[
383 for( i1 = 0; i1 < n; i1++ ) dist->pdf[
384 pdfXY = ptwXY_free( pdfXY );
385 cdfX = ptwX_free( cdfX );
386 }
387
388 energy->type = MCGIDI_energyType_MadlandNi
389
390 ptwXY_free( ptwXY_TM );
391 return( 0 );
392
393 err:
394 if( ptwXY_TM != NULL ) ptwXY_free( ptwXY_T
395 if( pdfXY != NULL ) ptwXY_free( pdfXY );
396 if( cdfX != NULL ) cdfX = ptwX_free( cdfX
397
398 return( 1 );
399 }
400 /*
401 **********************************************
402 */
403 static nfu_status MCGIDI_energy_parseMadlandNi
404
405 double *parameters = (double *) argList, E
406 nfu_status status = nfu_Okay;
407
408 EFL = parameters[0];
409 EFH = parameters[1];
410 T_M = parameters[2];
411 *y = MCGIDI_energy_parseMadlandNixFromTOM_
412 if( status == nfu_Okay ) *y += MCGIDI_ener
413 *y *= 0.5;
414 return( status );
415 }
416 /*
417 **********************************************
418 */
419 static double MCGIDI_energy_parseMadlandNixFro
420
421 double u1, u2, E1, E2 = 0., gamma1 = 0., g
422
423 u1 = std::sqrt( Ep ) - std::sqrt( E_F );
424 u1 *= u1 / T_M;
425 u2 = std::sqrt( Ep ) + std::sqrt( E_F );
426 u2 *= u2 / T_M;
427 E1 = 0; /* u1^3/2 * E
428 if( u1 != 0 ) E1 = nf_exponentialIntegral(
429 if( *status == nfu_Okay ) E2 = nf_exponent
430 if( *status != nfu_Okay ) return( 0. );
431 if( u1 > 2. ) {
432 signG = -1;
433 gamma1 = nf_incompleteGammaFunctionCom
434 if( *status == nfu_Okay ) gamma2 = nf_
435 else {
436 gamma1 = nf_incompleteGammaFunction( 1
437 if( *status == nfu_Okay ) gamma2 = nf_
438 }
439 if( *status != nfu_Okay ) return( 0. );
440 return( ( u2 * std::sqrt( u2 ) * E2 - u1 *
441 }
442 /*
443 **********************************************
444 */
445 static int MCGIDI_energy_parseNBodyPhaseSpaceF
446 MCGIDI_distribution *distribution ) {
447
448 int argList[1];
449 double xs[2] = { 0.0, 1.0 }, productMass_M
450 ptwXYPoints *pdf = NULL;
451 nfu_status status;
452 char const *mass;
453
454 if( xDataTOME_convertAttributeToInteger( N
455 if( ( mass = xDataTOM_getAttributesValueIn
456 smr_setReportError2( smr, smr_unknownI
457 goto err;
458 }
459 if( MCGIDI_misc_PQUStringToDouble( smr, ma
460 argList[0] = energy->NBodyPhaseSpace.numbe
461 if( ( pdf = ptwXY_createFromFunction( 2, x
462 smr_setReportError2( smr, smr_unknownI
463 status, nfu_statusMessage( status
464 goto err;
465 }
466 if( MCGIDI_fromTOM_pdfOfX( smr, pdf, &(ene
467 productMass_MeV = MCGIDI_product_getMass_M
468 if( !smr_isOk( smr ) ) goto err;
469 energy->NBodyPhaseSpace.massFactor = ( 1.
470 energy->NBodyPhaseSpace.Q_MeV = MCGIDI_out
471 if( !smr_isOk( smr ) ) goto err;
472
473 ptwXY_free( pdf );
474 energy->type = MCGIDI_energyType_NBodyPhas
475
476 return( 0 );
477
478 err:
479 if( pdf != NULL ) ptwXY_free( pdf );
480 return( 1 );
481 }
482 /*
483 **********************************************
484 */
485 static nfu_status MCGIDI_energy_NBodyPhaseSpac
486
487 int numberOfProducts = *((int *) argList);
488 double e = 0.5 * ( 3 * numberOfProducts -
489
490 *y = std::sqrt( x ) * G4Pow::GetInstance()
491 return( nfu_Okay );
492 }
493 /*
494 **********************************************
495 */
496 int MCGIDI_energy_sampleEnergy( statusMessageR
497 MCGIDI_decaySamplingInfo *decaySamplin
498 /*
499 * This function must be called before angula
500 */
501 double theta, randomEp, Watt_a, Watt_b, e_
502 MCGIDI_pdfsOfXGivenW_sampled sampled;
503
504 decaySamplingInfo->frame = energy->frame;
505 switch( energy->type ) {
506 case MCGIDI_energyType_primaryGamma :
507 decaySamplingInfo->Ep = energy->gammaE
508 break;
509 case MCGIDI_energyType_discreteGamma :
510 decaySamplingInfo->Ep = energy->gammaE
511 break;
512 case MCGIDI_energyType_linear :
513 randomEp = decaySamplingInfo->rng( dec
514 sampled.smr = smr;
515 sampled.w = e_in;
516 MCGIDI_sampling_sampleX_from_pdfsOfXGi
517 decaySamplingInfo->Ep = sampled.x;
518 break;
519 case MCGIDI_energyType_generalEvaporation
520 sampled.interpolationXY = energy->gInt
521 MCGIDI_sampling_sampleX_from_pdfOfX( &
522 theta = MCGIDI_sampling_ptwXY_getValue
523 decaySamplingInfo->Ep = theta * sample
524 break;
525 case MCGIDI_energyType_simpleMaxwellianFis
526 theta = MCGIDI_sampling_ptwXY_getValue
527 MCGIDI_energy_sampleSimpleMaxwellianFi
528 decaySamplingInfo->Ep *= theta;
529 break;
530 case MCGIDI_energyType_evaporation :
531 theta = MCGIDI_sampling_ptwXY_getValue
532 MCGIDI_energy_sampleEvaporation( smr,
533 decaySamplingInfo->Ep *= theta;
534 break;
535 case MCGIDI_energyType_Watt :
536 Watt_a = MCGIDI_sampling_ptwXY_getValu
537 Watt_b = MCGIDI_sampling_ptwXY_getValu
538 MCGIDI_energy_sampleWatt( smr, e_in -
539 break;
540 case MCGIDI_energyType_MadlandNix :
541 MCGIDI_sampling_sampleX_from_pdfsOfXGi
542 decaySamplingInfo->Ep = sampled.x;
543 break;
544 case MCGIDI_energyType_NBodyPhaseSpace :
545 MCGIDI_sampling_sampleX_from_pdfOfX( &
546 decaySamplingInfo->Ep = ( energy->e_in
547 break;
548 case MCGIDI_energyType_weightedFunctional
549 MCGIDI_energy_sampleWeightedFunctional
550 break;
551 default :
552 smr_setReportError2( smr, smr_unknownI
553 }
554
555 return( !smr_isOk( smr ) );
556 }
557 /*
558 **********************************************
559 */
560 static int MCGIDI_energy_sampleSimpleMaxwellia
561
562 int i1;
563 double a = e_in_U_theta, b, c, x, norm_a,
564
565 sqrt_x = std::sqrt( a );
566 norm_a = sqrt_pi_2 * erf( sqrt_x ) - sqrt_
567 b = norm_a * decaySamplingInfo->rng( decay
568 for( i1 = 0; i1 < 16; i1++ ) {
569 x = 0.5 * ( xMin + xMax );
570 sqrt_x = std::sqrt( x );
571 c = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x
572 if( b < c ) {
573 xMax = x; }
574 else {
575 xMin = x;
576 }
577 }
578 /* To order e, the correct x is x + e
579 decaySamplingInfo->Ep = x;
580
581 return( 0 );
582 }
583 /*
584 **********************************************
585 */
586 static int MCGIDI_energy_sampleEvaporation( st
587
588 int i1;
589 double a = e_in_U_theta, b, c, x, norm_a,
590
591 norm_a = 1 - ( 1 + a ) * G4Exp( -a );
592 b = 1. - norm_a * decaySamplingInfo->rng(
593 for( i1 = 0; i1 < 16; i1++ ) {
594 x = 0.5 * ( xMin + xMax );
595 c = ( 1 + x ) * G4Exp( -x );
596 if( b > c ) {
597 xMax = x; }
598 else {
599 xMin = x;
600 }
601 }
602 /* To order e, the correct x is x + e
603 decaySamplingInfo->Ep = x;
604
605 return( 0 );
606 }
607 /*
608 **********************************************
609 */
610 static int MCGIDI_energy_sampleWatt( statusMes
611 /*
612 * From MCAPM via Sample Watt Spectrum as in
613 */
614 double WattMin = 0., WattMax = e_in_U, x,
615
616 x = 1. + ( Watt_b / ( 8. * Watt_a ) );
617 y = ( x + std::sqrt( x * x - 1. ) ) / Watt
618 z = Watt_a * y - 1.;
619 G4int icounter=0;
620 G4int icounter_max=1024;
621 do
622 {
623 icounter++;
624 if ( icounter > icounter_max ) {
625 G4cout << "Loop-counter exceeded the
626 break;
627 }
628 rand1 = -G4Log( decaySamplingInfo->rng
629 rand2 = -G4Log( decaySamplingInfo->rng
630 energyOut = y * rand1;
631 }
632 while( ( ( rand2 - z * ( rand1 + 1. ) ) *
633 decaySamplingInfo->Ep = energyOut;
634
635 return( 0 );
636 }
637 /*
638 **********************************************
639 */
640 static int MCGIDI_energy_sampleWeightedFunctio
641 MCGIDI_quantitiesLookupModes &modes, M
642 /*
643 c This routine assumes that the weights sum
644 */
645 int iW;
646 double rW = decaySamplingInfo->rng( decayS
647 MCGIDI_energyWeightedFunctional *weightedF
648
649 for( iW = 0; iW < energy->weightedFunction
650 weightedFunctional = &(energy->weighte
651 weight = MCGIDI_sampling_ptwXY_getValu
652 cumulativeW += weight;
653 if( cumulativeW >= rW ) break;
654 }
655 return( MCGIDI_energy_sampleEnergy( smr, w
656 }
657
658 #if defined __cplusplus
659 }
660 #endif
661
662