#include <G4QMDReaction.hh>

Inheritance diagram for G4QMDReaction:

Public Member Functions
	G4QMDReaction ()

	~G4QMDReaction ()

std::vector< G4QMDSystem * >	GetFinalStates ()

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

G4ExcitationHandler *	GetExcitationHandler ()

void	UnUseGEM ()

void	UseFRAG ()

void	SetTMAX (G4int i)

void	SetDT (G4double t)

void	SetEF (G4double x)

virtual void	ModelDescription (std::ostream &outFile) const

Public Member Functions inherited from G4HadronicInteraction
	G4HadronicInteraction (const G4String &modelName="HadronicModel")

virtual	~G4HadronicInteraction ()

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

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	ModelDescription (std::ostream &outFile) const

virtual void	BuildPhysicsTable (const G4ParticleDefinition &)

virtual void	InitialiseModel ()

	G4HadronicInteraction (const G4HadronicInteraction &right)=delete

const G4HadronicInteraction &	operator= (const G4HadronicInteraction &right)=delete

G4bool	operator== (const G4HadronicInteraction &right) const =delete

G4bool	operator!= (const G4HadronicInteraction &right) const =delete

Additional Inherited Members
Protected Member Functions inherited from G4HadronicInteraction
void	SetModelName (const G4String &nam)

G4bool	IsBlocked () const

void	Block ()

Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState	theParticleChange

G4int	verboseLevel

G4double	theMinEnergy

G4double	theMaxEnergy

G4bool	isBlocked

Detailed Description

Definition at line 61 of file G4QMDReaction.hh.

Constructor & Destructor Documentation

◆ G4QMDReaction()

G4QMDReaction::G4QMDReaction ( )

Definition at line 49 of file G4QMDReaction.cc.

: G4HadronicInteraction("QMDModel")
, system ( NULL )
, deltaT ( 1 ) // in fsec (c=1)
, maxTime ( 100 ) // will have maxTime-th time step
, envelopF ( 1.05 ) // 10% for Peripheral reactions
, gem ( true )
, frag ( false )
{
 
   //090331
   shenXS = new G4IonsShenCrossSection();
   //genspaXS = new G4GeneralSpaceNNCrossSection();
 
   pipElNucXS = new G4BGGPionElasticXS(G4PionPlus::PionPlus() );
   pipElNucXS->BuildPhysicsTable(*(G4PionPlus::PionPlus() ) );
 
   pimElNucXS = new G4BGGPionElasticXS(G4PionMinus::PionMinus() );
   pimElNucXS->BuildPhysicsTable(*(G4PionMinus::PionMinus() ) );
 
   pipInelNucXS = new G4BGGPionInelasticXS(G4PionPlus::PionPlus() );
   pipInelNucXS->BuildPhysicsTable(*(G4PionPlus::PionPlus() ) );
 
   pimInelNucXS = new G4BGGPionInelasticXS(G4PionMinus::PionMinus() );
   pimInelNucXS->BuildPhysicsTable(*(G4PionMinus::PionMinus() ) );
 
   meanField = new G4QMDMeanField();
   collision = new G4QMDCollision();
 
   excitationHandler = new G4ExcitationHandler;
   excitationHandler->SetDeexChannelsType( fCombined );
   //fEvaporation - 8 default channels
   //fCombined    - 8 default + 60 GEM
   //fGEM         - 2 default + 66 GEM
   evaporation = new G4Evaporation;
   excitationHandler->SetEvaporation( evaporation );
   setEvaporationCh();
 
   coulomb_collision_gamma_proj = 0.0;
   coulomb_collision_rx_proj = 0.0;
   coulomb_collision_rz_proj = 0.0;
   coulomb_collision_px_proj = 0.0;
   coulomb_collision_pz_proj = 0.0;
 
   coulomb_collision_gamma_targ = 0.0;
   coulomb_collision_rx_targ = 0.0;
   coulomb_collision_rz_targ = 0.0;
   coulomb_collision_px_targ = 0.0;
   coulomb_collision_pz_targ = 0.0;
}

◆ ~G4QMDReaction()

G4QMDReaction::~G4QMDReaction ( )

Definition at line 101 of file G4QMDReaction.cc.

{
   delete evaporation; 
   delete excitationHandler;
   delete collision;
   delete meanField;
}

