G4GeneratorPrecompoundInterface Class Reference

#include <G4GeneratorPrecompoundInterface.hh>

Inheritance diagram for G4GeneratorPrecompoundInterface:

Public Member Functions
	G4GeneratorPrecompoundInterface (G4VPreCompoundModel *p=0)
virtual	~G4GeneratorPrecompoundInterface ()
virtual G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
virtual G4ReactionProductVector *	Propagate (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
virtual G4ReactionProductVector *	PropagateNuclNucl (G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus, G4V3DNucleus *theProjectileNucleus)
void	SetCaptureThreshold (G4double)
virtual void	PropagateModelDescription (std::ostream &) const
Public Member Functions inherited from G4VIntraNuclearTransportModel
	G4VIntraNuclearTransportModel (const G4String &modelName="CascadeModel", G4VPreCompoundModel *ptr=0)
virtual	~G4VIntraNuclearTransportModel ()
void	SetDeExcitation (G4VPreCompoundModel *ptr)
void	Set3DNucleus (G4V3DNucleus *const value)
void	SetPrimaryProjectile (const G4HadProjectile &aPrimary)
const G4String &	GetModelName () const
virtual void	ModelDescription (std::ostream &outFile) const
Public Member Functions inherited from G4HadronicInteraction
	G4HadronicInteraction (const G4String &modelName="HadronicModel")
virtual	~G4HadronicInteraction ()
virtual G4double	SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
virtual G4bool	IsApplicable (const G4HadProjectile &, G4Nucleus &)
G4double	GetMinEnergy () const
G4double	GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
void	SetMinEnergy (G4double anEnergy)
void	SetMinEnergy (G4double anEnergy, const G4Element *anElement)
void	SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
G4double	GetMaxEnergy () const
G4double	GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
void	SetMaxEnergy (const G4double anEnergy)
void	SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
void	SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
const G4HadronicInteraction *	GetMyPointer () const
virtual G4int	GetVerboseLevel () const
virtual void	SetVerboseLevel (G4int value)
const G4String &	GetModelName () const
void	DeActivateFor (const G4Material *aMaterial)
void	ActivateFor (const G4Material *aMaterial)
void	DeActivateFor (const G4Element *anElement)
void	ActivateFor (const G4Element *anElement)
G4bool	IsBlocked (const G4Material *aMaterial) const
G4bool	IsBlocked (const G4Element *anElement) const
void	SetRecoilEnergyThreshold (G4double val)
G4double	GetRecoilEnergyThreshold () const
G4bool	operator== (const G4HadronicInteraction &right) const
G4bool	operator!= (const G4HadronicInteraction &right) const
virtual const std::pair < G4double, G4double >	GetFatalEnergyCheckLevels () const
virtual std::pair< G4double, G4double >	GetEnergyMomentumCheckLevels () const
void	SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)

Additional Inherited Members
Protected Member Functions inherited from G4VIntraNuclearTransportModel
G4V3DNucleus *	Get3DNucleus () const
G4VPreCompoundModel *	GetDeExcitation () const
const G4HadProjectile *	GetPrimaryProjectile () const
Protected Member Functions inherited from G4HadronicInteraction
void	SetModelName (const G4String &nam)
G4bool	IsBlocked () const
void	Block ()
Protected Attributes inherited from G4VIntraNuclearTransportModel
G4String	theTransportModelName
G4V3DNucleus *	the3DNucleus
G4VPreCompoundModel *	theDeExcitation
const G4HadProjectile *	thePrimaryProjectile
Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState	theParticleChange
G4int	verboseLevel
G4double	theMinEnergy
G4double	theMaxEnergy
G4bool	isBlocked

Detailed Description

Definition at line 60 of file G4GeneratorPrecompoundInterface.hh.

Constructor & Destructor Documentation

G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface

(

G4VPreCompoundModel *

p = 0

)

Definition at line 73 of file G4GeneratorPrecompoundInterface.cc.

References G4Alpha::Alpha(), G4AntiAlpha::AntiAlpha(), G4AntiDeuteron::AntiDeuteron(), G4AntiHe3::AntiHe3(), G4AntiNeutron::AntiNeutron(), G4AntiProton::AntiProton(), G4AntiTriton::AntiTriton(), G4Deuteron::Deuteron(), G4HadronicInteractionRegistry::FindModel(), G4He3::He3(), G4HadronicInteractionRegistry::Instance(), G4Neutron::Neutron(), G4Proton::Proton(), G4VIntraNuclearTransportModel::SetDeExcitation(), and G4Triton::Triton().

: CaptureThreshold(10*MeV)
{
    proton = G4Proton::Proton();
    neutron = G4Neutron::Neutron();

    deuteron=G4Deuteron::Deuteron();
    triton  =G4Triton::Triton();
    He3     =G4He3::He3();
    He4     =G4Alpha::Alpha();

    ANTIproton=G4AntiProton::AntiProton();
    ANTIneutron=G4AntiNeutron::AntiNeutron();

    ANTIdeuteron=G4AntiDeuteron::AntiDeuteron();
    ANTItriton  =G4AntiTriton::AntiTriton();
    ANTIHe3     =G4AntiHe3::AntiHe3();
    ANTIHe4     =G4AntiAlpha::AntiAlpha();

    if(preModel) { SetDeExcitation(preModel); }
    else  { 
      G4HadronicInteraction* hadi =  
    G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
      G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(hadi);
      if(!pre) { pre = new G4PreCompoundModel(); }
      SetDeExcitation(pre); 
    }
}

G4GeneratorPrecompoundInterface::~G4GeneratorPrecompoundInterface ( )

virtual

Definition at line 102 of file G4GeneratorPrecompoundInterface.cc.

{
}

Member Function Documentation

G4HadFinalState * G4GeneratorPrecompoundInterface::ApplyYourself ( const G4HadProjectile &  aTrack, G4Nucleus &  targetNucleus  )

virtual

Implements G4HadronicInteraction.

Definition at line 264 of file G4GeneratorPrecompoundInterface.cc.

References G4cout, and G4endl.

{
    G4cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."
            << G4endl;
    G4cout << "This class is only a mediator between generator and precompound"<<G4endl;
    G4cout << "Please remove from your physics list."<<G4endl;
    throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
    return new G4HadFinalState;
}

G4ReactionProductVector * G4GeneratorPrecompoundInterface::Propagate ( G4KineticTrackVector *  theSecondaries, G4V3DNucleus *  theNucleus  )

virtual

Implements G4VIntraNuclearTransportModel.

Definition at line 113 of file G4GeneratorPrecompoundInterface.cc.

References G4Nucleon::AreYouHit(), G4INCL::ClusterDecay::decay(), G4VPreCompoundModel::DeExcite(), CLHEP::HepLorentzVector::e(), python.hepunit::eplus, G4cout, G4endl, G4HadProjectile::Get4Momentum(), G4Nucleon::Get4Momentum(), G4Nucleon::GetBindingEnergy(), G4V3DNucleus::GetCharge(), G4Nucleon::GetDefinition(), G4V3DNucleus::GetMassNumber(), G4V3DNucleus::GetNextNucleon(), G4NucleiProperties::GetNuclearMass(), G4V3DNucleus::GetNuclearRadius(), G4ParticleDefinition::GetPDGCharge(), G4VIntraNuclearTransportModel::GetPrimaryProjectile(), CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), CLHEP::HepLorentzVector::setE(), G4ReactionProduct::SetMomentum(), G4Fragment::SetNumberOfCharged(), G4Fragment::SetNumberOfHoles(), G4Fragment::SetNumberOfParticles(), G4ReactionProduct::SetTotalEnergy(), G4V3DNucleus::StartLoop(), G4VIntraNuclearTransportModel::theDeExcitation, and CLHEP::HepLorentzVector::vect().

{
    G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;

    // decay the strong resonances
    G4DecayKineticTracks decay(theSecondaries);

    // prepare the fragment
    G4int anA=theNucleus->GetMassNumber();
    G4int aZ=theNucleus->GetCharge();
    G4int numberOfEx = 0;
    G4int numberOfCh = 0;
    G4int numberOfHoles = 0;
    G4double exEnergy = 0.0;
    G4double R = theNucleus->GetNuclearRadius();
    G4ThreeVector exciton3Momentum(0.,0.,0.);
    G4ThreeVector captured3Momentum(0.,0.,0.);
    G4ThreeVector wounded3Momentum(0.,0.,0.);

    // loop over secondaries
#ifdef exactExcitationEnergy
    G4LorentzVector secondary4Momemtum(0,0,0,0);
#endif
    G4KineticTrackVector::iterator iter;
    for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
    {
        G4ParticleDefinition* part = (*iter)->GetDefinition();
        G4double e = (*iter)->Get4Momentum().e();
        G4double mass = (*iter)->Get4Momentum().mag();
        G4ThreeVector mom = (*iter)->Get4Momentum().vect();
        if((part != proton && part != neutron) ||
                (e > mass + CaptureThreshold) ||
                ((*iter)->GetPosition().mag() > R)) {
            G4ReactionProduct * theNew = new G4ReactionProduct(part);
            theNew->SetMomentum(mom);
            theNew->SetTotalEnergy(e);
            theTotalResult->push_back(theNew);
#ifdef exactExcitationEnergy
            secondary4Momemtum += (*iter)->Get4Momentum();
#endif
        } else {
            // within the nucleus, neutron or proton
            // now calculate  A, Z of the fragment, momentum, number of exciton states
            ++anA;
            ++numberOfEx;
            G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
            aZ += Z;
            numberOfCh += Z;
            captured3Momentum += mom;
            exEnergy += (e - mass);
        }
        delete (*iter);
    }
    delete theSecondaries;

    // loop over wounded nucleus
    G4Nucleon * theCurrentNucleon =
            theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
    while(theCurrentNucleon) {
        if(theCurrentNucleon->AreYouHit()) {
            ++numberOfHoles;
            ++numberOfEx;
            --anA;
            aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/eplus + 0.1);
            wounded3Momentum += theCurrentNucleon->Get4Momentum().vect();
            //G4cout << "hit nucleon " << theCurrentNucleon->Get4Momentum() << G4endl;
            exEnergy += theCurrentNucleon->GetBindingEnergy();
        }
        theCurrentNucleon = theNucleus->GetNextNucleon();
    }
    exciton3Momentum = captured3Momentum - wounded3Momentum;

    if(anA == 0) return theTotalResult;

    if(anA >= aZ) 
    {
        G4double fMass =  G4NucleiProperties::GetNuclearMass(anA, aZ);

#ifdef exactExcitationEnergy
       // recalculate exEnergy from Energy balance....
        const G4HadProjectile * primary = GetPrimaryProjectile();
        G4double Einitial= primary->Get4Momentum().e()
                + G4NucleiProperties::GetNuclearMass(theNucleus->GetMassNumber(),
                                                     theNucleus->GetCharge());
// Uzhi        G4double Efinal = fMass + secondary4Momemtum.e();
        G4double Efinal = std::sqrt(exciton3Momentum.mag2() + fMass*fMass)
                        + secondary4Momemtum.e();
        if ( (Einitial - Efinal) > 0 ) {
            // G4cout << "G4GPI::Propagate() : positive exact excitation Energy "
            //        << (Einitial - Efinal)/MeV << " MeV, exciton estimate " 
            //        << exEnergy/MeV << " MeV" << G4endl;

//          exEnergy=Einitial - Efinal;
            G4LorentzVector PrimMom=primary->Get4Momentum(); PrimMom.setE(Einitial);

            exEnergy=(PrimMom - secondary4Momemtum).mag() - fMass;
        }
        else {
            //  G4cout << "G4GeneratorPrecompoundInterface::Propagate() : " 
            //         << "negative exact excitation Energy "
            //         << (Einitial - Efinal)/MeV 
            //         << " MeV, setting  excitation to 0 MeV" << G4endl;
            exEnergy=0.;
        }
#endif

        if(exEnergy < 0.) exEnergy=0.;   // Uzhi 11 Dec. 2012

        fMass += exEnergy;

        G4ThreeVector balance=primary->Get4Momentum().vect() - 
                              secondary4Momemtum.vect() - exciton3Momentum;

#ifdef G4GPI_debug_excitation
        G4cout << "momentum balance" << balance 
               << " value " << balance.mag()                   <<G4endl
               << "primary         "<< primary->Get4Momentum() <<G4endl
               << "secondary       "<< secondary4Momemtum      <<G4endl
               << "captured        "<< captured3Momentum       <<G4endl
               << "wounded         "<< wounded3Momentum        <<G4endl
               << "exciton         "<< exciton3Momentum        <<G4endl
               << "second + exciton" 
               << secondary4Momemtum.vect() + exciton3Momentum << G4endl;
#endif
//#ifdef exactExcitationEnergy
//        G4LorentzVector exciton4Momentum(exciton3Momentum, fMass);
//        G4LorentzVector exciton4Momentum(exciton3Momentum,
//                std::sqrt(exciton3Momentum.mag2() + fMass*fMass));
//#else
        G4LorentzVector exciton4Momentum(exciton3Momentum,
                std::sqrt(exciton3Momentum.mag2() + fMass*fMass));
//#endif
//G4cout<<"exciton4Momentum "<<exciton4Momentum<<G4endl;
        // Need to de-excite the remnant nucleus only if excitation energy > 0.
        G4Fragment anInitialState(anA, aZ, exciton4Momentum);
        anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
        anInitialState.SetNumberOfCharged(numberOfCh);
        anInitialState.SetNumberOfHoles(numberOfHoles);

        G4ReactionProductVector * aPrecoResult = 
                                  theDeExcitation->DeExcite(anInitialState);
        // fill pre-compound part into the result, and return
        theTotalResult->insert(theTotalResult->end(),aPrecoResult->begin(),
                                            aPrecoResult->end() );
        delete aPrecoResult;
    }

    return theTotalResult;
}

