Geant4  10.01.p03
G4BigBanger.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: G4BigBanger.cc 71942 2013-06-28 19:08:11Z mkelsey $
27 //
28 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100301 M. Kelsey -- In generateBangInSCM(), restore old G4CascMom calcs.
30 // for (N-1)th outgoing nucleon.
31 // 20100319 M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff
32 // 20100407 M. Kelsey -- Replace std::vector<> returns with data members.
33 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
34 // 20100517 M. Kelsey -- Inherit from common base class, clean up code
35 // 20100628 M. Kelsey -- Use old "bindingEnergy" fn as wrapper, add balance
36 // checking after bang.
37 // 20100630 M. Kelsey -- Just do simple boost for target, instead of using
38 // G4LorentzConverter with dummy bullet.
39 // 20100701 M. Kelsey -- Re-throw momentum list, not just angles!
40 // 20100714 M. Kelsey -- Move conservation checking to base class
41 // 20100726 M. Kelsey -- Move std::vector<> buffer to .hh file
42 // 20100923 M. Kelsey -- Migrate to integer A and Z
43 // 20110214 M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
44 // 20110806 M. Kelsey -- Pre-allocate buffers to reduce memory churn
45 // 20110922 M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
46 // 20120608 M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
47 // 20130622 Inherit from G4CascadeDeexciteBase, move to deExcite() interface
48 // with G4Fragment
49 // 20130924 M. Kelsey -- Replace std::pow with G4Pow::powN() for CPU speed
50 
51 #include <algorithm>
52 
53 #include "G4BigBanger.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4CollisionOutput.hh"
56 #include "G4InuclNuclei.hh"
59 #include "G4ParticleLargerEkin.hh"
60 #include "G4Pow.hh"
61 
62 using namespace G4InuclSpecialFunctions;
63 
64 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
65 
67 
68 void G4BigBanger::deExcite(const G4Fragment& target,
69  G4CollisionOutput& globalOutput) {
70  if (verboseLevel) G4cout << " >>> G4BigBanger::deExcite" << G4endl;
71 
72  getTargetData(target);
73  G4ThreeVector toTheLabFrame = PEX.boostVector(); // From rest to lab
74 
75  // This "should" be difference between E-target and sum of m(nucleons)
76  G4double etot = (EEXS - bindingEnergy(A,Z)) * MeV/GeV; // To Bertini units
77  if (etot < 0.0) etot = 0.0;
78 
79  if (verboseLevel > 2) {
80  G4cout << " BigBanger: target\n" << target
81  << "\n etot " << etot << G4endl;
82  }
83 
84  if (verboseLevel > 3) {
85  G4LorentzVector PEXrest = PEX;
86  PEXrest.boost(-toTheLabFrame);
87  G4cout << " target rest frame: px " << PEXrest.px() << " py "
88  << PEXrest.py() << " pz " << PEXrest.pz() << " E " << PEXrest.e()
89  << G4endl;
90  }
91 
92  generateBangInSCM(etot, A, Z);
93 
94  if (verboseLevel > 2) {
95  G4cout << " particles " << particles.size() << G4endl;
96  for(G4int i = 0; i < G4int(particles.size()); i++)
97  G4cout << particles[i] << G4endl;
98  }
99 
100  if (particles.empty()) { // No bang! Don't know why...
101  G4cerr << " >>> G4BigBanger unable to process fragment "
102  << target << G4endl;
103 
104  // FIXME: This will violate baryon number, momentum, energy, etc.
105  return;
106  }
107 
108  // convert back to Lab
109  G4LorentzVector totscm;
110  G4LorentzVector totlab;
111 
112  if (verboseLevel > 2) G4cout << " BigBanger: boosting to lab" << G4endl;
113 
114  particleIterator ipart;
115  for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
116  G4LorentzVector mom = ipart->getMomentum();
117  if (verboseLevel > 2) totscm += mom;
118 
119  mom.boost(toTheLabFrame);
120  if (verboseLevel > 2) totlab += mom;
121 
122  ipart->setMomentum(mom);
123  if (verboseLevel > 2) G4cout << *ipart << G4endl;
124  }
125 
126  std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
127 
128  validateOutput(target, particles); // Checks <vector> directly
129 
130  if (verboseLevel > 2) {
131  G4cout << " In SCM: total outgoing momentum " << G4endl
132  << " E " << totscm.e() << " px " << totscm.x()
133  << " py " << totscm.y() << " pz " << totscm.z() << G4endl;
134  G4cout << " In Lab: mom cons " << G4endl
135  << " E " << PEX.e() - totlab.e() // PEX now includes EEXS
136  << " px " << PEX.x() - totlab.x()
137  << " py " << PEX.y() - totlab.y()
138  << " pz " << PEX.z() - totlab.z() << G4endl;
139  }
140 
141  globalOutput.addOutgoingParticles(particles);
142 }
143 
145  if (verboseLevel > 3) {
146  G4cout << " >>> G4BigBanger::generateBangInSCM" << G4endl;
147  }
148 
149  const G4double ang_cut = 0.9999;
150  const G4int itry_max = 1000;
151 
152  if (verboseLevel > 2) {
153  G4cout << " a " << a << " z " << z << G4endl;
154  }
155 
156  particles.clear(); // Reset output vector before filling
157 
158  if (a == 1) { // Special -- bare nucleon doesn't really "explode"
159  G4int knd = (z>0) ? 1 : 2;
160  particles.push_back(G4InuclElementaryParticle(knd)); // zero momentum
161  return;
162  }
163 
164  // NOTE: If distribution fails, need to regenerate magnitudes and angles!
165  //*** generateMomentumModules(etot, a, z);
166 
167  scm_momentums.reserve(a);
168  G4LorentzVector tot_mom;
169 
170  G4bool bad = true;
171  G4int itry = 0;
172  while(bad && itry < itry_max) {
173  itry++;
174  scm_momentums.clear();
175 
176  generateMomentumModules(etot, a, z);
177  if (a == 2) {
178  // This is only a three-vector, not a four-vector
180  scm_momentums.push_back(mom);
181  scm_momentums.push_back(-mom); // Only safe since three-vector!
182  bad = false;
183  } else {
184  tot_mom *= 0.; // Easy way to reset accumulator
185 
186  for(G4int i = 0; i < a-2; i++) { // All but last two are thrown
187  // This is only a three-vector, not a four-vector
189  scm_momentums.push_back(mom);
190  tot_mom += mom;
191  };
192 
193  // handle last two
194  G4double tot_mod = tot_mom.rho();
195  G4double ct = -0.5*(tot_mod*tot_mod + momModules[a-2]*momModules[a-2]
196  - momModules[a-1]*momModules[a-1]) / tot_mod
197  / momModules[a-2];
198 
199  if (verboseLevel > 2) G4cout << " ct last " << ct << G4endl;
200 
201  if(std::fabs(ct) < ang_cut) {
202  // This is only a three-vector, not a four-vector
204 
205  // rotate to the normal system
206  G4LorentzVector apr = tot_mom/tot_mod;
207  G4double a_tr = std::sqrt(apr.x()*apr.x() + apr.y()*apr.y());
208  G4LorentzVector mom;
209  mom.setX(mom2.z()*apr.x() + ( mom2.x()*apr.y() + mom2.y()*apr.z()*apr.x())/a_tr);
210  mom.setY(mom2.z()*apr.y() + (-mom2.x()*apr.x() + mom2.y()*apr.z()*apr.y())/a_tr);
211  mom.setZ(mom2.z()*apr.z() - mom2.y()*a_tr);
212 
213  scm_momentums.push_back(mom);
214 
215  // and the last one (again, not actually a four-vector!)
216  G4LorentzVector mom1 = -mom - tot_mom;
217 
218  scm_momentums.push_back(mom1);
219  bad = false;
220  } // if (abs(ct) < ang_cut)
221  } // (a > 2)
222  } // while (bad && itry<itry_max)
223 
224  if (!bad) {
225  particles.resize(a); // Use assignment to avoid temporaries
226  for(G4int i = 0; i < a; i++) {
227  G4int knd = i < z ? 1 : 2;
229  };
230  };
231 
232  if (verboseLevel > 2) {
233  if (itry == itry_max) G4cout << " BigBanger -> can not generate bang " << G4endl;
234  }
235 
236  return;
237 }
238 
240  if (verboseLevel > 3) {
241  G4cout << " >>> G4BigBanger::generateMomentumModules" << G4endl;
242  }
243 
244  // Proton and neutron masses
247 
248  momModules.clear(); // Reset buffer for filling
249 
250  G4double xtot = 0.0;
251 
252  if (a > 2) { // For "large" nuclei, energy is distributed
253  G4double promax = maxProbability(a);
254 
255  momModules.resize(a, 0.); // Pre-allocate to avoid memory churn
256  for(G4int i = 0; i < a; i++) {
257  momModules[i] = generateX(a, promax);
258  xtot += momModules[i];
259 
260  if (verboseLevel > 2) {
261  G4cout << " i " << i << " x " << momModules[i] << G4endl;
262  }
263  }
264  } else { // Two-body case is special, must be 50%
265  xtot = 1.;
266  momModules.push_back(0.5);
267  momModules.push_back(0.5);
268  }
269 
270  for(G4int i = 0; i < a; i++) {
271  G4double mass = i < z ? mp : mn;
272 
273  momModules[i] *= etot/xtot;
274  momModules[i] = std::sqrt(momModules[i] * (momModules[i] + 2.0 * mass));
275 
276  if (verboseLevel > 2) {
277  G4cout << " i " << i << " pmod " << momModules[i] << G4endl;
278  }
279  };
280 
281  return;
282 }
283 
285  if (verboseLevel > 3) G4cout << " >>> G4BigBanger::xProbability" << G4endl;
286 
287  G4Pow* theG4Pow = G4Pow::GetInstance(); // For convenience
288 
289  G4double ekpr = 0.0;
290  if(x < 1.0 || x > 0.0) {
291  ekpr = x * x;
292 
293  if (a%2 == 0) { // even A
294  ekpr *= std::sqrt(1.0 - x) * theG4Pow->powN((1.0 - x), (3*a-6)/2);
295  }
296  else {
297  ekpr *= theG4Pow->powN((1.0 - x), (3*a-5)/2);
298  };
299  };
300 
301  return ekpr;
302 }
303 
305  if (verboseLevel > 3) {
306  G4cout << " >>> G4BigBanger::maxProbability" << G4endl;
307  }
308 
309  return xProbability(2./3./(a-1.0), a);
310 }
311 
313  if (verboseLevel > 3) G4cout << " >>> G4BigBanger::generateX" << G4endl;
314 
315  const G4int itry_max = 1000;
316  G4int itry = 0;
317  G4double x;
318 
319  while(itry < itry_max) {
320  itry++;
321  x = inuclRndm();
322 
323  if(xProbability(x, a) >= promax * inuclRndm()) return x;
324  };
325  if (verboseLevel > 2) {
326  G4cout << " BigBanger -> can not generate x " << G4endl;
327  }
328 
329  return maxProbability(a);
330 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double xProbability(G4double x, G4int a) const
Definition: G4BigBanger.cc:284
static const double MeV
Definition: G4SIunits.hh:193
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:128
CLHEP::Hep3Vector G4ThreeVector
Definition: G4Pow.hh:56
std::vector< G4InuclElementaryParticle > particles
Definition: G4BigBanger.hh:63
G4double z
Definition: TRTMaterials.hh:39
G4double a
Definition: TRTMaterials.hh:39
static G4double getParticleMass(G4int type)
int G4int
Definition: G4Types.hh:78
void generateMomentumModules(G4double etot, G4int a, G4int z)
Definition: G4BigBanger.cc:239
void getTargetData(const G4Fragment &target)
G4GLOB_DLL std::ostream G4cout
std::vector< G4LorentzVector > scm_momentums
Definition: G4BigBanger.hh:65
std::vector< G4double > momModules
Definition: G4BigBanger.hh:64
bool G4bool
Definition: G4Types.hh:79
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
static const double GeV
Definition: G4SIunits.hh:196
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
G4LorentzVector generateWithFixedTheta(G4double ct, G4double p, G4double mass=0.)
G4double maxProbability(G4int a) const
Definition: G4BigBanger.cc:304
void generateBangInSCM(G4double etot, G4int a, G4int z)
Definition: G4BigBanger.cc:144
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:64
#define G4endl
Definition: G4ios.hh:61
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
G4double generateX(G4int ia, G4double promax) const
Definition: G4BigBanger.cc:312
double G4double
Definition: G4Types.hh:76
G4double bindingEnergy(G4int A, G4int Z)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
Definition: G4BigBanger.cc:68
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector