MCGIDI_KalbachMann.cc File Reference

#include <string.h>
#include <cmath>
#include "MCGIDI_fromTOM.h"
#include "MCGIDI_misc.h"
#include "MCGIDI_private.h"

Go to the source code of this file.

Functions
MCGIDI_KalbachMann *	MCGIDI_KalbachMann_new (statusMessageReporting *smr, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
int	MCGIDI_KalbachMann_initialize (statusMessageReporting *, MCGIDI_KalbachMann *KalbachMann, ptwXY_interpolation interpolationWY, ptwXY_interpolation interpolationXY)
MCGIDI_KalbachMann *	MCGIDI_KalbachMann_free (statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
int	MCGIDI_KalbachMann_release (statusMessageReporting *smr, MCGIDI_KalbachMann *KalbachMann)
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)

Variables
const double	C1 = 0.04
const double	C2 = 1.8e-6

Function Documentation

MCGIDI_KalbachMann_free()

MCGIDI_KalbachMann * MCGIDI_KalbachMann_free	(	statusMessageReporting *	smr,
		MCGIDI_KalbachMann *	KalbachMann
	)

Definition at line 61 of file MCGIDI_KalbachMann.cc.

                                                                                                            {
 
    MCGIDI_KalbachMann_release( smr, KalbachMann );
    smr_freeMemory( (void **) &KalbachMann );
    return( NULL );
}

Referenced by MCGIDI_distribution_release(), MCGIDI_KalbachMann_new(), and MCGIDI_KalbachMann_parseFromTOM().

MCGIDI_KalbachMann_initialize()

int MCGIDI_KalbachMann_initialize	(	statusMessageReporting *	smr,
		MCGIDI_KalbachMann *	KalbachMann,
		ptwXY_interpolation	interpolationWY,
		ptwXY_interpolation	interpolationXY
	)

Definition at line 51 of file MCGIDI_KalbachMann.cc.

                                                                                                                                                                           {
 
    memset( KalbachMann, 0, sizeof( MCGIDI_KalbachMann ) );
    KalbachMann->dists.interpolationWY = interpolationWY;
    KalbachMann->dists.interpolationXY = interpolationXY;
    return( 0 );
}

Referenced by MCGIDI_KalbachMann_new(), and MCGIDI_KalbachMann_release().

MCGIDI_KalbachMann_new()

MCGIDI_KalbachMann * MCGIDI_KalbachMann_new	(	statusMessageReporting *	smr,
		ptwXY_interpolation	interpolationWY,
		ptwXY_interpolation	interpolationXY
	)

Definition at line 39 of file MCGIDI_KalbachMann.cc.

                                              {
 
    MCGIDI_KalbachMann *KalbachMann;
 
    if( ( KalbachMann = (MCGIDI_KalbachMann *) smr_malloc2( smr, sizeof( MCGIDI_KalbachMann ), 0, "KalbachMann" ) ) == NULL ) return( NULL );
    if( MCGIDI_KalbachMann_initialize( smr, KalbachMann, interpolationWY, interpolationXY ) ) KalbachMann = MCGIDI_KalbachMann_free( smr, KalbachMann );
    return( KalbachMann );
}

Referenced by MCGIDI_KalbachMann_parseFromTOM().

MCGIDI_KalbachMann_parseFromTOM()

int MCGIDI_KalbachMann_parseFromTOM	(	statusMessageReporting *	smr,
		xDataTOM_element *	element,
		MCGIDI_distribution *	distribution
	)

Definition at line 89 of file MCGIDI_KalbachMann.cc.

                                                                                                                                 {
 
    MCGIDI_KalbachMann *KalbachMann = NULL;
    xDataTOM_element *KalbachMannElement;
    int index, dataPerEout = 3;
    double energyInFactor, energyOutFactor;
    xDataTOM_xDataInfo *xDataInfo;
    xDataTOM_KalbachMann *KalbachMannXData;
    ptwXY_interpolation interpolationXY, interpolationWY;
    char const *energyFromUnit, *energyToUnit = "MeV";
 
    MCGIDI_POP *productPOP = distribution->product->pop;
    double productZ = productPOP->Z, productA = productPOP->A, productN = productA - productZ;
    MCGIDI_target_heated *targetHeated = MCGIDI_product_getTargetHeated( smr, distribution->product );
    MCGIDI_POP *projectilePOP = MCGIDI_target_heated_getPOPForProjectile( smr, targetHeated );
    double projectileZ = projectilePOP->Z, projectileA = projectilePOP->A, projectileN = projectileA - projectileZ;
    MCGIDI_POP *targetPOP = MCGIDI_target_heated_getPOPForTarget( smr, targetHeated );
    double targetZ = targetPOP->Z, targetA = targetPOP->A, targetN = targetA - targetZ;
    double Ia = 0., Ib = 0., Ma = -1, mb = -1;
 
    if( ( targetA == 0 ) && ( targetZ == 6 ) ) {    /* Special case for C_000 evaluation. */
        targetN = 6;
        targetA = 12;
    }
    if( ( KalbachMannElement = xDataTOME_getOneElementByName( smr, element, "KalbachMann", 1 ) ) == NULL ) goto err;
 
    if( MCGIDI_fromTOM_interpolation( smr, KalbachMannElement, 0, &interpolationWY ) ) goto err;
    if( MCGIDI_fromTOM_interpolation( smr, KalbachMannElement, 1, &interpolationXY ) ) goto err;
 
    xDataInfo = &(KalbachMannElement->xDataInfo);
    KalbachMannXData = (xDataTOM_KalbachMann *) xDataInfo->data;
    if( KalbachMannXData->type == xDataTOM_KalbachMannType_fra ) dataPerEout = 4;
 
    energyFromUnit = xDataTOM_axes_getUnit( smr, &(xDataInfo->axes), 0 );
    if( !smr_isOk( smr ) ) goto err;
    energyInFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyFromUnit, energyToUnit );
    if( !smr_isOk( smr ) ) goto err;
 
    energyFromUnit = xDataTOM_axes_getUnit( smr, &(xDataInfo->axes), 1 );
    if( !smr_isOk( smr ) ) goto err;
    energyOutFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyFromUnit, energyToUnit );
    if( !smr_isOk( smr ) ) goto err;
 
    if( ( KalbachMann = distribution->KalbachMann = MCGIDI_KalbachMann_new( smr, interpolationWY, interpolationXY ) ) == NULL ) goto err;
 
