#include <G4GeneratorPrecompoundInterface.hh>

Inheritance diagram for G4GeneratorPrecompoundInterface:

Collaboration 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 &mName="CascadeModel", G4VPreCompoundModel *ptr=nullptr)

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 &aTrack, G4Nucleus &targetNucleus)

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)

G4int	GetVerboseLevel () const

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

virtual const std::pair < G4double, G4double >	GetFatalEnergyCheckLevels () const

virtual std::pair< G4double, G4double >	GetEnergyMomentumCheckLevels () const

void	SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)

virtual void	BuildPhysicsTable (const G4ParticleDefinition &)

virtual void	InitialiseModel ()

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 79 of file G4GeneratorPrecompoundInterface.cc.

 : CaptureThreshold(70*MeV)   // Uzhi 1.05.2015 10 ->70
 {
    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);
    }
 }

Here is the call graph for this function:

G4GeneratorPrecompoundInterface::~G4GeneratorPrecompoundInterface ( )

virtual

Definition at line 108 of file G4GeneratorPrecompoundInterface.cc.

109 {

110 }

Member Function Documentation

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

virtual

Implements G4HadronicInteraction.

Definition at line 326 of file G4GeneratorPrecompoundInterface.cc.

 {
    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.

 {
    #ifdef debugPrecoInt
       G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::Propagate"<<G4endl;
       G4cout<<"Target A and Z "<<theNucleus->GetMassNumber()<<" "<<theNucleus->GetCharge()<<G4endl;
       G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
    #endif
 
    G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
 
    // decay the strong resonances
    G4DecayKineticTracks decay(theSecondaries);
    #ifdef debugPrecoInt
       G4cout<<"Final stable particles number "<<theSecondaries->size()<<G4endl;
    #endif
 
    // prepare the fragment
    G4int anA=theNucleus->GetMassNumber();
    G4int aZ=theNucleus->GetCharge();
 // G4double TargetNucleusMass =  G4NucleiProperties::GetNuclearMass(anA, aZ);
 
    G4int numberOfEx = 0;
    G4int numberOfCh = 0;
    G4int numberOfHoles = 0;
 
    G4double R = theNucleus->GetNuclearRadius();
 
    G4LorentzVector captured4Momentum(0.,0.,0.,0.);
    G4LorentzVector Residual4Momentum(0.,0.,0.,0.);  // TargetNucleusMass is not need at the moment 
    G4LorentzVector Secondary4Momentum(0.,0.,0.,0.);
 
    // loop over secondaries
    G4KineticTrackVector::iterator iter;
    for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
    {
       const G4ParticleDefinition* part = (*iter)->GetDefinition();
       G4double e = (*iter)->Get4Momentum().e();
       G4double mass = (*iter)->Get4Momentum().mag();
       G4ThreeVector mom = (*iter)->Get4Momentum().vect();
       if((part != proton && part != neutron) ||
 // Uzhi 2.05.2015            (e > mass + CaptureThreshold) ||
             ((*iter)->GetPosition().mag() > R)) {
          G4ReactionProduct * theNew = new G4ReactionProduct(part);
          theNew->SetMomentum(mom);
          theNew->SetTotalEnergy(e);
          theTotalResult->push_back(theNew);
          Secondary4Momentum += (*iter)->Get4Momentum();                 // Uzhi 29 April
          #ifdef debugPrecoInt
             G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
                   <<(*iter)->Get4Momentum().mag()<<G4endl;
          #endif
       } else {
          if( e-mass > -CaptureThreshold*G4Log( G4UniformRand()) ) {  // Added by Uzhi 2.05.2015
             G4ReactionProduct * theNew = new G4ReactionProduct(part);
             theNew->SetMomentum(mom);
             theNew->SetTotalEnergy(e);
             theTotalResult->push_back(theNew);
             Secondary4Momentum += (*iter)->Get4Momentum();              // Uzhi 29 April
             #ifdef debugPrecoInt
                G4cout<<"Secondary 4Mom "<<part->GetParticleName()<<" "<<(*iter)->Get4Momentum()<<" "
                      <<(*iter)->Get4Momentum().mag()<<G4endl;
             #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;
             captured4Momentum += (*iter)->Get4Momentum();
             #ifdef debugPrecoInt
                G4cout<<"Captured  4Mom "<<part->GetParticleName()<<(*iter)->Get4Momentum()<<G4endl;
             #endif
          }
       }
       delete (*iter);
    }
    delete theSecondaries;
 
    // loop over wounded nucleus
    G4Nucleon * theCurrentNucleon =
          theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
    while(theCurrentNucleon)    /* Loop checking, 31.08.2015, G.Folger */
     {
       if(theCurrentNucleon->AreYouHit()) {
          ++numberOfHoles;
          ++numberOfEx;
          --anA;
          aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/eplus + 0.1);
 
          Residual4Momentum -= theCurrentNucleon->Get4Momentum();
       }
       theCurrentNucleon = theNucleus->GetNextNucleon();
    }
 
    #ifdef debugPrecoInt
       G4cout<<G4endl;
       G4cout<<"Secondary 4Mom "<<Secondary4Momentum<<G4endl;
       G4cout<<"Captured  4Mom "<<captured4Momentum<<G4endl;
       G4cout<<"Sec + Captured "<<Secondary4Momentum+captured4Momentum<<G4endl;
       G4cout<<"Residual4Mom   "<<Residual4Momentum<<G4endl;
       G4cout<<"Sum 4 momenta  "
             <<Secondary4Momentum + captured4Momentum + Residual4Momentum <<G4endl;
    #endif
 
    // Check that we use QGS model; loop over wounded nucleus
    G4bool QGSM(false);
    theCurrentNucleon = theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
    while(theCurrentNucleon)   /* Loop checking, 31.08.2015, G.Folger */
     {
       if(theCurrentNucleon->AreYouHit()) 
       {
        if(theCurrentNucleon->Get4Momentum().mag() < 
           theCurrentNucleon->GetDefinition()->GetPDGMass()) QGSM=true;
       }
       theCurrentNucleon = theNucleus->GetNextNucleon();
    }
 
    #ifdef debugPrecoInt
       if(!QGSM){
         G4cout<<G4endl;
         G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
         G4cout<<"Residual  4Mom "<<Residual4Momentum<<G4endl;
         if(numberOfEx == 0) 
         {G4cout<<"Residual  4Mom = 0 means that there were not wounded and captured nucleons"<<G4endl;}
       }
    #endif
 
    if(anA == 0) return theTotalResult;
 
    G4LorentzVector exciton4Momentum(0.,0.,0.,0.);                       // Uzhi 29 April
    if(anA >= aZ)
    {
     if(!QGSM)
     {          // FTF model was used  
       G4double fMass =  G4NucleiProperties::GetNuclearMass(anA, aZ);
 
 //      G4LorentzVector exciton4Momentum = Residual4Momentum + captured4Momentum;
       exciton4Momentum = Residual4Momentum + captured4Momentum;
 //exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass)));
       G4double ActualMass = exciton4Momentum.mag();      
       if(ActualMass <= fMass ) {                             //E*<=0,  Uzhi 5.05.2015
         exciton4Momentum.setE(std::sqrt(exciton4Momentum.vect().mag2()+sqr(fMass))); // Uzhi 13.05.2015
       }
 
       #ifdef debugPrecoInt
          G4double exEnergy = 0.0;
          if(ActualMass <= fMass ) {exEnergy = 0.;}                      // Uzhi 5.05.2015
          else                     {exEnergy = ActualMass - fMass;}
          G4cout<<"Ground state residual Mass "<<fMass<<" E* "<<exEnergy<<G4endl;
       #endif
     }
     else
     {          // QGS model was used 
      G4double InitialTargetMass = 
               G4NucleiProperties::GetNuclearMass(theNucleus->GetMassNumber(), theNucleus->GetCharge());
 
      exciton4Momentum = 
               GetPrimaryProjectile()->Get4Momentum() + G4LorentzVector(0.,0.,0.,InitialTargetMass)
              -Secondary4Momentum;
 
      G4double fMass =  G4NucleiProperties::GetNuclearMass(anA, aZ);
      G4double ActualMass = exciton4Momentum.mag();
 
      #ifdef debugPrecoInt
         G4cout<<G4endl;
         G4cout<<"Residual A and Z "<<anA<<" "<<aZ<<G4endl;
         G4cout<<"Residual4Momentum "<<exciton4Momentum<<G4endl;
         G4cout<<"ResidualMass, GroundStateMass and E* "<<ActualMass<<" "<<fMass<<" "
               <<ActualMass - fMass<<G4endl;
      #endif
 
      if(ActualMass - fMass < 0.)
      {
       G4double ResE = std::sqrt(exciton4Momentum.vect().mag2() + sqr(fMass+10*MeV));
       exciton4Momentum.setE(ResE);
       #ifdef debugPrecoInt
          G4cout<<"ActualMass - fMass < 0. "<<ActualMass<<" "<<fMass<<" "<<ActualMass - fMass<<G4endl;
          G4int Uzhi; G4cin>>Uzhi;
       #endif
      }
     }
       // 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
       #ifdef debugPrecoInt
          G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
       #endif
       for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
       {
          theTotalResult->push_back(aPrecoResult->operator[](ll));
          #ifdef debugPrecoInt
             G4cout<<"Fragment "<<ll<<" "
                   <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
                   <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
                   <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
                   <<aPrecoResult->operator[](ll)->GetDefinition()->GetPDGMass()<<G4endl;
          #endif
       }
       delete aPrecoResult;
    }
 
    return theTotalResult;
 }

