CexmcReimplementedGenbod Class Reference

#include <CexmcReimplementedGenbod.hh>

Inheritance diagram for CexmcReimplementedGenbod:

Collaboration diagram for CexmcReimplementedGenbod:

Public Member Functions
	CexmcReimplementedGenbod ()
G4double	Generate (void)
Public Member Functions inherited from CexmcPhaseSpaceGenerator
	CexmcPhaseSpaceGenerator ()
virtual	~CexmcPhaseSpaceGenerator ()
virtual G4bool	CheckKinematics (void)
void	SetParticles (const CexmcPhaseSpaceInVector &inVec_, const CexmcPhaseSpaceOutVector &outVec_)
void	SetFermiEnergyDependence (G4bool on=true)

Additional Inherited Members
Protected Attributes inherited from CexmcPhaseSpaceGenerator
CexmcPhaseSpaceInVector	inVec
CexmcPhaseSpaceOutVector	outVec
G4bool	fermiEnergyDepIsOn
G4double	totalEnergy
G4double	totalMass

Detailed Description

Definition at line 51 of file CexmcReimplementedGenbod.hh.

Constructor & Destructor Documentation

CexmcReimplementedGenbod::CexmcReimplementedGenbod

(

)

Definition at line 81 of file CexmcReimplementedGenbod.cc.

                                                   : maxWeight( 0. ),
    nmbOfOutputParticles( 0 )
{
}

Member Function Documentation

G4double CexmcReimplementedGenbod::Generate ( void  )

virtual

Implements CexmcPhaseSpaceGenerator.

Definition at line 87 of file CexmcReimplementedGenbod.cc.

{
    // Generate a random final state.
    // The function returns the weigth of the current event.
    // Note that Momentum, Energy units are Gev/C, GeV

    G4double  te_minus_tm( totalEnergy - totalMass );
    G4double  rno[ maxParticles ];
    rno[ 0 ] = 0;

    if ( nmbOfOutputParticles > 2 )
    {
        for ( G4int  i( 1 ); i < nmbOfOutputParticles - 1; ++i )
        {
            rno[ i ] = G4UniformRand();
        }
        qsort( rno + 1, nmbOfOutputParticles - 2, sizeof( G4double ),
               DoubleMax );
    }
    rno[ nmbOfOutputParticles - 1 ] = 1;

    G4double  invMas[ maxParticles ];
    G4double  sum( 0 );

    for ( int  i( 0 ); i < nmbOfOutputParticles; ++i )
    {
        sum += outVec[ i ].mass / GeV;
        invMas[ i ] = rno[ i ] * te_minus_tm / GeV + sum;
    }

    //
    //-----> compute the weight of the current event
    //
    G4double  wt( maxWeight );
    G4double  pd[ maxParticles ];

    for ( int  i( 0 ); i < nmbOfOutputParticles - 1; ++i )
    {
        pd[ i ] = PDK( invMas[ i + 1 ], invMas[ i ],
                       outVec[ i + 1 ].mass / GeV );
        wt *= pd[ i ];
    }

    //
    //-----> complete specification of event (Raubold-Lynch method)
    //
    outVec[ 0 ].lVec->setPx( 0. );
    outVec[ 0 ].lVec->setPy( pd[ 0 ] );
    outVec[ 0 ].lVec->setPz( 0. );
    outVec[ 0 ].lVec->setE( std::sqrt( pd[ 0 ] * pd[ 0 ] +
                                       outVec[ 0 ].mass / GeV *
                                       outVec[ 0 ].mass / GeV ) );

    G4int  i( 1 );

    while ( true )
    {
        outVec[ i ].lVec->setPx( 0. );
        outVec[ i ].lVec->setPy( -pd[ i - 1 ] );
        outVec[ i ].lVec->setPz( 0. );
        outVec[ i ].lVec->setE( std::sqrt( pd[ i - 1 ] * pd[ i - 1 ] +
                                           outVec[ i ].mass / GeV *
                                           outVec[ i ].mass / GeV ) );

        G4double  cZ( 2 * G4UniformRand() - 1 );
        G4double  sZ( std::sqrt( 1 - cZ * cZ ) );
        G4double  angY( 2 * pi * G4UniformRand() );
        G4double  cY( std::cos( angY ) );
        G4double  sY( std::sin( angY ) );

        for ( int  j( 0 ); j <= i; ++j )
        {
            G4LorentzVector *  v( outVec[ j ].lVec );
            G4double           x( v->px() );
            G4double           y( v->py() );
            v->setPx( cZ * x - sZ * y );
            v->setPy( sZ * x + cZ * y );   // rotation around Z
            x = v->px();
            G4double           z( v->pz() );
            v->setPx( cY * x - sY * z );
            v->setPz( sY * x + cY * z );   // rotation around Y
        }

        if ( i == nmbOfOutputParticles - 1 )
            break;

        G4double  beta( pd[ i ] / std::sqrt( pd[ i ] * pd[ i ] +
                                             invMas[ i ] * invMas[ i ] ) );
        for ( int  j( 0 ); j <= i; ++j )
            outVec[ j ].lVec->boost( 0, beta, 0 );

        ++i;
    }

    for ( int  i( 0 ); i < nmbOfOutputParticles; ++i )
        *outVec[ i ].lVec *= GeV;

    //
    //---> return the weigth of event
    //
    return wt;
}

Here is the call graph for this function:

---

The documentation for this class was generated from the following files:

• source/geant4.10.03.p02/examples/advanced/ChargeExchangeMC/include/CexmcReimplementedGenbod.hh
• source/geant4.10.03.p02/examples/advanced/ChargeExchangeMC/src/CexmcReimplementedGenbod.cc