/*
    double productMass MCGIDI_product_getMass_MeV( smr, distribution->product ), residualMass;
*/
    KalbachMann->energyToMeVFactor = MCGIDI_misc_getUnitConversionFactor( smr, energyToUnit, "MeV" );
    KalbachMann->massFactor = (double) productZ + productN;     /* This is not correct as masses are needed not Z and N. */
    KalbachMann->massFactor /= projectileN + projectileZ + targetZ + targetN - productZ + productN;
    KalbachMann->massFactor += 1.;
 
    if( projectileZ == 0 ) {
        if( projectileN == 1 ) Ma = 1; }
    else if( projectileZ == 1 ) {
        if( projectileN == 1 ) {
            Ma = 1; }
        else if( projectileN == 2 ) {
            Ia = 2.22;
            Ma = 1; } }
    else if( projectileZ == 2 ) {
        if( projectileN == 2 ) {
            Ia = 28.3;
            Ma = 0;
        }
    }
        
    if( productZ == 0 ) {
        if( productN == 1 ) mb = 0.5; }
    else if( productZ == 1 ) {
        if( productN == 1 ) {
            mb = 1; }
        else if( productN == 2 ) {
            Ia = 2.22;
            mb = 1; }
        else if( productN == 3 ) {
            Ib = 8.48;
            mb = 1; } }
    else if( productZ == 2 ) {
        if( productN == 1 ) {
            Ib = 7.72;
            mb = 1; }
        else if( productN == 2 ) {
            Ib = 28.3;
            mb = 2;
        }
    }
        
    KalbachMann->Ma = Ma;
    KalbachMann->mb = mb;
 
    KalbachMann->Sa = MCGIDI_KalbachMann_S_a_or_b( targetZ, targetN, targetZ + projectileZ, targetN + projectileN, Ia );
    KalbachMann->Sb = MCGIDI_KalbachMann_S_a_or_b( projectileZ + targetZ - productZ, projectileN + targetN - productN, 
        targetZ + projectileZ, targetN + projectileN, Ib );
 
    KalbachMann->dists.numberOfWs = 0;
    if( ( KalbachMann->dists.Ws = (double *) smr_malloc2( smr, KalbachMannXData->numberOfEnergies * sizeof( double ), 0, "KalbachMann->dists->Ws" ) ) == NULL ) goto err;
    if( ( KalbachMann->dists.dist = (MCGIDI_pdfOfX *) smr_malloc2( smr, KalbachMannXData->numberOfEnergies * sizeof( MCGIDI_pdfOfX ), 0, "KalbachMann->dists->dist" ) ) == NULL ) goto err;
    if( ( KalbachMann->ras = (MCGIDI_KalbachMann_ras *) smr_malloc2( smr, KalbachMannXData->numberOfEnergies * sizeof( MCGIDI_KalbachMann_ras ), 0, "KalbachMann->ras" ) ) == NULL ) goto err;
 
    for( index = 0; index < KalbachMannXData->numberOfEnergies; index++ ) {
        if( MCGIDI_KalbachMann_parseFromTOM2( smr, dataPerEout, index, &(KalbachMannXData->coefficients[index]), 
            energyInFactor, energyOutFactor, KalbachMann ) ) goto err;
    }
 
    if( ( KalbachMann->frame = MCGIDI_misc_getProductFrame( smr, KalbachMannElement ) ) == xDataTOM_frame_invalid ) goto err;
    distribution->type = MCGIDI_distributionType_KalbachMann_e;
 
    return( 0 );
 