Member Function Documentation

◆ ApplyYourself()

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

virtual

Reimplemented from G4HadronicInteraction.

Definition at line 110 of file G4QMDReaction.cc.

{
   //G4cout << "G4QMDReaction::ApplyYourself" << G4endl;
 
   theParticleChange.Clear();
 
   system = new G4QMDSystem;
 
   G4int proj_Z = 0;
   G4int proj_A = 0;
   const G4ParticleDefinition* proj_pd = ( const G4ParticleDefinition* ) projectile.GetDefinition();
   if ( proj_pd->GetParticleType() == "nucleus" )
   {
      proj_Z = proj_pd->GetAtomicNumber();
      proj_A = proj_pd->GetAtomicMass();
   }
   else
   {
      proj_Z = (int)( proj_pd->GetPDGCharge()/eplus );
      proj_A = 1;
   }
   //G4int targ_Z = int ( target.GetZ() + 0.5 );
   //G4int targ_A = int ( target.GetN() + 0.5 );
   //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
   G4int targ_Z = target.GetZ_asInt();
   G4int targ_A = target.GetA_asInt();
   const G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon( targ_Z , targ_A , 0.0 );
 
 
   //G4NistManager* nistMan = G4NistManager::Instance();
//   G4Element* G4NistManager::FindOrBuildElement( targ_Z );
 
   const G4DynamicParticle* proj_dp = new G4DynamicParticle ( proj_pd , projectile.Get4Momentum() );
   //const G4Element* targ_ele =  nistMan->FindOrBuildElement( targ_Z ); 
   //G4double aTemp = projectile.GetMaterial()->GetTemperature();
 
     //090331
  
   G4VCrossSectionDataSet* theXS = shenXS;
 
   G4double xs_0 = theXS->GetIsoCrossSection ( proj_dp , targ_Z , targ_A );
 
   // When the projectile is a pion
   if (proj_pd == G4PionPlus::PionPlus() ) {
     xs_0 = pipElNucXS->GetElementCrossSection(proj_dp, targ_Z, projectile.GetMaterial() ) +
            pipInelNucXS->GetElementCrossSection(proj_dp, targ_Z, projectile.GetMaterial() );      
   } else if (proj_pd == G4PionMinus::PionMinus() ) {
     xs_0 = pimElNucXS->GetElementCrossSection(proj_dp, targ_Z, projectile.GetMaterial() ) +
            pimInelNucXS->GetElementCrossSection(proj_dp, targ_Z, projectile.GetMaterial() );
   }
 
   //G4double xs_0 = genspaXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
   //G4double xs_0 = theXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
   //110822 
 
     G4double bmax_0 = std::sqrt( xs_0 / pi );
     //std::cout << "bmax_0 in fm (fermi) " <<  bmax_0/fermi << std::endl;
 
     //delete proj_dp; 
 
   G4bool elastic = true;
   
   std::vector< G4QMDNucleus* > nucleuses; // Secondary nuceluses 
   G4ThreeVector boostToReac; // ReactionSystem (CM or NN); 
   G4ThreeVector boostBackToLAB; // Reaction System to LAB; 
 
   G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ_pd->GetPDGMass()/GeV );
   G4ThreeVector boostLABtoCM = targ4p.findBoostToCM( proj_dp->Get4Momentum()/GeV ); // CM of target and proj; 
 
   G4double p1 = proj_dp->GetMomentum().mag()/GeV/proj_A; 
   G4double m1 = proj_dp->GetDefinition()->GetPDGMass()/GeV/proj_A;
   G4double e1 = std::sqrt( p1*p1 + m1*m1 ); 
   G4double e2 = targ_pd->GetPDGMass()/GeV/targ_A;
   G4double beta_nn = -p1 / ( e1+e2 );
 
   G4ThreeVector boostLABtoNN ( 0. , 0. , beta_nn ); // CM of NN; 
 
   G4double beta_nncm = ( - boostLABtoCM.beta() + boostLABtoNN.beta() ) / ( 1 - boostLABtoCM.beta() * boostLABtoNN.beta() ) ;  
 
   //std::cout << targ4p << std::endl; 
   //std::cout << proj_dp->Get4Momentum()<< std::endl; 
   //std::cout << beta_nncm << std::endl; 
   G4ThreeVector boostNNtoCM( 0. , 0. , beta_nncm ); // 
   G4ThreeVector boostCMtoNN( 0. , 0. , -beta_nncm ); // 
 
   boostToReac = boostLABtoNN; 
   boostBackToLAB = -boostLABtoNN; 
 
   delete proj_dp; 
 
   G4int icounter = 0;
   G4int icounter_max = 1024;
   while ( elastic ) // Loop checking, 11.03.2015, T. Koi
   {
      icounter++;
      if ( icounter > icounter_max ) { 
     G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
         break;
      }
 
// impact parameter 
      //G4double bmax = 1.05*(bmax_0/fermi);  // 10% for Peripheral reactions
      G4double bmax = envelopF*(bmax_0/fermi);
      G4double b = bmax * std::sqrt ( G4UniformRand() );
//071112
      //G4double b = 0;
      //G4double b = bmax;
      //G4double b = bmax/1.05 * 0.7 * G4UniformRand();
 
      //G4cout << "G4QMDRESULT bmax_0 = " << bmax_0/fermi << " fm, bmax = " << bmax << " fm , b = " << b  << " fm " << G4endl; 
 
      G4double plab = projectile.GetTotalMomentum()/GeV;
      G4double elab = ( projectile.GetKineticEnergy() + proj_pd->GetPDGMass() + targ_pd->GetPDGMass() )/GeV;
 
      calcOffSetOfCollision( b , proj_pd , targ_pd , plab , elab , bmax , boostCMtoNN );
 
// Projectile
      G4LorentzVector proj4pLAB = projectile.Get4Momentum()/GeV;
 
      G4QMDGroundStateNucleus* proj(NULL); 
      if ( projectile.GetDefinition()->GetParticleType() == "nucleus" 
        || projectile.GetDefinition()->GetParticleName() == "proton"
        || projectile.GetDefinition()->GetParticleName() == "neutron" )
      {
 
         proj_Z = proj_pd->GetAtomicNumber();
         proj_A = proj_pd->GetAtomicMass();
 
         proj = new G4QMDGroundStateNucleus( proj_Z , proj_A );
         //proj->ShowParticipants();
 
 
         meanField->SetSystem ( proj );
         proj->SetTotalPotential( meanField->GetTotalPotential() );
         proj->CalEnergyAndAngularMomentumInCM();
 
      }
 
// Target
      //G4int iz = int ( target.GetZ() );
      //G4int ia = int ( target.GetN() );
      //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
      G4int iz = int ( target.GetZ_asInt() );
      G4int ia = int ( target.GetA_asInt() );
 
      G4QMDGroundStateNucleus* targ = new G4QMDGroundStateNucleus( iz , ia );
 
      meanField->SetSystem (targ );
      targ->SetTotalPotential( meanField->GetTotalPotential() );
      targ->CalEnergyAndAngularMomentumInCM();
   
      //G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ->GetNuclearMass()/GeV );
// Boost Vector to CM
      //boostToCM = targ4p.findBoostToCM( proj4pLAB );
 
//    Target 
      for ( G4int i = 0 ; i < targ->GetTotalNumberOfParticipant() ; i++ )
      {
 
         G4ThreeVector p0 = targ->GetParticipant( i )->GetMomentum();
         G4ThreeVector r0 = targ->GetParticipant( i )->GetPosition();
 
         G4ThreeVector p ( p0.x() + coulomb_collision_px_targ 
                         , p0.y() 
                         , p0.z() * coulomb_collision_gamma_targ + coulomb_collision_pz_targ ); 
 
         G4ThreeVector r ( r0.x() + coulomb_collision_rx_targ 
                         , r0.y() 
                         , r0.z() / coulomb_collision_gamma_targ + coulomb_collision_rz_targ ); 
     
         system->SetParticipant( new G4QMDParticipant( targ->GetParticipant( i )->GetDefinition() , p , r ) );
         system->GetParticipant( i )->SetTarget();
 
      }
 
      G4LorentzVector proj4pCM = CLHEP::boostOf ( proj4pLAB , boostToReac );
      G4LorentzVector targ4pCM = CLHEP::boostOf ( targ4p , boostToReac );
 
//    Projectile
      //G4cout << "proj : " << proj << G4endl;
      //if ( proj != NULL )
      if ( proj_A != 1 )
      {
 
//    projectile is nucleus
 
         for ( G4int i = 0 ; i < proj->GetTotalNumberOfParticipant() ; i++ )
         {
 
            G4ThreeVector p0 = proj->GetParticipant( i )->GetMomentum();
            G4ThreeVector r0 = proj->GetParticipant( i )->GetPosition();
 
            G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 
                            , p0.y() 
                            , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 
 
            G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 
                            , r0.y() 
                            , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 
     
            system->SetParticipant( new G4QMDParticipant( proj->GetParticipant( i )->GetDefinition() , p  , r ) );
            system->GetParticipant ( i + targ->GetTotalNumberOfParticipant() )->SetProjectile();
         }
 
      }
      else
      {
 
//       projectile is particle
 
         // avoid multiple set in "elastic" loop
         //G4cout << "system Total Participants : " << system->GetTotalNumberOfParticipant() << ", target : " << targ->GetTotalNumberOfParticipant() << G4endl;
         if ( system->GetTotalNumberOfParticipant() == targ->GetTotalNumberOfParticipant() )
         {
 
            G4int i = targ->GetTotalNumberOfParticipant(); 
      
            G4ThreeVector p0( 0 ); 
            G4ThreeVector r0( 0 );
 
            G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 
                            , p0.y() 
                            , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 
 
            G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 
                            , r0.y() 
                            , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 
 
            system->SetParticipant( new G4QMDParticipant( (G4ParticleDefinition*)projectile.GetDefinition() , p , r ) );
            // This is not important becase only 1 projectile particle.
            system->GetParticipant ( i )->SetProjectile();
         }
 
      }
      //system->ShowParticipants();
 
      delete targ;
      delete proj;
 
   meanField->SetSystem ( system );
   collision->SetMeanField ( meanField );
 
// Time Evolution 
   //std::cout << "Start time evolution " << std::endl;
   //system->ShowParticipants();
   for ( G4int i = 0 ; i < maxTime ; i++ )
   {
      //G4cout << " do Paropagate " << i << " th time step. " << G4endl;
      meanField->DoPropagation( deltaT );
      //system->ShowParticipants();
      collision->CalKinematicsOfBinaryCollisions( deltaT );
 
      if ( i / 10 * 10 == i ) 
      {
         //G4cout << i << " th time step. " << G4endl;
         //system->ShowParticipants();
      } 
      //system->ShowParticipants();
   }
   //system->ShowParticipants();
 
 
   //std::cout << "Doing Cluster Judgment " << std::endl;
 
   nucleuses = meanField->DoClusterJudgment();
 
// Elastic Judgment  
 
   G4int numberOfSecondary = int ( nucleuses.size() ) + system->GetTotalNumberOfParticipant(); 
 
   G4int sec_a_Z = 0;
   G4int sec_a_A = 0;
   const G4ParticleDefinition* sec_a_pd = NULL;
   G4int sec_b_Z = 0;
   G4int sec_b_A = 0;
   const G4ParticleDefinition* sec_b_pd = NULL;
 
   if ( numberOfSecondary == 2 )
   {
 
      G4bool elasticLike_system = false;
      if ( nucleuses.size() == 2 ) 
      {
 
         sec_a_Z = nucleuses[0]->GetAtomicNumber();
         sec_a_A = nucleuses[0]->GetMassNumber();
         sec_b_Z = nucleuses[1]->GetAtomicNumber();
         sec_b_A = nucleuses[1]->GetMassNumber();
 
         if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_Z == targ_Z && sec_b_A == targ_A )
           || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_Z == proj_Z && sec_b_A == proj_A ) )
         {
            elasticLike_system = true;
         } 
 
      }
      else if ( nucleuses.size() == 1 ) 
      {
 
         sec_a_Z = nucleuses[0]->GetAtomicNumber();
         sec_a_A = nucleuses[0]->GetMassNumber();
         sec_b_pd = system->GetParticipant( 0 )->GetDefinition();
 
         if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_pd == targ_pd )
           || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_pd == proj_pd ) )
         {
            elasticLike_system = true;
         } 
 
      }  
      else
      {
 
         sec_a_pd = system->GetParticipant( 0 )->GetDefinition();
         sec_b_pd = system->GetParticipant( 1 )->GetDefinition();
 
         if ( ( sec_a_pd == proj_pd && sec_b_pd == targ_pd ) 
           || ( sec_a_pd == targ_pd && sec_b_pd == proj_pd ) ) 
         {
            elasticLike_system = true;
         } 
 
      } 
 
      if ( elasticLike_system == true )
      {
 
         G4bool elasticLike_energy = true;
//    Cal ExcitationEnergy 
         for ( G4int i = 0 ; i < int ( nucleuses.size() ) ; i++ )
         { 
 
            //meanField->SetSystem( nucleuses[i] );
            meanField->SetNucleus( nucleuses[i] );
            //nucleuses[i]->SetTotalPotential( meanField->GetTotalPotential() );
            //nucleuses[i]->CalEnergyAndAngularMomentumInCM();
 
            if ( nucleuses[i]->GetExcitationEnergy()*GeV > 1.0*MeV ) elasticLike_energy = false;  
 
         } 
 
//    Check Collision 
         G4bool withCollision = true;
         if ( system->GetNOCollision() == 0 ) withCollision = false;
 
//    Final judegement for Inelasitc or Elastic;
//
//       ElasticLike without Collision 
         //if ( elasticLike_energy == true && withCollision == false ) elastic = true;  // ielst = 0
//       ElasticLike with Collision 
         //if ( elasticLike_energy == true && withCollision == true ) elastic = true;   // ielst = 1 
//       InelasticLike without Collision 
         //if ( elasticLike_energy == false ) elastic = false;                          // ielst = 2                
         if ( frag == true )
            if ( elasticLike_energy == false ) elastic = false;
//       InelasticLike with Collision 
         if ( elasticLike_energy == false && withCollision == true ) elastic = false; // ielst = 3
 
      }
 
      }
      else
      {
 
//       numberOfSecondary != 2 
         elastic = false;
 
      }
 
