Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 // 20150608 M. Kelsey -- Label all while loops as terminating.
51 
52 #include <algorithm>
53 
54 #include "G4BigBanger.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4CollisionOutput.hh"
57 #include "G4InuclNuclei.hh"
60 #include "G4ParticleLargerEkin.hh"
61 #include "G4Pow.hh"
62 
63 using namespace G4InuclSpecialFunctions;
64 
65 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
66 
68 
70  G4CollisionOutput& globalOutput) {
71  if (verboseLevel) G4cout << " >>> G4BigBanger::deExcite" << G4endl;
72 
73  getTargetData(target);
74  G4ThreeVector toTheLabFrame = PEX.boostVector(); // From rest to lab
75 
76  // This "should" be difference between E-target and sum of m(nucleons)
77  G4double etot = (EEXS - bindingEnergy(A,Z)) * MeV/GeV; // To Bertini units
78  if (etot < 0.0) etot = 0.0;
79 
80  if (verboseLevel > 2) {
81  G4cout << " BigBanger: target\n" << target
82  << "\n etot " << etot << G4endl;
83  }
84 
85  if (verboseLevel > 3) {
86  G4LorentzVector PEXrest = PEX;
87  PEXrest.boost(-toTheLabFrame);
88  G4cout << " target rest frame: px " << PEXrest.px() << " py "
89  << PEXrest.py() << " pz " << PEXrest.pz() << " E " << PEXrest.e()
90  << G4endl;
91  }
92 
93  generateBangInSCM(etot, A, Z);
94 
95  if (verboseLevel > 2) {
96  G4cout << " particles " << particles.size() << G4endl;
97  for(G4int i = 0; i < G4int(particles.size()); i++)
98  G4cout << particles[i] << G4endl;
99  }
100 
101  if (particles.empty()) { // No bang! Don't know why...
102  G4cerr << " >>> G4BigBanger unable to process fragment "
103  << target << G4endl;
104 
105  // FIXME: This will violate baryon number, momentum, energy, etc.
106  return;
107  }
108 
109  // convert back to Lab
110  G4LorentzVector totscm;
111  G4LorentzVector totlab;
112 
113  if (verboseLevel > 2) G4cout << " BigBanger: boosting to lab" << G4endl;
114 
115  particleIterator ipart;
116  for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
117  G4LorentzVector mom = ipart->getMomentum();
118  if (verboseLevel > 2) totscm += mom;
119 
120  mom.boost(toTheLabFrame);
121  if (verboseLevel > 2) totlab += mom;
122 
123  ipart->setMomentum(mom);
124  if (verboseLevel > 2) G4cout << *ipart << G4endl;
125  }
126 
127  std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
128 
129  validateOutput(target, particles); // Checks <vector> directly
130 
131  if (verboseLevel > 2) {
132  G4cout << " In SCM: total outgoing momentum " << G4endl
133  << " E " << totscm.e() << " px " << totscm.x()
134  << " py " << totscm.y() << " pz " << totscm.z() << G4endl;
135  G4cout << " In Lab: mom cons " << G4endl
136  << " E " << PEX.e() - totlab.e() // PEX now includes EEXS
137  << " px " << PEX.x() - totlab.x()
138  << " py " << PEX.y() - totlab.y()
139  << " pz " << PEX.z() - totlab.z() << G4endl;
140  }
141 
142  globalOutput.addOutgoingParticles(particles);
143 }
144 
145 void G4BigBanger::generateBangInSCM(G4double etot, G4int a, G4int z) {
146  if (verboseLevel > 3) {
147  G4cout << " >>> G4BigBanger::generateBangInSCM" << G4endl;
148  }
149 
150  const G4double ang_cut = 0.9999;
151  const G4int itry_max = 1000;
152 
153  if (verboseLevel > 2) {
154  G4cout << " a " << a << " z " << z << G4endl;
155  }
156 
157  particles.clear(); // Reset output vector before filling
158 
159  if (a == 1) { // Special -- bare nucleon doesn't really "explode"
160  G4int knd = (z>0) ? 1 : 2;
161  particles.push_back(G4InuclElementaryParticle(knd)); // zero momentum
162  return;
163  }
164 
165  // NOTE: If distribution fails, need to regenerate magnitudes and angles!
166  //*** generateMomentumModules(etot, a, z);
167 
168  scm_momentums.reserve(a);
169  G4LorentzVector tot_mom;
170 
171  G4bool bad = true;
172  G4int itry = 0;
173  while(bad && itry < itry_max) { /* Loop checking 08.06.2015 MHK */
174  itry++;
175  scm_momentums.clear();
176 
177  generateMomentumModules(etot, a, z);
178  if (a == 2) {
179  // This is only a three-vector, not a four-vector
180  G4LorentzVector mom = generateWithRandomAngles(momModules[0]);
181  scm_momentums.push_back(mom);
182  scm_momentums.push_back(-mom); // Only safe since three-vector!
183  bad = false;
184  } else {
185  tot_mom *= 0.; // Easy way to reset accumulator
186 
187  for(G4int i = 0; i < a-2; i++) { // All but last two are thrown
188  // This is only a three-vector, not a four-vector
189  G4LorentzVector mom = generateWithRandomAngles(momModules[i]);
190  scm_momentums.push_back(mom);
191  tot_mom += mom;
192  };
193 
194  // handle last two
195  G4double tot_mod = tot_mom.rho();
196  G4double ct = -0.5*(tot_mod*tot_mod + momModules[a-2]*momModules[a-2]
197  - momModules[a-1]*momModules[a-1]) / tot_mod
198  / momModules[a-2];
199 
200  if (verboseLevel > 2) G4cout << " ct last " << ct << G4endl;
201 
202  if(std::fabs(ct) < ang_cut) {
203  // This is only a three-vector, not a four-vector
204  G4LorentzVector mom2 = generateWithFixedTheta(ct, momModules[a - 2]);
205 
206  // rotate to the normal system
207  G4LorentzVector apr = tot_mom/tot_mod;
208  G4double a_tr = std::sqrt(apr.x()*apr.x() + apr.y()*apr.y());
209  G4LorentzVector mom;
210  mom.setX(mom2.z()*apr.x() + ( mom2.x()*apr.y() + mom2.y()*apr.z()*apr.x())/a_tr);
211  mom.setY(mom2.z()*apr.y() + (-mom2.x()*apr.x() + mom2.y()*apr.z()*apr.y())/a_tr);
212  mom.setZ(mom2.z()*apr.z() - mom2.y()*a_tr);
213 
214  scm_momentums.push_back(mom);
215 
216  // and the last one (again, not actually a four-vector!)
217  G4LorentzVector mom1 = -mom - tot_mom;
218 
219  scm_momentums.push_back(mom1);
220  bad = false;
221  } // if (abs(ct) < ang_cut)
222  } // (a > 2)
223  } // while (bad && itry<itry_max)
224 
225  if (!bad) {
226  particles.resize(a); // Use assignment to avoid temporaries
227  for(G4int i = 0; i < a; i++) {
228  G4int knd = i < z ? 1 : 2;
229  particles[i].fill(scm_momentums[i], knd, G4InuclParticle::BigBanger);
230  };
231  };
232 
233  if (verboseLevel > 2) {
234  if (itry == itry_max) G4cout << " BigBanger -> can not generate bang " << G4endl;
235  }
236 
237  return;
238 }
239 
240 void G4BigBanger::generateMomentumModules(G4double etot, G4int a, G4int z) {
241  if (verboseLevel > 3) {
242  G4cout << " >>> G4BigBanger::generateMomentumModules" << G4endl;
243  }
244 
245  // Proton and neutron masses
248 
249  momModules.clear(); // Reset buffer for filling
250 
251  G4double xtot = 0.0;
252 
253  if (a > 2) { // For "large" nuclei, energy is distributed
254  G4double promax = maxProbability(a);
255 
256  momModules.resize(a, 0.); // Pre-allocate to avoid memory churn
257  for(G4int i = 0; i < a; i++) {
258  momModules[i] = generateX(a, promax);
259  xtot += momModules[i];
260 
261  if (verboseLevel > 2) {
262  G4cout << " i " << i << " x " << momModules[i] << G4endl;
263  }
264  }
265  } else { // Two-body case is special, must be 50%
266  xtot = 1.;
267  momModules.push_back(0.5);
268  momModules.push_back(0.5);
269  }
270 
271  for(G4int i = 0; i < a; i++) {
272  G4double mass = i < z ? mp : mn;
273 
274  momModules[i] *= etot/xtot;
275  momModules[i] = std::sqrt(momModules[i] * (momModules[i] + 2.0 * mass));
276 
277  if (verboseLevel > 2) {
278  G4cout << " i " << i << " pmod " << momModules[i] << G4endl;
279  }
280  };
281 
282  return;
283 }
284 
285 G4double G4BigBanger::xProbability(G4double x, G4int a) const {
286  if (verboseLevel > 3) G4cout << " >>> G4BigBanger::xProbability" << G4endl;
287 
288  G4Pow* theG4Pow = G4Pow::GetInstance(); // For convenience
289 
290  G4double ekpr = 0.0;
291  if(x < 1.0 || x > 0.0) {
292  ekpr = x * x;
293 
294  if (a%2 == 0) { // even A
295  ekpr *= std::sqrt(1.0 - x) * theG4Pow->powN((1.0 - x), (3*a-6)/2);
296  }
297  else {
298  ekpr *= theG4Pow->powN((1.0 - x), (3*a-5)/2);
299  };
300  };
301 
302  return ekpr;
303 }
304 
305 G4double G4BigBanger::maxProbability(G4int a) const {
306  if (verboseLevel > 3) {
307  G4cout << " >>> G4BigBanger::maxProbability" << G4endl;
308  }
309 
310  return xProbability(2./3./(a-1.0), a);
311 }
312 
313 G4double G4BigBanger::generateX(G4int a, G4double promax) const {
314  if (verboseLevel > 3) G4cout << " >>> G4BigBanger::generateX" << G4endl;
315 
316  const G4int itry_max = 1000;
317  G4int itry = 0;
318  G4double x;
319 
320  while(itry < itry_max) { /* Loop checking 08.06.2015 MHK */
321  itry++;
322  x = inuclRndm();
323 
324  if(xProbability(x, a) >= promax * inuclRndm()) return x;
325  };
326  if (verboseLevel > 2) {
327  G4cout << " BigBanger -> can not generate x " << G4endl;
328  }
329 
330  return maxProbability(a);
331 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
Hep3Vector boostVector() const
const XML_Char * target
Definition: expat.h:268
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:128
Definition: G4Pow.hh:56
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple x
Definition: test.py:50
static G4double getParticleMass(G4int type)
int G4int
Definition: G4Types.hh:78
void getTargetData(const G4Fragment &target)
G4GLOB_DLL std::ostream G4cout
double py() const
bool G4bool
Definition: G4Types.hh:79
HepLorentzVector & boost(double, double, double)
double px() const
double rho() const
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
G4LorentzVector generateWithFixedTheta(G4double ct, G4double p, G4double mass=0.)
static constexpr double GeV
Definition: G4SIunits.hh:217
tuple z
Definition: test.py:28
double pz() const
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
double G4double
Definition: G4Types.hh:76
G4double bindingEnergy(G4int A, G4int Z)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
Definition: G4BigBanger.cc:69
G4GLOB_DLL std::ostream G4cerr