G4HadronElastic Class Reference

#include <G4HadronElastic.hh>

Inheritance diagram for G4HadronElastic:

G4HadronicInteraction G4AntiNuclElastic G4CHIPSElastic G4ChipsElasticModel G4DiffuseElastic G4ElasticHadrNucleusHE G4NuclNuclDiffuseElastic

Public Member Functions

 G4HadronElastic (const G4String &name="hElasticLHEP")
virtual ~G4HadronElastic ()
virtual G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
void SetLowestEnergyLimit (G4double value)
G4double LowestEnergyLimit () const
G4double ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
virtual void Description () const

Detailed Description

Definition at line 50 of file G4HadronElastic.hh.


Constructor & Destructor Documentation

G4HadronElastic::G4HadronElastic ( const G4String name = "hElasticLHEP"  ) 

Definition at line 45 of file G4HadronElastic.cc.

References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4Neutron::Neutron(), G4Proton::Proton(), G4HadronicInteraction::SetMaxEnergy(), and G4HadronicInteraction::SetMinEnergy().

00046   : G4HadronicInteraction(name)
00047 {
00048   SetMinEnergy( 0.0*GeV );
00049   SetMaxEnergy( 100.*TeV );
00050   lowestEnergyLimit= 1.e-6*eV;  
00051 
00052   theProton   = G4Proton::Proton();
00053   theNeutron  = G4Neutron::Neutron();
00054   theDeuteron = G4Deuteron::Deuteron();
00055   theAlpha    = G4Alpha::Alpha();
00056   //Description();
00057 }

G4HadronElastic::~G4HadronElastic (  )  [virtual]

Definition at line 60 of file G4HadronElastic.cc.

00061 {}


Member Function Documentation

G4HadFinalState * G4HadronElastic::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus targetNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 93 of file G4HadronElastic.cc.