Here is the call graph for this function:

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

virtual

Reimplemented from G4VIntraNuclearTransportModel.

Definition at line 335 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 350 of file G4GeneratorPrecompoundInterface.cc.

 {
 #ifdef debugPrecoInt
    G4cout<<G4endl<<"G4GeneratorPrecompoundInterface::PropagateNuclNucl "<<G4endl;
    G4cout<<"Projectile A and Z "<<theProjectileNucleus->GetMassNumber()<<" "
                                 <<theProjectileNucleus->GetCharge()<<G4endl;
    G4cout<<"Target     A and Z "<<theNucleus->GetMassNumber()<<" "
                                 <<theNucleus->GetCharge()<<G4endl;
    G4cout<<"Directly produced particles number "<<theSecondaries->size()<<G4endl;
    G4cout<<"Projectile 4Mom and mass "<<GetPrimaryProjectile()->Get4Momentum()<<" "
                                       <<GetPrimaryProjectile()->Get4Momentum().mag()<<G4endl<<G4endl;
 #endif
 
    // 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.);
 
    // loop over wounded target nucleus
    G4Nucleon * theCurrentNucleon =
          theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
    while(theCurrentNucleon)   /* Loop checking, 31.08.2015, G.Folger */
     {
       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
 
    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.);
 
    // loop over wounded projectile nucleus
    theCurrentNucleon =
          theProjectileNucleus->StartLoop() ? theProjectileNucleus->GetNextNucleon() : 0;
    while(theCurrentNucleon)    /* Loop checking, 31.08.2015, G.Folger */
     {
       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
    G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
    G4DecayKineticTracks decay(theSecondaries);
    #ifdef debugPrecoInt
       G4cout<<"Secondary stable particles number "<<theSecondaries->size()<<G4endl;
    #endif
 
 #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)
    {
       const 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()<<" "
                <<theNew->GetMomentum()<<" "<<theNew->GetTotalEnergy()<<G4endl;
 
 #endif
          delete (*iter);
          continue;
       }
 
       G4bool CanBeCapturedByTarget = false;
       if( part == proton || part == neutron)
       {
          CanBeCapturedByTarget = ExistTargetRemnant    &&
                (-CaptureThreshold*G4Log( G4UniformRand()) >
               (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*G4Log( G4UniformRand()) >
               (aTrack4Momentum       + Projectile4Momentum).mag() -
                aTrack4Momentum.mag() - Projectile4Momentum.mag())    &&
                    (Position.vect().mag() < Rb);
       }
 
       if( ProjectileIsAntiNucleus &&
             ( part == ANTIproton || part == ANTIneutron))
       {
          CanBeCapturedByProjectile = ExistProjectileRemnant &&
                (-CaptureThreshold*G4Log( G4UniformRand()) >
               (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 "<<" "<<part->GetParticleName()<<" "
                <<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"<<" "<<part->GetParticleName()<<" "
                <<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);
 
       #ifdef debugPrecoInt
          G4cout<<"Target fragment number "<<aPrecoResult->size()<<G4endl;
       #endif
 
       // 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<<"Target fragment "<<ll<<" "
                   <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
                   <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
                   <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
                   <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
          #endif
       }
       delete aPrecoResult;
    }
 
    //-----------------------------------------------------
    if((anAb == theProjectileNucleus->GetMassNumber())&& (exEnergyB <= 0.))
    {Projectile4Momentum = GetPrimaryProjectile()->Get4Momentum();}
 
    #ifdef debugPrecoInt
       G4cout<<"Final projectile residual A Z E* Pmom Pmag2 "<<anAb<<" "<<aZb<<" "
             <<exEnergyB<<" "<<Projectile4Momentum<<" "
                             <<Projectile4Momentum.mag2()<<G4endl;
    #endif
 
    if(0!=anAb)
    {
   //    G4ThreeVector bstToCM =Projectile4Momentum.findBoostToCM();     // Uzhi Apr. 2015
   //    Projectile4Momentum.boost(bstToCM);                             // Uzhi Apr. 2015
 
       G4double fMass =  G4NucleiProperties::GetNuclearMass(anAb, aZb);
       G4double RemnMass=Projectile4Momentum.mag();
 
       if(RemnMass < fMass)
       {
          RemnMass=fMass + exEnergyB;
          Projectile4Momentum.setE(std::sqrt(Projectile4Momentum.vect().mag2() + // Uzhi 8.05.2015
                RemnMass*RemnMass));                                             // Uzhi 8.05.2015
       } else
       { exEnergyB=RemnMass-fMass;}
 
       if( exEnergyB < 0.) exEnergyB=0.;
 
       G4ThreeVector bstToCM =Projectile4Momentum.findBoostToCM();       // Uzhi Apr. 2015
       Projectile4Momentum.boost(bstToCM);                               // Uzhi Apr. 2015
 
       // 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);
 
       #ifdef debugPrecoInt
          G4cout<<"Projectile fragment number "<<aPrecoResult->size()<<G4endl;
       #endif
 
       // fill pre-compound part into the result, and return
       for(unsigned int ll=0; ll<aPrecoResult->size(); ++ll)
       {
          G4LorentzVector tmp=G4LorentzVector(aPrecoResult->operator[](ll)->GetMomentum(),    // Uzhi 2015
                                              aPrecoResult->operator[](ll)->GetTotalEnergy());// Uzhi 2015
          tmp.boost(-bstToCM); // Transformation to the system of original remnant            // Uzhi 2015
          aPrecoResult->operator[](ll)->SetMomentum(tmp.vect());                              // Uzhi 2015
          aPrecoResult->operator[](ll)->SetTotalEnergy(tmp.e());                              // Uzhi 2015
 
          if(ProjectileIsAntiNucleus)
          {
             const G4ParticleDefinition * aFragment=aPrecoResult->operator[](ll)->GetDefinition();
             const 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<<"Projectile fragment "<<ll<<" "
                   <<aPrecoResult->operator[](ll)->GetDefinition()->GetParticleName()<<" "
                   <<aPrecoResult->operator[](ll)->GetMomentum()<<" "
                   <<aPrecoResult->operator[](ll)->GetTotalEnergy()<<" "
                   <<aPrecoResult->operator[](ll)->GetMass()<<G4endl;
          #endif
 //Uzhi
          theTotalResult->push_back(aPrecoResult->operator[](ll));
       }
 
       delete aPrecoResult;
    }
 
    return theTotalResult;
 }

Here is the call graph for this function:

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:

source/geant4.10.03.p03/source/processes/hadronic/models/binary_cascade/include/G4GeneratorPrecompoundInterface.hh
source/geant4.10.03.p03/source/processes/hadronic/models/binary_cascade/src/G4GeneratorPrecompoundInterface.cc

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation