Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HadronElastic.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id$
27 //
28 // Geant4 Header : G4HadronElastic
29 //
30 // Author : V.Ivanchenko 29 June 2009 (redesign old elastic model)
31 //
32 
33 #include "G4HadronElastic.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4ParticleTable.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4IonTable.hh"
38 #include "Randomize.hh"
39 #include "G4Proton.hh"
40 #include "G4Neutron.hh"
41 #include "G4Deuteron.hh"
42 #include "G4Alpha.hh"
43 #include "G4Pow.hh"
44 
45 G4HadronElastic::G4HadronElastic(const G4String& name)
46  : G4HadronicInteraction(name)
47 {
48  SetMinEnergy( 0.0*GeV );
49  SetMaxEnergy( 100.*TeV );
50  lowestEnergyLimit= 1.e-6*eV;
51 
52  theProton = G4Proton::Proton();
53  theNeutron = G4Neutron::Neutron();
54  theDeuteron = G4Deuteron::Deuteron();
55  theAlpha = G4Alpha::Alpha();
56  //Description();
57 }
58 
59 
60 G4HadronElastic::~G4HadronElastic()
61 {}
62 
63 
64 void G4HadronElastic::Description() const
65 {
66  char* dirName = getenv("G4PhysListDocDir");
67  if (dirName) {
68  std::ofstream outFile;
69  G4String outFileName = GetModelName() + ".html";
70  G4String pathName = G4String(dirName) + "/" + outFileName;
71  outFile.open(pathName);
72  outFile << "<html>\n";
73  outFile << "<head>\n";
74 
75  outFile << "<title>Description of G4HadronElastic Model</title>\n";
76  outFile << "</head>\n";
77  outFile << "<body>\n";
78 
79  outFile << "G4HadronElastic is a hadron-nucleus elastic scattering\n"
80  << "model which uses the Gheisha two-exponential momentum\n"
81  << "transfer parameterization. The model is fully relativistic\n"
82  << "as opposed to the original Gheisha model which was not.\n"
83  << "This model may be used for all long-lived hadrons at all\n"
84  << "incident energies.\n";
85 
86  outFile << "</body>\n";
87  outFile << "</html>\n";
88  outFile.close();
89  }
90 }
91 
92 
93 G4HadFinalState* G4HadronElastic::ApplyYourself(
94  const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
95 {
96  theParticleChange.Clear();
97 
98  const G4HadProjectile* aParticle = &aTrack;
99  G4double ekin = aParticle->GetKineticEnergy();
100  if(ekin <= lowestEnergyLimit) {
101  theParticleChange.SetEnergyChange(ekin);
102  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
103  return &theParticleChange;
104  }
105 
106  G4int A = targetNucleus.GetA_asInt();
107  G4int Z = targetNucleus.GetZ_asInt();
108 
109  G4double plab = aParticle->GetTotalMomentum();
110 
111  // Scattered particle referred to axis of incident particle
112  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
113  G4double m1 = theParticle->GetPDGMass();
114 
115  if (verboseLevel>1) {
116  G4cout << "G4HadronElastic: "
117  << aParticle->GetDefinition()->GetParticleName()
118  << " Plab(GeV/c)= " << plab/GeV
119  << " Ekin(MeV) = " << ekin/MeV
120  << " scattered off Z= " << Z
121  << " A= " << A
122  << G4endl;
123  }
124 
125  G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
126  G4LorentzVector lv1 = aParticle->Get4Momentum();
127  G4LorentzVector lv(0.0,0.0,0.0,mass2);
128  lv += lv1;
129 
130  G4ThreeVector bst = lv.boostVector();
131  lv1.boost(-bst);
132 
133  G4ThreeVector p1 = lv1.vect();
134  G4double momentumCMS = p1.mag();
135  G4double tmax = 4.0*momentumCMS*momentumCMS;
136 
137  // Sampling in CM system
138  G4double t = SampleInvariantT(theParticle, plab, Z, A);
139  G4double phi = G4UniformRand()*CLHEP::twopi;
140  G4double cost = 1. - 2.0*t/tmax;
141  G4double sint;
142 
143  // problem in sampling
144  if(cost > 1.0 || cost < -1.0) {
145  //if(verboseLevel > 0) {
146  G4cout << "G4HadronElastic WARNING (1 - cost)= " << 1 - cost
147  << " after scattering of "
148  << aParticle->GetDefinition()->GetParticleName()
149  << " p(GeV/c)= " << plab/GeV
150  << " on an ion Z= " << Z << " A= " << A
151  << G4endl;
152  //}
153  cost = 1.0;
154  sint = 0.0;
155 
156  // normal situation
157  } else {
158  sint = std::sqrt((1.0-cost)*(1.0+cost));
159  }
160  if (verboseLevel>1) {
161  G4cout << " t= " << t << " tmax(GeV^2)= " << tmax/(GeV*GeV)
162  << " Pcms(GeV)= " << momentumCMS/GeV << " cos(t)=" << cost
163  << " sin(t)=" << sint << G4endl;
164  }
165  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
166  v1 *= momentumCMS;
167  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),
168  std::sqrt(momentumCMS*momentumCMS + m1*m1));
169 
170  nlv1.boost(bst);
171 
172  G4double eFinal = nlv1.e() - m1;
173  if (verboseLevel > 1) {
174  G4cout <<" m= " << m1 << " Efin(MeV)= " << eFinal
175  << " Proj: 4-mom " << lv1 << " Final: " << nlv1
176  << G4endl;
177  }
178  if(eFinal <= lowestEnergyLimit) {
179  if(eFinal < 0.0 && verboseLevel > 0) {
180  G4cout << "G4HadronElastic WARNING Efinal= " << eFinal
181  << " after scattering of "
182  << aParticle->GetDefinition()->GetParticleName()
183  << " p(GeV/c)= " << plab/GeV
184  << " on an ion Z= " << Z << " A= " << A
185  << G4endl;
186  }
187  theParticleChange.SetEnergyChange(0.0);
188  nlv1 = G4LorentzVector(0.0,0.0,0.0,m1);
189 
190  } else {
191  theParticleChange.SetMomentumChange(nlv1.vect().unit());
192  theParticleChange.SetEnergyChange(eFinal);
193  }
194 
195  lv -= nlv1;
196  G4double erec = lv.e() - mass2;
197  if (verboseLevel > 1) {
198  G4cout << "Recoil: " <<" m= " << mass2 << " Erec(MeV)= " << erec
199  << " 4-mom: " << lv
200  << G4endl;
201  }
202 
203  if(erec > GetRecoilEnergyThreshold()) {
204  G4ParticleDefinition * theDef = 0;
205  if(Z == 1 && A == 1) { theDef = theProton; }
206  else if (Z == 1 && A == 2) { theDef = theDeuteron; }
207  else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
208  else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
209  else if (Z == 2 && A == 4) { theDef = theAlpha; }
210  else {
211  theDef =
212  G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
213  }
214  G4DynamicParticle * aSec = new G4DynamicParticle(theDef, lv);
215  theParticleChange.AddSecondary(aSec);
216  } else if(erec > 0.0) {
217  theParticleChange.SetLocalEnergyDeposit(erec);
218  }
219 
220  return &theParticleChange;
221 }
222 
223 // sample momentum transfer in the CMS system
224 G4double
225 G4HadronElastic::SampleInvariantT(const G4ParticleDefinition* p,
226  G4double plab,
227  G4int Z, G4int A)
228 {
229  static const G4double GeV2 = GeV*GeV;
230  G4double momentumCMS = ComputeMomentumCMS(p,plab,Z,A);
231  G4double tmax = 4.0*momentumCMS*momentumCMS/GeV2;
232  G4double aa, bb, cc;
233  G4double dd = 10.;
234  G4Pow* g4pow = G4Pow::GetInstance();
235  if (A <= 62) {
236  bb = 14.5*g4pow->Z23(A);
237  aa = g4pow->powZ(A, 1.63)/bb;
238  cc = 1.4*g4pow->Z13(A)/dd;
239  } else {
240  bb = 60.*g4pow->Z13(A);
241  aa = g4pow->powZ(A, 1.33)/bb;
242  cc = 0.4*g4pow->powZ(A, 0.4)/dd;
243  }
244  G4double q1 = 1.0 - std::exp(-bb*tmax);
245  G4double q2 = 1.0 - std::exp(-dd*tmax);
246  G4double s1 = q1*aa;
247  G4double s2 = q2*cc;
248  if((s1 + s2)*G4UniformRand() < s2) {
249  q1 = q2;
250  bb = dd;
251  }
252  return -GeV2*std::log(1.0 - G4UniformRand()*q1)/bb;
253 }
254 
Generated on Sat May 25 2013 14:33:51 for Geant4 by   doxygen 1.8.4