void G4GeneratorPrecompoundInterface::PropagateModelDescription ( std::ostream &  outFile ) const

virtual

Reimplemented from G4VIntraNuclearTransportModel.

Definition at line 273 of file G4GeneratorPrecompoundInterface.cc.

{
    outFile << "G4GeneratorPrecompoundInterface interfaces a high\n"
            << "energy model through the wounded nucleus to precompound de-excition.\n"
            << "Low energy protons and neutron present among secondaries produced by \n"
            << "the high energy generator and within the nucleus are captured. The wounded\n"
            << "nucleus and the captured particles form an excited nuclear fragment. This\n"
            << "fragment is passed to the Geant4 pre-compound model for de-excitation.\n"
            << "Nuclear de-excitation:\n";
    // preco

}

G4ReactionProductVector * G4GeneratorPrecompoundInterface::PropagateNuclNucl ( G4KineticTrackVector *  theSecondaries, G4V3DNucleus *  theNucleus, G4V3DNucleus *  theProjectileNucleus  )

virtual

Reimplemented from G4VIntraNuclearTransportModel.

Definition at line 288 of file G4GeneratorPrecompoundInterface.cc.

References G4AntiAlpha::AntiAlphaDefinition(), G4AntiDeuteron::AntiDeuteronDefinition(), G4AntiHe3::AntiHe3Definition(), G4AntiNeutron::AntiNeutronDefinition(), G4AntiProton::AntiProtonDefinition(), G4AntiTriton::AntiTritonDefinition(), G4Nucleon::AreYouHit(), CLHEP::HepLorentzVector::boost(), CLHEP::HepLorentzVector::boostVector(), G4INCL::ClusterDecay::decay(), G4VPreCompoundModel::DeExcite(), CLHEP::HepLorentzVector::e(), python.hepunit::eplus, G4cout, G4endl, G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleon::Get4Momentum(), G4ParticleDefinition::GetBaryonNumber(), G4Nucleon::GetBindingEnergy(), G4V3DNucleus::GetCharge(), G4HadProjectile::GetDefinition(), G4Nucleon::GetDefinition(), G4ReactionProduct::GetDefinition(), G4V3DNucleus::GetMassNumber(), G4V3DNucleus::GetNextNucleon(), G4NucleiProperties::GetNuclearMass(), G4V3DNucleus::GetNuclearRadius(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGCharge(), G4VIntraNuclearTransportModel::GetPrimaryProjectile(), CLHEP::HepLorentzVector::mag(), CLHEP::Hep3Vector::mag2(), CLHEP::HepLorentzVector::setE(), G4ReactionProduct::SetMomentum(), G4Fragment::SetNumberOfCharged(), G4Fragment::SetNumberOfHoles(), G4Fragment::SetNumberOfParticles(), G4ReactionProduct::SetTotalEnergy(), G4V3DNucleus::StartLoop(), G4VIntraNuclearTransportModel::theDeExcitation, and CLHEP::HepLorentzVector::vect().

{
#ifdef debugPrecoInt
  G4cout<<"G4GeneratorPrecompoundInterface::PropagateNuclNucl "<<G4endl;
#endif

    G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;

    // prepare the target residual
    G4int anA=theNucleus->GetMassNumber();
    G4int aZ=theNucleus->GetCharge();
    G4int numberOfEx = 0;
    G4int numberOfCh = 0;
    G4int numberOfHoles = 0;
    G4double exEnergy = 0.0;
    G4double R = theNucleus->GetNuclearRadius();
    G4LorentzVector Target4Momentum(0,0,0,0);

#ifdef debugPrecoInt
  G4cout<<"Target A Z "<<anA<<" "<<aZ<<G4endl;
#endif

    // loop over wounded target nucleus
    G4Nucleon * theCurrentNucleon =
            theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
    while(theCurrentNucleon) {
        if(theCurrentNucleon->AreYouHit()) {
            ++numberOfHoles;
            ++numberOfEx;
            --anA;
            aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
                                                                        eplus + 0.1);
            exEnergy += theCurrentNucleon->GetBindingEnergy();
                        Target4Momentum -=theCurrentNucleon->Get4Momentum();
        } 
        theCurrentNucleon = theNucleus->GetNextNucleon();
    }

#ifdef debugPrecoInt
  G4cout<<"Residual Target A Z E* 4mom "<<anA<<" "<<aZ<<" "<<exEnergy<<" "
        <<Target4Momentum<<G4endl;
#endif

    // prepare the projectile residual
#ifdef debugPrecoInt
  G4cout<<"Primary BaryonNumber "
        <<GetPrimaryProjectile()->GetDefinition()->GetBaryonNumber()<<G4endl;
#endif

        G4bool ProjectileIsAntiNucleus=
          GetPrimaryProjectile()->GetDefinition()->GetBaryonNumber() < -1;

        G4ThreeVector bst = GetPrimaryProjectile()->Get4Momentum().boostVector();

    G4int anAb=theProjectileNucleus->GetMassNumber();
    G4int aZb=theProjectileNucleus->GetCharge();
    G4int numberOfExB = 0;
    G4int numberOfChB = 0;
    G4int numberOfHolesB = 0;
    G4double exEnergyB = 0.0;
    G4double Rb = theProjectileNucleus->GetNuclearRadius();
    G4LorentzVector Projectile4Momentum(0,0,0,0);

#ifdef debugPrecoInt
  G4cout<<"Projectile A Z "<<anAb<<" "<<aZb<<G4endl;
#endif

    // loop over wounded projectile nucleus
    theCurrentNucleon =
      theProjectileNucleus->StartLoop() ? theProjectileNucleus->GetNextNucleon() : 0;
    while(theCurrentNucleon) {
        if(theCurrentNucleon->AreYouHit()) {
            ++numberOfHolesB;
            ++numberOfExB;
            --anAb;
                     if(!ProjectileIsAntiNucleus)
                     {
            aZb -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
                                                                         eplus + 0.1);
                     } else
                     {
            aZb += G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/
                                                                         eplus - 0.1);
                     }
            exEnergyB += theCurrentNucleon->GetBindingEnergy();
                        Projectile4Momentum -=theCurrentNucleon->Get4Momentum();
        } 
        theCurrentNucleon = theProjectileNucleus->GetNextNucleon();
    }

        G4bool ExistTargetRemnant   =  G4double (numberOfHoles)        <
                                  0.3* G4double (numberOfHoles + anA);  
        G4bool ExistProjectileRemnant= G4double (numberOfHolesB)       <
                                   0.3*G4double (numberOfHolesB + anAb);  