References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), G4cout, G4endl, G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetDefinition(), G4IonTable::GetIon(), G4ParticleTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetParticleName(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4HadronicInteraction::GetRecoilEnergyThreshold(), G4HadProjectile::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4He3::He3(), SampleInvariantT(), G4HadFinalState::SetEnergyChange(), G4HadFinalState::SetLocalEnergyDeposit(), G4HadFinalState::SetMomentumChange(), G4HadronicInteraction::theParticleChange, G4Triton::Triton(), and G4HadronicInteraction::verboseLevel.

00095 {
00096   theParticleChange.Clear();
00097 
00098   const G4HadProjectile* aParticle = &aTrack;
00099   G4double ekin = aParticle->GetKineticEnergy();
00100   if(ekin <= lowestEnergyLimit) {
00101     theParticleChange.SetEnergyChange(ekin);
00102     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00103     return &theParticleChange;
00104   }
00105 
00106   G4int A = targetNucleus.GetA_asInt();
00107   G4int Z = targetNucleus.GetZ_asInt();
00108 
00109   G4double plab = aParticle->GetTotalMomentum();
00110 
00111   // Scattered particle referred to axis of incident particle
00112   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00113   G4double m1 = theParticle->GetPDGMass();
00114 
00115   if (verboseLevel>1) {
00116     G4cout << "G4HadronElastic: " 
00117            << aParticle->GetDefinition()->GetParticleName() 
00118            << " Plab(GeV/c)= " << plab/GeV  
00119            << " Ekin(MeV) = " << ekin/MeV 
00120            << " scattered off Z= " << Z 
00121            << " A= " << A 
00122            << G4endl;
00123   }
00124 
00125   G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
00126   G4LorentzVector lv1 = aParticle->Get4Momentum();
00127   G4LorentzVector lv(0.0,0.0,0.0,mass2);   
00128   lv += lv1;
00129 
00130   G4ThreeVector bst = lv.boostVector();
00131   lv1.boost(-bst);
00132 
00133   G4ThreeVector p1 = lv1.vect();
00134   G4double momentumCMS = p1.mag();
00135   G4double tmax = 4.0*momentumCMS*momentumCMS;
00136 
00137   // Sampling in CM system
00138   G4double t    = SampleInvariantT(theParticle, plab, Z, A);
00139   G4double phi  = G4UniformRand()*CLHEP::twopi;
00140   G4double cost = 1. - 2.0*t/tmax;
00141   G4double sint;
00142 
00143   // problem in sampling
00144   if(cost > 1.0 || cost < -1.0) {
00145     //if(verboseLevel > 0) {
00146       G4cout << "G4HadronElastic WARNING (1 - cost)= " << 1 - cost
00147              << " after scattering of " 
00148              << aParticle->GetDefinition()->GetParticleName()
00149              << " p(GeV/c)= " << plab/GeV
00150              << " on an ion Z= " << Z << " A= " << A
00151              << G4endl;
00152       //}
00153     cost = 1.0;
00154     sint = 0.0;
00155 
00156     // normal situation
00157   } else  {
00158     sint = std::sqrt((1.0-cost)*(1.0+cost));
00159   }    
00160   if (verboseLevel>1) {
00161     G4cout << " t= " << t << " tmax(GeV^2)= " << tmax/(GeV*GeV) 
00162            << " Pcms(GeV)= " << momentumCMS/GeV << " cos(t)=" << cost 
00163            << " sin(t)=" << sint << G4endl;
00164   }
00165   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
00166   v1 *= momentumCMS;
00167   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),
00168                        std::sqrt(momentumCMS*momentumCMS + m1*m1));
00169 
00170   nlv1.boost(bst); 
00171 
00172   G4double eFinal = nlv1.e() - m1;
00173   if (verboseLevel > 1) {
00174     G4cout <<" m= " << m1 << " Efin(MeV)= " << eFinal 
00175            << " Proj: 4-mom " << lv1 << " Final: " << nlv1 
00176            << G4endl;
00177   }
00178   if(eFinal <= lowestEnergyLimit) {
00179     if(eFinal < 0.0 && verboseLevel > 0) {
00180       G4cout << "G4HadronElastic WARNING Efinal= " << eFinal
00181              << " after scattering of " 
00182              << aParticle->GetDefinition()->GetParticleName()
00183              << " p(GeV/c)= " << plab/GeV
00184              << " on an ion Z= " << Z << " A= " << A
00185              << G4endl;
00186     }
00187     theParticleChange.SetEnergyChange(0.0);
00188     nlv1 = G4LorentzVector(0.0,0.0,0.0,m1);
00189 
00190   } else {
00191     theParticleChange.SetMomentumChange(nlv1.vect().unit());
00192     theParticleChange.SetEnergyChange(eFinal);
00193   }  
00194 
00195   lv -= nlv1;
00196   G4double erec =  lv.e() - mass2;
00197   if (verboseLevel > 1) {
00198     G4cout << "Recoil: " <<" m= " << mass2 << " Erec(MeV)= " << erec
00199            << " 4-mom: " << lv 
00200            << G4endl;
00201   }
00202  
00203   if(erec > GetRecoilEnergyThreshold()) {
00204     G4ParticleDefinition * theDef = 0;
00205     if(Z == 1 && A == 1)       { theDef = theProton; }
00206     else if (Z == 1 && A == 2) { theDef = theDeuteron; }
00207     else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
00208     else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
00209     else if (Z == 2 && A == 4) { theDef = theAlpha; }
00210     else {
00211       theDef = 
00212         G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
00213     }
00214     G4DynamicParticle * aSec = new G4DynamicParticle(theDef, lv);
00215     theParticleChange.AddSecondary(aSec);
00216   } else if(erec > 0.0) {
00217     theParticleChange.SetLocalEnergyDeposit(erec);
00218   }
00219 
00220   return &theParticleChange;
00221 }

G4double G4HadronElastic::ComputeMomentumCMS ( const G4ParticleDefinition p,
G4double  plab,
G4int  Z,
G4int  A 
) [inline]

Definition at line 98 of file G4HadronElastic.hh.

References G4NucleiProperties::GetNuclearMass(), and G4ParticleDefinition::GetPDGMass().

Referenced by SampleInvariantT().

00100 {
00101   G4double m1 = p->GetPDGMass();
00102   G4double m12= m1*m1;
00103   G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
00104   return plab*mass2/std::sqrt(m12 + mass2*mass2 + 2.*mass2*std::sqrt(m12 + plab*plab));
00105 }