//071115
      //G4cout << elastic << G4endl;
      // if elastic is true try again from sampling of impact parameter 
 
      if ( elastic == true )
      {
         // delete this nucleues
         for ( std::vector< G4QMDNucleus* >::iterator
               it = nucleuses.begin() ; it != nucleuses.end() ; it++ )
         {
            delete *it;
         }
         nucleuses.clear();
      }
 
   } 
 
 
// Statical Decay Phase
 
   for ( std::vector< G4QMDNucleus* >::iterator it
       = nucleuses.begin() ; it != nucleuses.end() ; it++ )
   {
 
/*
      G4cout << "G4QMDRESULT "
             << (*it)->GetAtomicNumber() 
             << " " 
             << (*it)->GetMassNumber() 
             << " " 
             << (*it)->Get4Momentum() 
             << " " 
             << (*it)->Get4Momentum().vect() 
             << " " 
             << (*it)->Get4Momentum().restMass() 
             << " " 
             << (*it)->GetNuclearMass()/GeV 
             << G4endl;
*/
 
      meanField->SetNucleus ( *it );
 
      if ( (*it)->GetAtomicNumber() == 0  // neutron cluster
        || (*it)->GetAtomicNumber() == (*it)->GetMassNumber() ) // proton cluster
      {
         // push back system 
         for ( G4int i = 0 ; i < (*it)->GetTotalNumberOfParticipant() ; i++ )
         {
            G4QMDParticipant* aP = new G4QMDParticipant( ( (*it)->GetParticipant( i ) )->GetDefinition() , ( (*it)->GetParticipant( i ) )->GetMomentum() , ( (*it)->GetParticipant( i ) )->GetPosition() );  
            system->SetParticipant ( aP );  
         } 
         continue;  
      }
 
      G4double nucleus_e = std::sqrt ( G4Pow::GetInstance()->powN ( (*it)->GetNuclearMass()/GeV , 2 ) + G4Pow::GetInstance()->powN ( (*it)->Get4Momentum().vect().mag() , 2 ) );
      G4LorentzVector nucleus_p4CM ( (*it)->Get4Momentum().vect() , nucleus_e ); 
 
//      std::cout << "G4QMDRESULT nucleus deltaQ " << deltaQ << std::endl;
 
      G4int ia = (*it)->GetMassNumber();
      G4int iz = (*it)->GetAtomicNumber();
 
      G4LorentzVector lv ( G4ThreeVector( 0.0 ) , (*it)->GetExcitationEnergy()*GeV + G4IonTable::GetIonTable()->GetIonMass( iz , ia ) );
 
      G4Fragment* aFragment = new G4Fragment( ia , iz , lv );
 
      G4ReactionProductVector* rv;
      rv = excitationHandler->BreakItUp( *aFragment );
      G4bool notBreak = true;
      for ( G4ReactionProductVector::iterator itt
          = rv->begin() ; itt != rv->end() ; itt++ )
      {
 
          notBreak = false;
          // Secondary from this nucleus (*it) 
          const G4ParticleDefinition* pd = (*itt)->GetDefinition();
 
          G4LorentzVector p4 ( (*itt)->GetMomentum()/GeV , (*itt)->GetTotalEnergy()/GeV );  //in nucleus(*it) rest system
          G4LorentzVector p4_CM = CLHEP::boostOf( p4 , -nucleus_p4CM.findBoostToCM() );  // Back to CM
          G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB  
 
 
//090122
          //theParticleChange.AddSecondary( dp ); 
          if ( !( pd->GetAtomicNumber() == 4 && pd->GetAtomicMass() == 8 ) )
          {
             //G4cout << "pd out of notBreak loop : " << pd->GetParticleName() << G4endl;
             G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
             theParticleChange.AddSecondary( dp ); 
          }
          else
          {
             //Be8 -> Alpha + Alpha + Q
             G4ThreeVector randomized_direction( G4UniformRand() , G4UniformRand() , G4UniformRand() );
             randomized_direction = randomized_direction.unit();
             G4double q_decay = (*itt)->GetMass() - 2*G4Alpha::Alpha()->GetPDGMass();
             G4double p_decay = std::sqrt ( G4Pow::GetInstance()->powN(G4Alpha::Alpha()->GetPDGMass()+q_decay/2,2) - G4Pow::GetInstance()->powN(G4Alpha::Alpha()->GetPDGMass() , 2 ) ); 
             G4LorentzVector p4_a1 ( p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 );  //in Be8 rest system
             
             G4LorentzVector p4_a1_Be8 = CLHEP::boostOf ( p4_a1/GeV , -p4.findBoostToCM() );
             G4LorentzVector p4_a1_CM = CLHEP::boostOf ( p4_a1_Be8 , -nucleus_p4CM.findBoostToCM() );
             G4LorentzVector p4_a1_LAB = CLHEP::boostOf ( p4_a1_CM , boostBackToLAB );
 
             G4LorentzVector p4_a2 ( -p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 );  //in Be8 rest system
             
             G4LorentzVector p4_a2_Be8 = CLHEP::boostOf ( p4_a2/GeV , -p4.findBoostToCM() );
             G4LorentzVector p4_a2_CM = CLHEP::boostOf ( p4_a2_Be8 , -nucleus_p4CM.findBoostToCM() );
             G4LorentzVector p4_a2_LAB = CLHEP::boostOf ( p4_a2_CM , boostBackToLAB );
             
             G4DynamicParticle* dp1 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a1_LAB*GeV );  
             G4DynamicParticle* dp2 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a2_LAB*GeV );  
             theParticleChange.AddSecondary( dp1 ); 
             theParticleChange.AddSecondary( dp2 ); 
          }