#ifdef debugPrecoInt
  G4cout<<"Projectile residual A Z E* 4mom "<<anAb<<" "<<aZb<<" "<<exEnergyB<<" "
        <<Projectile4Momentum<<G4endl;
  G4cout<<" ExistTargetRemnant ExistProjectileRemnant "
        <<ExistTargetRemnant<<" "<< ExistProjectileRemnant<<G4endl;
#endif
//-----------------------------------------------------------------------------
    // decay the strong resonances
    G4DecayKineticTracks decay(theSecondaries);

#ifdef debugPrecoInt
  G4LorentzVector secondary4Momemtum(0,0,0,0);
  G4int SecondrNum(0);
#endif

    // loop over secondaries
        G4KineticTrackVector::iterator iter;
        for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
    {
        G4ParticleDefinition* part = (*iter)->GetDefinition();
            G4LorentzVector aTrack4Momentum=(*iter)->Get4Momentum();

        if( part != proton && part != neutron &&
                   (part != ANTIproton  && ProjectileIsAntiNucleus) && 
                   (part != ANTIneutron && ProjectileIsAntiNucleus)   )  
                {
         G4ReactionProduct * theNew = new G4ReactionProduct(part);
         theNew->SetMomentum(aTrack4Momentum.vect());
         theNew->SetTotalEnergy(aTrack4Momentum.e());
         theTotalResult->push_back(theNew);
#ifdef debugPrecoInt
  SecondrNum++;
  secondary4Momemtum += (*iter)->Get4Momentum();
  G4cout<<"Secondary  "<<SecondrNum<<" "
        <<theNew->GetDefinition()->GetParticleName()<<" "
        <<secondary4Momemtum<<G4endl;
#endif
         delete (*iter);
                 continue;
                }

                G4bool CanBeCapturedByTarget = false;
        if( part == proton || part == neutron)
                {
                   CanBeCapturedByTarget = ExistTargetRemnant    &&
                                 (CaptureThreshold >
                   (aTrack4Momentum       + Target4Momentum).mag() -
                    aTrack4Momentum.mag() - Target4Momentum.mag())   &&
                             ((*iter)->GetPosition().mag() < R);
                }
// ---------------------------
                G4LorentzVector Position((*iter)->GetPosition(),
                                        (*iter)->GetFormationTime());
                Position.boost(bst);

                G4bool CanBeCapturedByProjectile = false;

        if( !ProjectileIsAntiNucleus && 
                    ( part == proton || part == neutron))
                {
                   CanBeCapturedByProjectile = ExistProjectileRemnant &&
                                 (CaptureThreshold >
                   (aTrack4Momentum       + Projectile4Momentum).mag() -
                    aTrack4Momentum.mag() - Projectile4Momentum.mag())    &&
                             (Position.vect().mag() < Rb);
                }

        if( ProjectileIsAntiNucleus && 
                    ( part == ANTIproton || part == ANTIneutron))
                {
                   CanBeCapturedByProjectile = ExistProjectileRemnant &&
                                 (CaptureThreshold >
                   (aTrack4Momentum       + Projectile4Momentum).mag() -
                    aTrack4Momentum.mag() - Projectile4Momentum.mag())    &&
                             (Position.vect().mag() < Rb);
                }

                if(CanBeCapturedByTarget && CanBeCapturedByProjectile)
                {
                 if(G4UniformRand() < 0.5) 
                 { CanBeCapturedByTarget = true; CanBeCapturedByProjectile = false;} 
                 else
                 { CanBeCapturedByTarget = false; CanBeCapturedByProjectile = true;}
                }

                if(CanBeCapturedByTarget)
                {
             // within the target nucleus, neutron or proton
         // now calculate  A, Z of the fragment, momentum, 
                 // number of exciton states
#ifdef debugPrecoInt
  G4cout<<"Track is CapturedByTarget "<<" "
        <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
#endif
         ++anA;
         ++numberOfEx;
         G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
         aZ += Z;
         numberOfCh += Z;
                 Target4Momentum +=aTrack4Momentum;
         delete (*iter);
        } else if(CanBeCapturedByProjectile)
                {
             // within the projectile nucleus, neutron or proton
         // now calculate  A, Z of the fragment, momentum, 
                 // number of exciton states
#ifdef debugPrecoInt
 G4cout<<"Track is CapturedByProjectile"<<" "
       <<aTrack4Momentum<<" "<<aTrack4Momentum.mag()<<G4endl;
#endif
         ++anAb;
         ++numberOfExB;
         G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
                 if( ProjectileIsAntiNucleus ) Z=-Z;
         aZb += Z;
         numberOfChB += Z;
                 Projectile4Momentum +=aTrack4Momentum;
         delete (*iter);
                } else
                { // the track is not captured
         G4ReactionProduct * theNew = new G4ReactionProduct(part);
         theNew->SetMomentum(aTrack4Momentum.vect());
         theNew->SetTotalEnergy(aTrack4Momentum.e());
         theTotalResult->push_back(theNew);

#ifdef debugPrecoInt
  SecondrNum++;
  secondary4Momemtum += (*iter)->Get4Momentum();
  G4cout<<"Secondary  "<<SecondrNum<<" "
        <<theNew->GetDefinition()->GetParticleName()<<" "
        <<secondary4Momemtum<<G4endl;
#endif
         delete (*iter);
                 continue;
                }
    }
    delete theSecondaries;
