8 #if defined __cplusplus
15 const double C1 = 0.04,
C2 = 1.8e-6;
19 #if defined __cplusplus
23 #include "MCGIDI_fromTOM.h"
24 #include "MCGIDI_misc.h"
25 #include "MCGIDI_private.h"
27 #if defined __cplusplus
34 double energyInFactor,
double energyOutFactor, MCGIDI_KalbachMann *KalbachMann );
40 ptwXY_interpolation interpolationXY ) {
42 MCGIDI_KalbachMann *KalbachMann;
44 if( ( KalbachMann = (MCGIDI_KalbachMann *) smr_malloc2( smr,
sizeof( MCGIDI_KalbachMann ), 0,
"KalbachMann" ) ) == NULL )
return( NULL );
46 return( KalbachMann );
51 int MCGIDI_KalbachMann_initialize( statusMessageReporting * , MCGIDI_KalbachMann *KalbachMann, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY ) {
53 memset( KalbachMann, 0,
sizeof( MCGIDI_KalbachMann ) );
54 KalbachMann->dists.interpolationWY = interpolationWY;
55 KalbachMann->dists.interpolationXY = interpolationXY;
73 MCGIDI_pdfsOfXGivenW *dists = &(KalbachMann->dists);
75 for( i = 0; i < dists->numberOfWs; i++ ) {
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";
100 MCGIDI_POP *productPOP = distribution->product->pop;
101 double productZ = productPOP->Z, productA = productPOP->A, productN = productA - productZ;
104 double projectileZ = projectilePOP->Z, projectileA = projectilePOP->A, projectileN = projectileA - projectileZ;
106 double targetZ = targetPOP->Z, targetA = targetPOP->A, targetN = targetA - targetZ;
107 double Ia = 0., Ib = 0., Ma = -1, mb = -1;
109 if( ( targetA == 0 ) && ( targetZ == 6 ) ) {
118 xDataInfo = &(KalbachMannElement->xDataInfo);
119 KalbachMannXData = (xDataTOM_KalbachMann *) xDataInfo->data;
120 if( KalbachMannXData->type == xDataTOM_KalbachMannType_fra ) dataPerEout = 4;
132 if( ( KalbachMann = distribution->KalbachMann =
MCGIDI_KalbachMann_new( smr, interpolationWY, interpolationXY ) ) == NULL )
goto err;
138 KalbachMann->massFactor = (double) productZ + productN;
139 KalbachMann->massFactor /= projectileN + projectileZ + targetZ + targetN - productZ + productN;
140 KalbachMann->massFactor += 1.;
142 if( projectileZ == 0 ) {
143 if( projectileN == 1 ) Ma = 1; }
144 else if( projectileZ == 1 ) {
145 if( projectileN == 1 ) {
147 else if( projectileN == 2 ) {
150 else if( projectileZ == 2 ) {
151 if( projectileN == 2 ) {
157 if( productZ == 0 ) {
158 if( productN == 1 ) mb = 0.5; }
159 else if( productZ == 1 ) {
160 if( productN == 1 ) {
162 else if( productN == 2 ) {
165 else if( productN == 3 ) {
168 else if( productZ == 2 ) {
169 if( productN == 1 ) {
172 else if( productN == 2 ) {
178 KalbachMann->Ma = Ma;
179 KalbachMann->mb = mb;
183 targetZ + projectileZ, targetN + projectileN, Ib );
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;
190 for( index = 0; index < KalbachMannXData->numberOfEnergies; index++ ) {
192 energyInFactor, energyOutFactor, KalbachMann ) )
goto err;
196 distribution->type = MCGIDI_distributionType_KalbachMann_e;
208 double energyInFactor,
double energyOutFactor, MCGIDI_KalbachMann *KalbachMann ) {
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;
215 ptwXYPoints *pdfXY = NULL;
217 ptwXPoints *cdfX = NULL;
218 char const *ptwFunc =
"";
220 if( ( Xs = (
double *) smr_malloc2( smr, 3 * n *
sizeof(
double ), 0,
"Xs" ) ) == NULL )
goto err;
224 if( ( rs = (
double *) smr_malloc2( smr, ( dataPerEout - 2 ) * n *
sizeof(
double ), 0,
"rs" ) ) == NULL )
goto err;
225 if( dataPerEout == 4 ) as_ = &(rs[
n]);
227 ptwFunc =
"ptwXY_new";
228 if( ( pdfXY =
ptwXY_new( KalbachMann->dists.interpolationXY, NULL, 2., 1e-3, n, 10, &status, 0 ) ) == NULL )
goto errXY;
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;
234 if( dataPerEout == 4 ) as_[i] = p[3];
237 for( j = 0; j <
n; j++ ) {
239 Xs[j] = energyOutFactor * point->x;
240 pdf[j] = point->y / energyOutFactor;
243 ptwFunc =
"ptwXY_runningIntegral";
246 if( std::fabs( 1. - norm ) > 0.99 ) {
247 smr_setReportError2( smr, smr_unknownID, 1,
"bad norm = %e for angular.linear data", norm );
251 for( j = 0; j <
n; j++ ) pdf[j] /= norm;
254 dists->Ws[index] = energyInFactor * coefficientsXData->value;
255 dist->numberOfXs =
n;
259 KalbachMann->ras[index].rs = rs;
260 KalbachMann->ras[index].as = as_;
267 smr_setReportError2( smr, smr_unknownID, 1,
"%s error = %d: %s\n", ptwFunc, status,
nfu_statusMessage( status ) );
273 if( cdfX != NULL ) cdfX =
ptwX_free( cdfX );
281 double A_AB = Z_AB + N_AB, A_C = Z_C + N_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;
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 )
295 MCGIDI_decaySamplingInfo *decaySamplingInfo ) {
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;
303 sampled.w = modes.getProjectileEnergy( );
306 interpolationWY = sampled.interpolationWY;
307 if( sampled.iW < 0 ) {
308 interpolationWY = ptwXY_interpolationFlat;
309 if( sampled.iW == -2 ) {
311 else if( sampled.iW == -1 ) {
312 sampled.iW = dists->numberOfWs - 1;
317 if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
318 r = KalbachMann->ras[sampled.iW].rs[sampled.iX1]; }
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;
326 if( interpolationWY == ptwXY_interpolationLinLin ) {
327 if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
328 r2 = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2]; }
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;
336 r = sampled.frac * r + ( 1. - sampled.frac ) * r2;
339 if( KalbachMann->ras[0].as == NULL ) {
341 double eb = KalbachMann->massFactor * KalbachMann->energyToMeVFactor * Ep + KalbachMann->Sb;
345 a = X1 * (
C1 +
C2 * X1 * X1 ) +
C2 * KalbachMann->Ma * KalbachMann->mb * X3_2 * X3_2; }
347 if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
348 a = KalbachMann->ras[sampled.iW].as[sampled.iX1]; }
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;
357 if( interpolationWY == ptwXY_interpolationLinLin ) {
358 if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
359 a2 = KalbachMann->ras[sampled.iW+1].as[sampled.iX2]; }
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;
368 a = sampled.frac * a + ( 1. - sampled.frac ) * a2;
372 if( decaySamplingInfo->rng( decaySamplingInfo->rngState ) >= r ) {
373 double T = ( 2. * decaySamplingInfo->rng( decaySamplingInfo->rngState ) - 1. ) * std::sinh( a );
375 mu =
G4Log( T + std::sqrt( T * T + 1. ) ) /
a; }
377 double rng1 = decaySamplingInfo->rng( decaySamplingInfo->rngState ), exp_a =
G4Exp( a );
379 mu =
G4Log( rng1 * exp_a + ( 1. - rng1 ) / exp_a ) /
a;
387 decaySamplingInfo->frame = KalbachMann->frame;
388 decaySamplingInfo->Ep = Ep;
389 decaySamplingInfo->mu = mu;
393 #if defined __cplusplus
xDataTOM_element * xDataTOME_getOneElementByName(statusMessageReporting *smr, xDataTOM_element *element, char const *name, int required)
static G4Pow * GetInstance()
MCGIDI_KalbachMann * MCGIDI_KalbachMann_free(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
MCGIDI_POP * MCGIDI_target_heated_getPOPForProjectile(statusMessageReporting *, MCGIDI_target_heated *target)
int MCGIDI_KalbachMann_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution)
int MCGIDI_KalbachMann_sampleEp(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
MCGIDI_POP * MCGIDI_target_heated_getPOPForTarget(statusMessageReporting *, MCGIDI_target_heated *target)
MCGIDI_KalbachMann * MCGIDI_KalbachMann_new(statusMessageReporting *smr, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
int MCGIDI_KalbachMann_release(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
double ptwX_getPointAtIndex_Unsafely(ptwXPoints *ptwX, int64_t index)
int smr_isOk(statusMessageReporting *smr)
int MCGIDI_KalbachMann_initialize(statusMessageReporting *, MCGIDI_KalbachMann *KalbachMann, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
void * smr_freeMemory(void **p)
double MCGIDI_misc_getUnitConversionFactor(statusMessageReporting *smr, char const *fromUnit, char const *toUnit)
char const * xDataTOM_axes_getUnit(statusMessageReporting *smr, xDataTOM_axes *axes, int index)
MCGIDI_target_heated * MCGIDI_product_getTargetHeated(statusMessageReporting *smr, MCGIDI_product *product)
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
G4double G4Log(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
G4double A13(G4double A) const
static double MCGIDI_KalbachMann_S_a_or_b(double Z_AB, double N_AB, double Z_C, double N_C, double I_ab)
const char * nfu_statusMessage(nfu_status status)
int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW(MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double rngValue)
ptwXPoints * ptwXY_runningIntegral(ptwXYPoints *ptwXY, nfu_status *status)
static int MCGIDI_KalbachMann_parseFromTOM2(statusMessageReporting *smr, int dataPerEout, int index, xDataTOM_KalbachMannCoefficients *coefficientsXData, double energyInFactor, double energyOutFactor, MCGIDI_KalbachMann *KalbachMann)
enum xDataTOM_frame MCGIDI_misc_getProductFrame(statusMessageReporting *smr, xDataTOM_element *frameElement)
ptwXYPoints * ptwXY_new(ptwXY_interpolation interpolation, ptwXY_interpolationOtherInfo const *interpolationOtherInfo, double biSectionMax, double accuracy, int64_t primarySize, int64_t secondarySize, nfu_status *status, int userFlag)
int MCGIDI_fromTOM_interpolation(statusMessageReporting *smr, xDataTOM_element *element, int index, ptwXY_interpolation *interpolation)