//090122
 
/*
          G4cout
                << "Regist Secondary "
                << (*itt)->GetDefinition()->GetParticleName()
                << " "
                << (*itt)->GetMomentum()/GeV
                << " "
                << (*itt)->GetKineticEnergy()/GeV
                << " "
                << (*itt)->GetMass()/GeV
                << " "
                << (*itt)->GetTotalEnergy()/GeV
                << " "
                << (*itt)->GetTotalEnergy()/GeV * (*itt)->GetTotalEnergy()/GeV
                 - (*itt)->GetMomentum()/GeV * (*itt)->GetMomentum()/GeV
                << " "
                << nucleus_p4CM.findBoostToCM() 
                << " "
                << p4
                << " "
                << p4_CM
                << " "
                << p4_LAB
                << G4endl;
*/
 
      }
      if ( notBreak == true )
      {
 
         const G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon( (*it)->GetAtomicNumber() , (*it)->GetMassNumber(), (*it)->GetExcitationEnergy()*GeV );
             //G4cout << "pd in notBreak loop : " << pd->GetParticleName() << G4endl;
         G4LorentzVector p4_CM = nucleus_p4CM;
         G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB  
         G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
         theParticleChange.AddSecondary( dp ); 
 
      }
 
      for ( G4ReactionProductVector::iterator itt
          = rv->begin() ; itt != rv->end() ; itt++ )
      {
          delete *itt;
      }
      delete rv;
 
      delete aFragment;
 
   }
 
 
 
   for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i++ )
   {
      // Secondary particles 
 
      const G4ParticleDefinition* pd = system->GetParticipant( i )->GetDefinition();
      G4LorentzVector p4_CM = system->GetParticipant( i )->Get4Momentum();
      G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB );
      G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
      theParticleChange.AddSecondary( dp ); 
      //G4cout << "In the last theParticleChange loop : " << pd->GetParticleName() << G4endl;
 
