8 #if defined __cplusplus
15 const double C1 = 0.04,
C2 = 1.8e-6;
19 #if defined __cplusplus
27 #if defined __cplusplus
46 return( KalbachMann );
93 int index, dataPerEout = 3;
94 double energyInFactor, energyOutFactor;
98 char const *energyFromUnit, *energyToUnit =
"MeV";
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);
138 KalbachMann->
massFactor = (double) productZ + productN;
139 KalbachMann->
massFactor /= projectileN + projectileZ + targetZ + targetN - productZ + productN;
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 );
192 energyInFactor, energyOutFactor, KalbachMann ) )
goto err;
210 int i, j,
n = coefficientsXData->
length / dataPerEout;
213 double norm, *
p, *rs = NULL, *as_ = NULL, *Xs = NULL, *pdf, *cdf;
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";
230 ptwFunc =
"ptwXY_setXYPairAtIndex";
231 for( i = 0, p = coefficientsXData->
coefficients; i < n; i++, p += dataPerEout ) {
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 ) {
251 for( j = 0; j <
n; j++ ) pdf[j] /= norm;
254 dists->
Ws[index] = energyInFactor * coefficientsXData->
value;
259 KalbachMann->
ras[index].
rs = rs;
260 KalbachMann->
ras[index].
as = as_;
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 )
297 double Epl, Epu, Ep, r, r2, rl, ru,
a, a2, al, au, mu, randomEp = decaySamplingInfo->
rng( decaySamplingInfo->
rngState );
307 if( sampled.
iW < 0 ) {
309 if( sampled.
iW == -2 ) {
311 else if( sampled.
iW == -1 ) {
318 r = KalbachMann->
ras[sampled.
iW].
rs[sampled.
iX1]; }
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;
328 r2 = KalbachMann->
ras[sampled.
iW+1].
rs[sampled.
iX2]; }
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 ) {
345 a = X1 * (
C1 +
C2 * X1 *
X1 ) +
C2 * KalbachMann->
Ma * KalbachMann->
mb * X3_2 * X3_2; }
348 a = KalbachMann->
ras[sampled.
iW].
as[sampled.
iX1]; }
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;
359 a2 = KalbachMann->
ras[sampled.
iW+1].
as[sampled.
iX2]; }
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
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
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)
ptwXPoints * ptwX_free(ptwXPoints *ptwX)
ptwXY_interpolation interpolationXY
ptwXY_interpolation interpolationWY
std::vector< ExP01TrackerHit * > a
MCGIDI_target_heated * MCGIDI_product_getTargetHeated(statusMessageReporting *smr, MCGIDI_product *product)
G4double A13(G4double A) const
int MCGIDI_KalbachMann_release(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
double MCGIDI_misc_getUnitConversionFactor(statusMessageReporting *smr, char const *fromUnit, char const *toUnit)
#define smr_setReportError2(smr, libraryID, code, fmt,...)
enum ptwXY_interpolation_e ptwXY_interpolation
int MCGIDI_fromTOM_interpolation(statusMessageReporting *smr, xDataTOM_element *element, int index, enum ptwXY_interpolation_e *interpolation)
MCGIDI_POP * MCGIDI_target_heated_getPOPForProjectile(statusMessageReporting *smr, MCGIDI_target_heated *target)
statusMessageReporting * smr
MCGIDI_pdfsOfXGivenW dists
ptwXPoints * ptwXY_runningIntegral(ptwXYPoints *ptwXY, nfu_status *status)
double ptwX_getPointAtIndex_Unsafely(ptwXPoints *ptwX, int64_t index)
static double MCGIDI_KalbachMann_S_a_or_b(double Z_AB, double N_AB, double Z_C, double N_C, double I_ab)
enum MCGIDI_distributionType type
G4double G4Log(G4double x)
int MCGIDI_KalbachMann_initialize(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
MCGIDI_KalbachMann * MCGIDI_KalbachMann_free(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
xDataTOM_xDataInfo xDataInfo
xDataTOM_KalbachMannCoefficients * coefficients
ptwXY_interpolation interpolationWY
#define smr_malloc2(smr, size, zero, forItem)
xDataTOM_element * xDataTOME_getOneElementByName(statusMessageReporting *smr, xDataTOM_element *element, char const *name, int required)
enum xDataTOM_frame MCGIDI_misc_getProductFrame(statusMessageReporting *smr, xDataTOM_element *frameElement)
static G4Pow * GetInstance()
MCGIDI_POP * MCGIDI_target_heated_getPOPForTarget(statusMessageReporting *smr, MCGIDI_target_heated *target)
int MCGIDI_sampling_sampleX_from_pdfsOfXGivenW(MCGIDI_pdfsOfXGivenW *dists, MCGIDI_pdfsOfXGivenW_sampled *sampled, double r)
int smr_isOk(statusMessageReporting *smr)
enum nfu_status_e nfu_status
MCGIDI_KalbachMann * KalbachMann
enum xDataTOM_frame frame
int MCGIDI_KalbachMann_sampleEp(statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann, MCGIDI_quantitiesLookupModes &modes, MCGIDI_decaySamplingInfo *decaySamplingInfo)
void * smr_freeMemory(void **p)
ptwXY_interpolation interpolationXY
ptwXYPoint * ptwXY_getPointAtIndex_Unsafely(ptwXYPoints *ptwXY, int64_t index)
ptwXYPoints * ptwXY_free(ptwXYPoints *ptwXY)
static int MCGIDI_KalbachMann_parseFromTOM2(statusMessageReporting *smr, int dataPerEout, int index, xDataTOM_KalbachMannCoefficients *coefficientsXData, double energyInFactor, double energyOutFactor, MCGIDI_KalbachMann *KalbachMann)
double getProjectileEnergy(void) const
int MCGIDI_KalbachMann_parseFromTOM(statusMessageReporting *smr, xDataTOM_element *element, MCGIDI_distribution *distribution)
const char * nfu_statusMessage(nfu_status status)
enum xDataTOM_KalbachMannType type
nfu_status ptwXY_setValueAtX(ptwXYPoints *ptwXY, double x, double y)
enum xDataTOM_frame frame
MCGIDI_KalbachMann_ras * ras
MCGIDI_KalbachMann * MCGIDI_KalbachMann_new(statusMessageReporting *smr, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
char const * xDataTOM_axes_getUnit(statusMessageReporting *smr, xDataTOM_axes *axes, int index)