err:
    if( KalbachMann != NULL ) MCGIDI_KalbachMann_free( smr, KalbachMann );
    return( 1 );
}

Referenced by MCGIDI_energyAngular_parseFromTOM().

MCGIDI_KalbachMann_release()

int MCGIDI_KalbachMann_release	(	statusMessageReporting *	smr,
		MCGIDI_KalbachMann *	KalbachMann
	)

Definition at line 70 of file MCGIDI_KalbachMann.cc.

                                                                                               {
 
    int i;
    MCGIDI_pdfsOfXGivenW *dists = &(KalbachMann->dists);
 
    for( i = 0; i < dists->numberOfWs; i++ ) {
        smr_freeMemory( (void **) &(KalbachMann->ras[i].rs) );
        smr_freeMemory( (void **) &(dists->dist[i].Xs) );
    }
    smr_freeMemory( (void **) &(KalbachMann->ras) );
    smr_freeMemory( (void **) &(dists->Ws) );
    smr_freeMemory( (void **) &(dists->dist) );
 
    MCGIDI_KalbachMann_initialize( smr, KalbachMann, ptwXY_interpolationLinLin, ptwXY_interpolationLinLin );
    return( 0 );
}

Referenced by MCGIDI_KalbachMann_free().

MCGIDI_KalbachMann_sampleEp()

int MCGIDI_KalbachMann_sampleEp	(	statusMessageReporting *	smr,
		MCGIDI_KalbachMann *	KalbachMann,
		MCGIDI_quantitiesLookupModes &	modes,
		MCGIDI_decaySamplingInfo *	decaySamplingInfo
	)