//-----------------------------------------------------

#ifdef debugPrecoInt
 G4cout<<"Final target residual A Z E* 4mom "<<anA<<" "<<aZ<<" "
       <<exEnergy<<" "<<Target4Momentum<<G4endl;
#endif

    if(0!=anA ) 
        {
         G4double fMass =  G4NucleiProperties::GetNuclearMass(anA, aZ);

         if((anA == theNucleus->GetMassNumber()) && (exEnergy <= 0.))
         {Target4Momentum.setE(fMass);}

         G4double RemnMass=Target4Momentum.mag();
         if(RemnMass < fMass) 
         {
          RemnMass=fMass + exEnergy;
          Target4Momentum.setE(std::sqrt(Target4Momentum.vect().mag2() +
                                                RemnMass*RemnMass));
         } else
         { exEnergy=RemnMass-fMass;}

         if( exEnergy < 0.) exEnergy=0.;

         // Need to de-excite the remnant nucleus 
         G4Fragment anInitialState(anA, aZ, Target4Momentum);
     anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
     anInitialState.SetNumberOfCharged(numberOfCh);
     anInitialState.SetNumberOfHoles(numberOfHoles);

     G4ReactionProductVector * aPrecoResult = 
                                   theDeExcitation->DeExcite(anInitialState);

     // fill pre-compound part into the result, and return
     for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll) 
         {
            theTotalResult->push_back(aPrecoResult->operator[](ll));
#ifdef debugPrecoInt
 G4cout<<"Tr frag "<<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()
       <<" "<<aPrecoResult->operator[](ll)->GetMomentum()<<G4endl;
#endif
     }
     delete aPrecoResult;
    }

//-----------------------------------------------------
        if((anAb == theProjectileNucleus->GetMassNumber())&& (exEnergyB <= 0.))
        {Projectile4Momentum = GetPrimaryProjectile()->Get4Momentum();}

#ifdef debugPrecoInt
  G4cout<<"Final projectile residual A Z E* Pmom "<<anAb<<" "<<aZb<<" "
        <<exEnergyB<<" "<<Projectile4Momentum<<G4endl;
#endif

    if(0!=anAb) 
        {
     G4double fMass =  G4NucleiProperties::GetNuclearMass(anAb, aZb);
         G4double RemnMass=Projectile4Momentum.mag();

         if(RemnMass < fMass) 
         {
          RemnMass=fMass + exEnergyB;
          Projectile4Momentum.setE(std::sqrt(Projectile4Momentum.vect().mag2() +
                                                    RemnMass*RemnMass));
         } else
         { exEnergyB=RemnMass-fMass;}

         if( exEnergyB < 0.) exEnergyB=0.;

         // Need to de-excite the remnant nucleus 
         G4Fragment anInitialState(anAb, aZb, Projectile4Momentum);
     anInitialState.SetNumberOfParticles(numberOfExB-numberOfHolesB);
     anInitialState.SetNumberOfCharged(numberOfChB);
     anInitialState.SetNumberOfHoles(numberOfHolesB);

     G4ReactionProductVector * aPrecoResult = 
                                   theDeExcitation->DeExcite(anInitialState);

            // fill pre-compound part into the result, and return
     for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll) 
         {
          if(ProjectileIsAntiNucleus)
          {

#ifdef debugPrecoInt
  G4cout<<"aPrecoRes  "<<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()
        <<" "<<aPrecoResult->operator[](ll)->GetMomentum()
        <<" "<<aPrecoResult->operator[](ll)->GetTotalEnergy()
        <<" "<<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
#endif

      G4ParticleDefinition * aFragment=aPrecoResult->operator[](ll)->GetDefinition();
      G4ParticleDefinition * LastFragment=aFragment;
      if     (aFragment == proton)  {LastFragment=G4AntiProton::AntiProtonDefinition();}
      else if(aFragment == neutron) {LastFragment=G4AntiNeutron::AntiNeutronDefinition();}
      else if(aFragment == deuteron){LastFragment=G4AntiDeuteron::AntiDeuteronDefinition();}
      else if(aFragment == triton)  {LastFragment=G4AntiTriton::AntiTritonDefinition();}
      else if(aFragment == He3)     {LastFragment=G4AntiHe3::AntiHe3Definition();}
      else if(aFragment == He4)     {LastFragment=G4AntiAlpha::AntiAlphaDefinition();}
      else {}

           aPrecoResult->operator[](ll)->SetDefinitionAndUpdateE(LastFragment);
          }

#ifdef debugPrecoInt
  G4cout<<"aPrecoResA "<<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()
        <<" "<<aPrecoResult->operator[](ll)->GetMomentum()
        <<" "<<aPrecoResult->operator[](ll)->GetTotalEnergy()
        <<" "<<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
#endif
      theTotalResult->push_back(aPrecoResult->operator[](ll));
     }

         delete aPrecoResult;
    }

        return theTotalResult;
}

void G4GeneratorPrecompoundInterface::SetCaptureThreshold ( G4double  value )

inline

Definition at line 107 of file G4GeneratorPrecompoundInterface.hh.

{
  CaptureThreshold=value;
}

---

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

• G4GeneratorPrecompoundInterface.hh
• G4GeneratorPrecompoundInterface.cc