G4LightIonQMDReaction Class Reference

#include <G4LightIonQMDReaction.hh>

Inheritance diagram for G4LightIonQMDReaction:

Public Member Functions
	G4LightIonQMDReaction ()
	~G4LightIonQMDReaction () override
std::vector< G4QMDSystem * >	GetFinalStates ()
G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) override
G4ExcitationHandler *	GetExcitationHandler ()
void	UnUseGEM ()
void	UseFRAG ()
void	SetTMAX (G4int i)
void	SetDT (G4double t)
void	SetEF (G4double x)
void	ModelDescription (std::ostream &outFile) const override
	G4LightIonQMDReaction (const G4LightIonQMDReaction &right)=delete
const G4LightIonQMDReaction &	operator= (const G4LightIonQMDReaction &right)=delete
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 ()
	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 59 of file G4LightIonQMDReaction.hh.

Constructor & Destructor Documentation

G4LightIonQMDReaction() [1/2]

G4LightIonQMDReaction::G4LightIonQMDReaction

(

)

Definition at line 70 of file G4LightIonQMDReaction.cc.

: G4HadronicInteraction("LightIonQMDModel")
, 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 )
, secID( -1 )
{
   G4cout << "G4LightIonQMDReaction::G4LightIonQMDReaction" << G4endl;
   G4cout << "Recommended Energy of LightIonQMD: 30MeV/u - 500MeV/u" << G4endl;
    
   theXS = new G4CrossSectionInelastic( new G4ComponentGGNuclNuclXsc );
   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 G4LightIonQMDMeanField();
   collision = new G4LightIonQMDCollision();
 
   excitationHandler = new G4ExcitationHandler();
   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;
 
   secID = G4PhysicsModelCatalog::GetModelID( "model_LightIonQMDModel" );
}

Referenced by G4LightIonQMDReaction(), and operator=().

~G4LightIonQMDReaction()

G4LightIonQMDReaction::~G4LightIonQMDReaction ( )

override

Definition at line 118 of file G4LightIonQMDReaction.cc.

{
   delete excitationHandler;
   delete collision;
   delete meanField;
}

G4LightIonQMDReaction() [2/2]

G4LightIonQMDReaction::G4LightIonQMDReaction ( const G4LightIonQMDReaction & right )

delete

Member Function Documentation

ApplyYourself()

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

overridevirtual

Reimplemented from G4HadronicInteraction.

Definition at line 126 of file G4LightIonQMDReaction.cc.

{
   //G4cout << "G4LightIonQMDReaction::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();
 
   // Glauber-Gribov nucleus-nucleus cross section does not have GetIsoCrossSection,
   // therefore call GetElementCrossSection instead.
   //G4double xs_0 = theXS->GetIsoCrossSection ( proj_dp , targ_Z , targ_A );
   G4double xs_0 = theXS->GetElementCrossSection( proj_dp , targ_Z , projectile.GetMaterial() );
 
   // 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< G4LightIonQMDNucleus* > 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;
 
      G4LightIonQMDGroundStateNucleus* 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 G4LightIonQMDGroundStateNucleus( proj_Z , proj_A );
          //proj->ShowParticipants();
          
          
          meanField->SetSystem ( proj );
          if ( proj_A != 1 )
          {
              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() );
      G4LightIonQMDGroundStateNucleus* targ = new G4LightIonQMDGroundStateNucleus( iz , ia );
 
      meanField->SetSystem (targ );
      if ( ia != 1 )
      {
          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;
         }
          // QMD should be inelastic collision, so that nucleon-nucleon collision should also be inelastic in this phase. by Y-H. S and A. H, Mar. 6, 2023.
          if ( (proj_pd->GetParticleName() == "proton" && targ_pd->GetParticleName()  == "proton")
              || (proj_pd->GetParticleName() == "neutron" && targ_pd->GetParticleName()  == "proton")
              || (proj_pd->GetParticleName() == "pi+" && targ_pd->GetParticleName()  == "proton")
              || (proj_pd->GetParticleName() == "pi-" && targ_pd->GetParticleName()  == "proton"))
          {
              elasticLike_system = false;
              //G4cout << "elasticLike_system = false proton NOCollision " << system->GetNOCollision() << G4endl;
              if ( system->GetNOCollision() == 1 || icounter+900 > icounter_max) elastic = false;
          }
          // Addition -- end
      } 
 
      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< G4LightIonQMDNucleus* >::iterator
               it = nucleuses.begin() ; it != nucleuses.end() ; it++ )
         {
            delete *it;
         }
         nucleuses.clear();
         // system->Clear() should be included here. Otherwise, the nucleon is repeatedly regstered if the nucleon is the projectile. by Y-H. S. and A. H, Mar. 6, 2023.
         system->Clear();
      }
 
   } 
 
// random rotation
   G4double randAzimuthAngle = CLHEP::twopi*G4UniformRand();
   G4ThreeVector zAxis(0.,0.,1.);
 
// Statical Decay Phase
 
   for ( std::vector< G4LightIonQMDNucleus* >::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  
          p4_LAB.rotate(randAzimuthAngle, zAxis); // random rotation
 
//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 );
             p4_a1_LAB.rotate(randAzimuthAngle, zAxis); // random rotation
 
             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 );
             p4_a2_LAB.rotate(randAzimuthAngle, zAxis); // random rotation
             
             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  
         p4_LAB.rotate(randAzimuthAngle, zAxis); // random rotation
         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 );
      p4_LAB.rotate(randAzimuthAngle, zAxis); // random rotation
      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< G4LightIonQMDNucleus* >::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 << "modelID : " << theParticleChange.GetSecondary(i)->GetCreatorModelID() << G4endl;
     theParticleChange.GetSecondary(i)->SetCreatorModelID(secID);
   }
 
   return &theParticleChange;
 
}

GetExcitationHandler()

G4ExcitationHandler * G4LightIonQMDReaction::GetExcitationHandler ( )

inline

Definition at line 70 of file G4LightIonQMDReaction.hh.

{return excitationHandler;};

GetFinalStates()

std::vector< G4QMDSystem * > G4LightIonQMDReaction::GetFinalStates

(

)

ModelDescription()

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

overridevirtual

Reimplemented from G4HadronicInteraction.

Definition at line 882 of file G4LightIonQMDReaction.cc.

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

operator=()

const G4LightIonQMDReaction & G4LightIonQMDReaction::operator= ( const G4LightIonQMDReaction & right )

delete

SetDT()

void G4LightIonQMDReaction::SetDT ( G4double t )

inline

Definition at line 76 of file G4LightIonQMDReaction.hh.

{ deltaT = t; };

SetEF()

void G4LightIonQMDReaction::SetEF ( G4double x )

inline

Definition at line 77 of file G4LightIonQMDReaction.hh.

{ envelopF = x; };

SetTMAX()

void G4LightIonQMDReaction::SetTMAX ( G4int i )

inline

Definition at line 75 of file G4LightIonQMDReaction.hh.

{ maxTime = i; };

UnUseGEM()

void G4LightIonQMDReaction::UnUseGEM ( )

inline

Definition at line 72 of file G4LightIonQMDReaction.hh.

{gem = false; setEvaporationCh();};

UseFRAG()

void G4LightIonQMDReaction::UseFRAG ( )

inline

Definition at line 73 of file G4LightIonQMDReaction.hh.

{frag = true;};

---

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

• G4LightIonQMDReaction.hh
• G4LightIonQMDReaction.cc