/*
      G4cout << "G4QMDRESULT "
      << "r" << i << " " << system->GetParticipant ( i ) -> GetPosition() << " "
      << "p" << i << " " << system->GetParticipant ( i ) -> Get4Momentum()
      << G4endl;
*/
 
   }
 
   for ( std::vector< G4QMDNucleus* >::iterator it
       = nucleuses.begin() ; it != nucleuses.end() ; it++ )
   {
      delete *it;  // delete nulceuse 
   }
   nucleuses.clear();
 
   system->Clear();
   delete system; 
 
   theParticleChange.SetStatusChange( stopAndKill );
 
   for (G4int i = 0; i < G4int(theParticleChange.GetNumberOfSecondaries() ); i++)
   {
     //G4cout << "Particle : " << theParticleChange.GetSecondary(i)->GetParticle()->GetParticleDefinition()->GetParticleName() << G4endl;
     //G4cout << "KEnergy : " << theParticleChange.GetSecondary(i)->GetParticle()->GetKineticEnergy() << G4endl;
     //G4cout << "KEnergy : " << theParticleChange.GetSecondary(i)->GetCreatorModelType() << G4endl;
     theParticleChange.GetSecondary(i)->SetCreatorModelType(1111);
   }
 
   return &theParticleChange;
 
}

