G4LowEIonFragmentation Class Reference

#include <G4LowEIonFragmentation.hh>

Inheritance diagram for G4LowEIonFragmentation:

G4HadronicInteraction

Public Member Functions

 G4LowEIonFragmentation (G4ExcitationHandler *const value)
 G4LowEIonFragmentation ()
virtual ~G4LowEIonFragmentation ()
virtual G4HadFinalStateApplyYourself (const G4HadProjectile &thePrimary, G4Nucleus &theNucleus)

Static Public Member Functions

static G4double GetCrossSection ()

Detailed Description

Definition at line 50 of file G4LowEIonFragmentation.hh.


Constructor & Destructor Documentation

G4LowEIonFragmentation::G4LowEIonFragmentation ( G4ExcitationHandler *const   value  ) 

Definition at line 53 of file G4LowEIonFragmentation.cc.

References G4Proton::Proton().

00054 {
00055   theHandler = value;
00056   theModel = new G4PreCompoundModel(theHandler);
00057   proton = G4Proton::Proton();
00058 }

G4LowEIonFragmentation::G4LowEIonFragmentation (  ) 

Definition at line 60 of file G4LowEIonFragmentation.cc.

References G4Proton::Proton().

00061 {
00062   theHandler = new G4ExcitationHandler;
00063   theModel = new G4PreCompoundModel(theHandler);
00064   proton = G4Proton::Proton();
00065 }

G4LowEIonFragmentation::~G4LowEIonFragmentation (  )  [virtual]

Definition at line 67 of file G4LowEIonFragmentation.cc.

00068 {
00069   delete theModel;
00070 }


Member Function Documentation

G4HadFinalState * G4LowEIonFragmentation::ApplyYourself ( const G4HadProjectile thePrimary,
G4Nucleus theNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 73 of file G4LowEIonFragmentation.cc.

References G4HadFinalState::AddSecondary(), G4ExcitationHandler::BreakItUp(), G4HadFinalState::Clear(), G4PreCompoundModel::DeExcite(), G4lrint(), G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4HadProjectile::GetGlobalTime(), G4Fancy3DNucleus::GetNextNucleon(), G4NucleiProperties::GetNuclearMass(), G4Fancy3DNucleus::GetOuterRadius(), G4Nucleon::GetParticleType(), G4ParticleDefinition::GetPDGCharge(), G4Nucleon::GetPosition(), G4HadProjectile::GetTotalEnergy(), G4Nucleus::GetZ_asInt(), G4Fancy3DNucleus::Init(), G4INCL::Math::pi, G4Fragment::SetCreationTime(), G4HadFinalState::SetEnergyChange(), G4Fragment::SetNumberOfExcitedParticle(), G4Fragment::SetNumberOfHoles(), G4HadFinalState::SetStatusChange(), G4Fancy3DNucleus::StartLoop(), and stopAndKill.

00074 {
00075   area = 0;
00076   // initialize the particle change
00077   theResult.Clear();
00078   theResult.SetStatusChange( stopAndKill );
00079   theResult.SetEnergyChange( 0.0 );
00080 
00081   // Get Target A, Z
00082   G4int aTargetA = theNucleus.GetA_asInt();
00083   G4int aTargetZ = theNucleus.GetZ_asInt();
00084 
00085   // Get Projectile A, Z
00086   G4int aProjectileA = thePrimary.GetDefinition()->GetBaryonNumber();
00087   G4int aProjectileZ = G4lrint(thePrimary.GetDefinition()->GetPDGCharge()/eplus);
00088 
00089   // Get Maximum radius of both
00090   
00091   G4Fancy3DNucleus aPrim;
00092   aPrim.Init(aProjectileA, aProjectileZ);
00093   G4double projectileOuterRadius = aPrim.GetOuterRadius();
00094   
00095   G4Fancy3DNucleus aTarg;
00096   aTarg.Init(aTargetA, aTargetZ);
00097   G4double targetOuterRadius = aTarg.GetOuterRadius();
00098 
00099   // Get the Impact parameter
00100   G4int particlesFromProjectile = 0;
00101   G4int chargedFromProjectile = 0;
00102   G4double impactParameter = 0;
00103   G4double x,y;
00104   G4Nucleon * pNucleon;
00105   // need at lease one particle from the projectile model beyond the 
00106   // projectileHorizon.
00107   while(0==particlesFromProjectile)
00108   {
00109     do
00110     {
00111       x = 2*G4UniformRand() - 1;
00112       y = 2*G4UniformRand() - 1;
00113     }
00114     while(x*x + y*y > 1);
00115     impactParameter = std::sqrt(x*x+y*y)*(targetOuterRadius+projectileOuterRadius);
00116     ++totalTries;
00117     area = pi*(targetOuterRadius+projectileOuterRadius)*
00118               (targetOuterRadius+projectileOuterRadius);
00119     G4double projectileHorizon = impactParameter-targetOuterRadius; 
00120     
00121     // Empirical boundary transparency.
00122     G4double empirical = G4UniformRand();
00123     if(projectileHorizon > empirical*projectileOuterRadius) { continue; }
00124     
00125     // Calculate the number of nucleons involved in collision
00126     // From projectile
00127     aPrim.StartLoop();
00128     while((pNucleon = aPrim.GetNextNucleon()))
00129     {
00130       if(pNucleon->GetPosition().y()>projectileHorizon)
00131       {
00132         // We have one
00133         ++particlesFromProjectile;
00134         if(pNucleon->GetParticleType() == proton) 
00135         {
00136           ++chargedFromProjectile;
00137         } 
00138       }
00139     }
00140   }
00141   ++hits;
00142 
00143   // From target:
00144   G4double targetHorizon = impactParameter-projectileOuterRadius;
00145   G4int chargedFromTarget = 0;
00146   G4int particlesFromTarget = 0;
00147   aTarg.StartLoop();  
00148   while((pNucleon = aTarg.GetNextNucleon()))
00149   {
00150     if(pNucleon->GetPosition().y()>targetHorizon)
00151     {
00152       // We have one
00153       ++particlesFromTarget;
00154       if(pNucleon->GetParticleType() == proton) 
00155       {
00156         ++chargedFromTarget;
00157       }
00158     }
00159   }
00160   
00161   // Energy sharing between projectile and target. 
00162   // Note that this is a quite simplistic kinetically.
00163   G4ThreeVector momentum = thePrimary.Get4Momentum().vect();
00164   G4double w = (G4double)particlesFromProjectile/(G4double)aProjectileA;
00165   
00166   G4double projTotEnergy = thePrimary.GetTotalEnergy();  
00167   G4double targetMass = G4NucleiProperties::GetNuclearMass(aTargetA, aTargetZ);
00168   G4LorentzVector fragment4Momentum(momentum*w, projTotEnergy*w + targetMass);
00169  
00170   // take the nucleons and fill the Fragments
00171   G4Fragment anInitialState(aTargetA+particlesFromProjectile,
00172                             aTargetZ+chargedFromProjectile,
00173                             fragment4Momentum);
00174   // M.A. Cortes fix
00175   //anInitialState.SetNumberOfParticles(particlesFromProjectile);
00176   anInitialState.SetNumberOfExcitedParticle(particlesFromProjectile+particlesFromTarget,
00177                                             chargedFromProjectile + chargedFromTarget);
00178   anInitialState.SetNumberOfHoles(particlesFromProjectile+particlesFromTarget,
00179                                   chargedFromProjectile + chargedFromTarget);
00180   G4double time = thePrimary.GetGlobalTime();
00181   anInitialState.SetCreationTime(time);
00182 
00183   // Fragment the Fragment using Pre-compound
00184   G4ReactionProductVector* thePreCompoundResult = theModel->DeExcite(anInitialState);
00185   
00186   // De-excite the projectile using ExcitationHandler
00187   G4ReactionProductVector * theExcitationResult = 0; 
00188   if(particlesFromProjectile < aProjectileA)
00189   {
00190     G4LorentzVector residual4Momentum(momentum*(1.0-w), projTotEnergy*(1.0-w));  
00191  
00192     G4Fragment initialState2(aProjectileA-particlesFromProjectile,
00193                              aProjectileZ-chargedFromProjectile,
00194                              residual4Momentum );
00195 
00196     // half of particles are excited (?!)
00197     G4int pinit = (aProjectileA-particlesFromProjectile)/2;
00198     G4int cinit = (aProjectileZ-chargedFromProjectile)/2;
00199 
00200     initialState2.SetNumberOfExcitedParticle(pinit,cinit);
00201     initialState2.SetNumberOfHoles(pinit,cinit);
00202     initialState2.SetCreationTime(time);
00203 
00204     theExcitationResult = theHandler->BreakItUp(initialState2);
00205   }
00206 
00207   // Fill the particle change and clear intermediate vectors
00208   G4int nexc = 0;
00209   G4int npre = 0;
00210   if(theExcitationResult)  { nexc = theExcitationResult->size(); }
00211   if(thePreCompoundResult) { npre = thePreCompoundResult->size();}
00212   
00213   if(nexc > 0) {
00214     for(G4int k=0; k<nexc; ++k) {
00215       G4ReactionProduct* p = (*theExcitationResult)[k];
00216       theResult.AddSecondary(new G4DynamicParticle(p->GetDefinition(),p->GetMomentum()));
00217       delete p;
00218     }
00219   }
00220   
00221   if(npre > 0) {
00222     for(G4int k=0; k<npre; ++k) {
00223       G4ReactionProduct* p = (*thePreCompoundResult)[k];
00224       theResult.AddSecondary(new G4DynamicParticle(p->GetDefinition(),p->GetMomentum()));
00225       delete p;
00226     }
00227   }
00228   
00229   delete thePreCompoundResult;
00230   delete theExcitationResult;
00231 
00232   // return the particle change
00233   return &theResult;
00234   
00235 }

static G4double G4LowEIonFragmentation::GetCrossSection (  )  [inline, static]

Definition at line 63 of file G4LowEIonFragmentation.hh.

00064   {
00065     //    clog << "area/millibarn = "<<area/millibarn<<G4endl;
00066     //    clog << "hits = "<<hits<<G4endl;
00067     //    clog << "totalTries = "<<totalTries<<G4endl;
00068     return area*hits/(static_cast<G4double>(totalTries)*CLHEP::millibarn);
00069   }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:27 2013 for Geant4 by  doxygen 1.4.7