void G4HadronElastic::Description (  )  const [virtual]

Reimplemented in G4CHIPSElastic, G4ChipsElasticModel, and G4ElasticHadrNucleusHE.

Definition at line 64 of file G4HadronElastic.cc.

References G4HadronicInteraction::GetModelName().

00065 {
00066   char* dirName = getenv("G4PhysListDocDir");
00067   if (dirName) {
00068     std::ofstream outFile;
00069     G4String outFileName = GetModelName() + ".html";
00070     G4String pathName = G4String(dirName) + "/" + outFileName;
00071     outFile.open(pathName);
00072     outFile << "<html>\n";
00073     outFile << "<head>\n";
00074 
00075     outFile << "<title>Description of G4HadronElastic Model</title>\n";
00076     outFile << "</head>\n";
00077     outFile << "<body>\n";
00078 
00079     outFile << "G4HadronElastic is a hadron-nucleus elastic scattering\n"
00080             << "model which uses the Gheisha two-exponential momentum\n"
00081             << "transfer parameterization.  The model is fully relativistic\n"
00082             << "as opposed to the original Gheisha model which was not.\n"
00083             << "This model may be used for all long-lived hadrons at all\n"
00084             << "incident energies.\n";
00085 
00086     outFile << "</body>\n";
00087     outFile << "</html>\n";
00088     outFile.close();
00089   }
00090 }

G4double G4HadronElastic::LowestEnergyLimit (  )  const [inline]

Definition at line 92 of file G4HadronElastic.hh.

00093 {
00094   return lowestEnergyLimit;
00095 }

G4double G4HadronElastic::SampleInvariantT ( const G4ParticleDefinition p,
G4double  plab,
G4int  Z,
G4int  A 
) [virtual]

Reimplemented from G4HadronicInteraction.

Reimplemented in G4AntiNuclElastic, G4CHIPSElastic, G4ChipsElasticModel, G4DiffuseElastic, G4ElasticHadrNucleusHE, and G4NuclNuclDiffuseElastic.

Definition at line 225 of file G4HadronElastic.cc.

References ComputeMomentumCMS(), G4UniformRand, G4Pow::GetInstance(), G4Pow::powZ(), G4Pow::Z13(), and G4Pow::Z23().

Referenced by ApplyYourself(), G4ChipsElasticModel::SampleInvariantT(), and G4CHIPSElastic::SampleInvariantT().

00228 {
00229   static const G4double GeV2 = GeV*GeV;
00230   G4double momentumCMS = ComputeMomentumCMS(p,plab,Z,A);
00231   G4double tmax = 4.0*momentumCMS*momentumCMS/GeV2;
00232   G4double aa, bb, cc;
00233   G4double dd = 10.;
00234   G4Pow* g4pow = G4Pow::GetInstance();
00235   if (A <= 62) {
00236     bb = 14.5*g4pow->Z23(A);
00237     aa = g4pow->powZ(A, 1.63)/bb;
00238     cc = 1.4*g4pow->Z13(A)/dd;
00239   } else {
00240     bb = 60.*g4pow->Z13(A);
00241     aa = g4pow->powZ(A, 1.33)/bb;
00242     cc = 0.4*g4pow->powZ(A, 0.4)/dd;
00243   }
00244   G4double q1 = 1.0 - std::exp(-bb*tmax);
00245   G4double q2 = 1.0 - std::exp(-dd*tmax);
00246   G4double s1 = q1*aa;
00247   G4double s2 = q2*cc;
00248   if((s1 + s2)*G4UniformRand() < s2) {
00249     q1 = q2;
00250     bb = dd;
00251   }
00252   return -GeV2*std::log(1.0 - G4UniformRand()*q1)/bb;
00253 }

void G4HadronElastic::SetLowestEnergyLimit ( G4double  value  )  [inline]

Reimplemented in G4DiffuseElastic, and G4NuclNuclDiffuseElastic.

Definition at line 87 of file G4HadronElastic.hh.

00088 {
00089   lowestEnergyLimit = value;
00090 }


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