◆ GetExcitationHandler()

G4ExcitationHandler * G4QMDReaction::GetExcitationHandler ( )

inline

Definition at line 72 of file G4QMDReaction.hh.

72{return excitationHandler;};

◆ GetFinalStates()

std::vector< G4QMDSystem * > G4QMDReaction::GetFinalStates ( )

◆ ModelDescription()

void G4QMDReaction::ModelDescription ( std::ostream & outFile ) const

virtual

Reimplemented from G4HadronicInteraction.

Definition at line 849 of file G4QMDReaction.cc.

{
   outFile << "Lorentz covarianted Quantum Molecular Dynamics model for nucleus (particle) vs nucleus reactions\n";
}

◆ SetDT()

void G4QMDReaction::SetDT ( G4double t )

inline

Definition at line 78 of file G4QMDReaction.hh.

78{ deltaT = t; };

◆ SetEF()

void G4QMDReaction::SetEF ( G4double x )

inline

Definition at line 79 of file G4QMDReaction.hh.

79{ envelopF = x; };

◆ SetTMAX()

void G4QMDReaction::SetTMAX ( G4int i )

inline

Definition at line 77 of file G4QMDReaction.hh.

77{ maxTime = i; };

◆ UnUseGEM()

void G4QMDReaction::UnUseGEM ( )

inline

Definition at line 74 of file G4QMDReaction.hh.

74{gem = false; setEvaporationCh();};

◆ UseFRAG()

void G4QMDReaction::UseFRAG ( )

inline

Definition at line 75 of file G4QMDReaction.hh.

75{frag = true;};

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

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ G4QMDReaction()

◆ ~G4QMDReaction()

Member Function Documentation

◆ ApplyYourself()

◆ GetExcitationHandler()

◆ GetFinalStates()

◆ ModelDescription()

◆ SetDT()

◆ SetEF()

◆ SetTMAX()

◆ UnUseGEM()

◆ UseFRAG()