Definition at line 294 of file MCGIDI_KalbachMann.cc.

                                                      {
 
    double Epl, Epu, Ep, r, r2, rl, ru, a, a2, al, au, mu, randomEp = decaySamplingInfo->rng( decaySamplingInfo->rngState );
    MCGIDI_pdfsOfXGivenW_sampled sampled;
    MCGIDI_pdfsOfXGivenW *dists = &(KalbachMann->dists);
    ptwXY_interpolation interpolationWY;
 
    sampled.smr = smr;
    sampled.w = modes.getProjectileEnergy( );
    MCGIDI_sampling_sampleX_from_pdfsOfXGivenW( dists, &sampled, randomEp );
 
    interpolationWY = sampled.interpolationWY;
    if( sampled.iW < 0 ) {
        interpolationWY = ptwXY_interpolationFlat;
        if( sampled.iW == -2 ) {        /* ???????????? This should probably report a warning. */
            sampled.iW = 0; }
        else if( sampled.iW == -1 ) {
            sampled.iW = dists->numberOfWs - 1;
        }
    }
 
    Ep = sampled.x;                                                 /* Sampled Ep. */
    if( sampled.interpolationXY == ptwXY_interpolationFlat ) {      /* Now sample r. */
        r = KalbachMann->ras[sampled.iW].rs[sampled.iX1]; }
    else {
        Epl = dists->dist[sampled.iW].Xs[sampled.iX1];
        Epu = dists->dist[sampled.iW].Xs[sampled.iX1+1];
        rl = KalbachMann->ras[sampled.iW].rs[sampled.iX1];
        ru = KalbachMann->ras[sampled.iW].rs[sampled.iX1+1];
        r = ( ru - rl ) / ( Epu - Epl ) * ( Ep - Epl ) + rl;
    }
    if( interpolationWY == ptwXY_interpolationLinLin ) {
        if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
            r2 = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2]; }
        else {
            Epl = dists->dist[sampled.iW+1].Xs[sampled.iX2];
            Epu = dists->dist[sampled.iW+1].Xs[sampled.iX2+1];
            rl = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2];
            ru = KalbachMann->ras[sampled.iW+1].rs[sampled.iX2+1];
            r2 = ( ru - rl ) / ( Epu - Epl ) * ( Ep - Epl ) + rl;
        }
        r = sampled.frac * r + ( 1. - sampled.frac ) * r2;
    }
 
    if( KalbachMann->ras[0].as == NULL ) {                          /* Now determine a. */
        double X1, X3_2;
        double eb = KalbachMann->massFactor * KalbachMann->energyToMeVFactor * Ep + KalbachMann->Sb;
 
        X1 = eb;             /* Not valid for ea > Et1. */
        X3_2 = eb * eb;       /* Not valid for ea > Et3. */
        a = X1 * ( C1 + C2 * X1 * X1 ) + C2 * KalbachMann->Ma * KalbachMann->mb * X3_2 * X3_2; }
    else {
        if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
            a = KalbachMann->ras[sampled.iW].as[sampled.iX1]; }
        else {
            Epl = dists->dist[sampled.iW].Xs[sampled.iX1];
            Epu = dists->dist[sampled.iW].Xs[sampled.iX1+1];
            al = KalbachMann->ras[sampled.iW].as[sampled.iX1];
            au = KalbachMann->ras[sampled.iW].as[sampled.iX1+1];
            a = ( au - al ) / ( Epu - Epl ) * ( Ep - Epl ) + al;
        }
        a2 = 0.;
        if( interpolationWY == ptwXY_interpolationLinLin ) {
            if( sampled.interpolationXY == ptwXY_interpolationFlat ) {
                a2 = KalbachMann->ras[sampled.iW+1].as[sampled.iX2]; }
            else {
                Epl = dists->dist[sampled.iW+1].Xs[sampled.iX2];
                Epu = dists->dist[sampled.iW+1].Xs[sampled.iX2+1];
                al = KalbachMann->ras[sampled.iW+1].as[sampled.iX2];
                au = KalbachMann->ras[sampled.iW+1].as[sampled.iX2+1];
                a2 = ( au - al ) / ( Epu - Epl ) * ( Ep - Epl ) + al;
            }
        }
        a = sampled.frac * a + ( 1. - sampled.frac ) * a2;
    }
 
        /* In the following: Cosh[ a mu ] + r Sinh[ a mu ] = ( 1 - r ) Cosh[ a mu ] + r ( Cosh[ a mu ] + Sinh[ a mu ] ). */
    if( decaySamplingInfo->rng( decaySamplingInfo->rngState ) >= r ) {   /* Sample the '( 1 - r ) Cosh[ a mu ]' term. */
        double T = ( 2. * decaySamplingInfo->rng( decaySamplingInfo->rngState ) - 1. ) * std::sinh( a );
 
        mu = G4Log( T + std::sqrt( T * T + 1. ) ) / a; }
    else {                                                              /* Sample the 'r ( Cosh[ a mu ] + Sinh[ a mu ]' term. */
        double rng1 = decaySamplingInfo->rng( decaySamplingInfo->rngState ), exp_a = G4Exp( a );
 
        mu = G4Log( rng1 * exp_a + ( 1. - rng1 ) / exp_a ) / a;
    }
    if( mu < -1 ) {
        mu = -1;}
    else if( mu >  1 ) {
        mu = 1;
    }
 
    decaySamplingInfo->frame = KalbachMann->frame;
    decaySamplingInfo->Ep = Ep;
    decaySamplingInfo->mu = mu;
    return( !smr_isOk( smr ) );
}

Referenced by MCGIDI_outputChannel_sampleProductsAtE().

Variable Documentation

C1

const double C1 = 0.04

Definition at line 15 of file MCGIDI_KalbachMann.cc.

C2

const double C2 = 1.8e-6

Definition at line 15 of file MCGIDI_KalbachMann.cc.

Referenced by G4ParticleHPKallbachMannSyst::A(), G4ecpssrBaseKxsModel::CalculateCrossSection(), G4StrawTubeXTRadiator::GetStackFactor(), MCGIDI_KalbachMann_sampleEp(), G4GammaConversionToMuons::PostStepDoIt(), and G4Abla::width().