Geant4  10.01.p02
G4CascadeRecoilMaker.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: G4CascadeRecoilMaker.cc 66241 2012-12-13 18:34:42Z gunter $
27 //
28 // Collects generated cascade data (using Collider::collide() interface)
29 // and computes the nuclear recoil kinematics needed to balance the event.
30 //
31 // 20100909 M. Kelsey -- Inspired by G4CascadeCheckBalance
32 // 20100909 M. Kelsey -- Move G4IntraNucleiCascader::goodCase() here, add
33 // tolerance for "almost zero" excitation energy
34 // 20100910 M. Kelsey -- Drop getRecoilFragment() in favor of user calling
35 // makeRecoilFragment() with returned non-const pointer. Drop
36 // handling of excitons.
37 // 20100921 M. Kelsey -- Return G4InuclNuclei using "makeRecoilNuclei()".
38 // Repurpose "makeRecoilFragment()" to return G4Fragment.
39 // 20100924 M. Kelsey -- Remove unusable G4Fragment::SetExcitationEnergy().
40 // Add deltaM to compute mass difference, use excitationEnergy
41 // to force G4Fragment four-vector to match.
42 // 20110214 M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
43 // 20110303 M. Kelsey -- Add diagnostic messages to goodNucleus().
44 // 20110308 M. Kelsey -- Follow new G4Fragment interface for hole types
45 // 20110722 M. Kelsey -- For IntraNucleiCascader, take G4CollOut as argument
46 
47 #include <vector>
48 
49 #include "G4CascadeRecoilMaker.hh"
50 #include "globals.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4CascadParticle.hh"
53 #include "G4CascadeCheckBalance.hh"
54 #include "G4CollisionOutput.hh"
55 #include "G4Fragment.hh"
57 #include "G4InuclNuclei.hh"
58 #include "G4InuclParticle.hh"
60 #include "G4LorentzVector.hh"
61 
62 using namespace G4InuclSpecialFunctions;
63 
64 
65 // Constructor and destructor
66 
68  : G4VCascadeCollider("G4CascadeRecoilMaker"),
69  excTolerance(tolerance), inputEkin(0.),
70  recoilA(0), recoilZ(0), excitationEnergy(0.) {
71  balance = new G4CascadeCheckBalance(tolerance, tolerance, theName);
72 }
73 
75  delete balance;
76 }
77 
78 
79 // Standard Collider interface (non-const output "buffer")
80 
82  G4InuclParticle* target,
83  G4CollisionOutput& output) {
84  if (verboseLevel > 1)
85  G4cout << " >>> G4CascadeRecoilMaker::collide" << G4endl;
86 
87  // Available energy needed for "goodNucleus()" test at end
88  inputEkin = bullet ? bullet->getKineticEnergy() : 0.;
89 
91  balance->collide(bullet, target, output);
92  fillRecoil();
93 }
94 
95 // This is for use with G4IntraNucleiCascader
96 
98  G4InuclParticle* target,
99  G4CollisionOutput& output,
100  const std::vector<G4CascadParticle>& cparticles) {
101  if (verboseLevel > 1)
102  G4cout << " >>> G4CascadeRecoilMaker::collide(<EP>,<CP>)" << G4endl;
103 
104  // Available energy needed for "goodNucleus()" test at end
105  inputEkin = bullet ? bullet->getKineticEnergy() : 0.;
106 
108  balance->collide(bullet, target, output, cparticles);
109  fillRecoil();
110 }
111 
112 
113 // Used current event configuration to construct recoil nucleus
114 // NOTE: CheckBalance uses "final-initial", we want "initial-final"
115 
117  recoilZ = -(balance->deltaQ()); // Charge "non-conservation"
118  recoilA = -(balance->deltaB()); // Baryon "non-conservation"
120 
121  theExcitons.clear(); // Discard previous exciton configuraiton
122 
123  // Bertini uses MeV for excitation energy
124  if (!goodFragment()) excitationEnergy = 0.;
125  else excitationEnergy = deltaM() * GeV;
126 
127  // Allow for very small negative mass difference, and round to zero
128  if (std::abs(excitationEnergy) < excTolerance) excitationEnergy = 0.;
129 
130  if (verboseLevel > 2) {
131  G4cout << " recoil px " << recoilMomentum.px()
132  << " py " << recoilMomentum.py() << " pz " << recoilMomentum.pz()
133  << " E " << recoilMomentum.e() << " baryon " << recoilA
134  << " charge " << recoilZ
135  << "\n recoil mass " << recoilMomentum.m()
136  << " 'excitation' energy " << excitationEnergy << G4endl;
137  }
138 }
139 
140 
141 // Construct physical nucleus from recoil parameters, if reasonable
142 
145  if (verboseLevel > 1)
146  G4cout << " >>> G4CascadeRecoilMaker::makeRecoilNuclei" << G4endl;
147 
148  if (!goodRecoil()) {
149  if (verboseLevel > 2 && !wholeEvent())
150  G4cout << theName << ": event recoil is not a physical nucleus" << G4endl;
151 
152  return 0; // Null pointer means no fragment
153  }
154 
156  excitationEnergy, model);
158 
159  return &theRecoilNuclei;
160 }
161 
162 
163 // Construct pre-compound nuclear fragment from recoil parameters
164 
166  if (verboseLevel > 1)
167  G4cout << " >>> G4CascadeRecoilMaker::makeRecoilFragment" << G4endl;
168 
169  if (!goodRecoil()) {
170  if (verboseLevel > 2 && !wholeEvent())
171  G4cout << theName << ": event recoil is not a physical nucleus" << G4endl;
172 
173  return 0; // Null pointer means no fragment
174  }
175 
176  theRecoilFragment.SetZandA_asInt(recoilZ, recoilA); // Note convention!
177 
178  // User may have overridden excitation energy; force four-momentum to match
179  G4double fragMass =
181 
182  G4LorentzVector fragMom; fragMom.setVectM(recoilMomentum.vect(), fragMass);
183  theRecoilFragment.SetMomentum(fragMom*GeV); // Bertini uses GeV!
184 
185  // Note: exciton configuration has to be set piece by piece
186  // (arguments are Ntotal,Nproton in both cases)
189 
194 
195  return &theRecoilFragment;
196 }
197 
198 
199 // Compute raw mass difference from recoil parameters
200 
203  return (recoilMomentum.m() - nucMass);
204 }
205 
206 
207 // Data quality checks
208 
210  return (recoilA>0 && recoilZ>=0 && recoilA >= recoilZ);
211 }
212 
214  return (goodFragment() && excitationEnergy > -excTolerance);
215 }
216 
218  if (verboseLevel > 2) {
219  G4cout << " >>> G4CascadeRecoilMaker::wholeEvent:"
220  << " A " << recoilA << " Z " << recoilZ
221  << " P " << recoilMomentum.rho() << " E " << recoilMomentum.e()
222  << "\n wholeEvent returns "
223  << (recoilA==0 && recoilZ==0 &&
224  recoilMomentum.rho() < excTolerance/GeV &&
225  std::abs(recoilMomentum.e()) < excTolerance/GeV) << G4endl;
226  }
227 
228  return (recoilA==0 && recoilZ==0 &&
229  recoilMomentum.rho() < excTolerance/GeV &&
230  std::abs(recoilMomentum.e()) < excTolerance/GeV);
231 }
232 
233 // Determine whether desired nuclear fragment is constructable outcome
234 
236  if (verboseLevel > 2) {
237  G4cout << " >>> G4CascadeRecoilMaker::goodNucleus" << G4endl;
238  }
239 
240  const G4double minExcitation = 0.1*keV;
241  const G4double reasonableExcitation = 7.0; // Multiple of binding energy
242  const G4double fractionalExcitation = 0.2; // Fraction of input to excite
243 
244  if (!goodRecoil()) {
245  if (verboseLevel>2) {
246  if (!goodFragment()) G4cerr << " goodNucleus: invalid A/Z" << G4endl;
247  else if (excitationEnergy < -excTolerance)
248  G4cerr << " goodNucleus: negative excitation" << G4endl;
249  }
250  return false; // Not a sensible nucleus
251  }
252 
253  if (excitationEnergy <= minExcitation) return true; // Effectively zero
254 
255  // Maximum possible excitation energy determined by initial energy
257  G4double exc_max0z = fractionalExcitation * inputEkin*GeV;
258  G4double exc_dm = reasonableExcitation * dm;
259  G4double exc_max = (exc_max0z > exc_dm) ? exc_max0z : exc_dm;
260 
261  if (verboseLevel > 3) {
262  G4cout << " eexs " << excitationEnergy << " max " << exc_max
263  << " dm " << dm << G4endl;
264  }
265 
266  if (verboseLevel > 2 && excitationEnergy >= exc_max)
267  G4cerr << " goodNucleus: too much excitation" << G4endl;
268 
269  return (excitationEnergy < exc_max); // Below maximum possible
270 }
void fill(G4int a, G4int z, G4double exc=0., Model model=DefaultModel)
G4LorentzVector deltaLV() const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
Definition: G4Fragment.hh:360
static const G4double tolerance
int G4int
Definition: G4Types.hh:78
virtual void setVerboseLevel(G4int verbose=0)
G4double getKineticEnergy() const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
void SetNumberOfExcitedParticle(G4int valueTot, G4int valueP)
Definition: G4Fragment.hh:341
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:281
G4double getNucleiMass() const
static const double GeV
Definition: G4SIunits.hh:196
G4InuclNuclei * makeRecoilNuclei(G4InuclParticle::Model model=G4InuclParticle::DefaultModel)
void setExitonConfiguration(const G4ExitonConfiguration &config)
G4ExitonConfiguration theExcitons
void SetZandA_asInt(G4int Znew, G4int Anew)
Definition: G4Fragment.hh:253
#define G4endl
Definition: G4ios.hh:61
G4CascadeCheckBalance * balance
static const double keV
Definition: G4SIunits.hh:195
double G4double
Definition: G4Types.hh:76
G4LorentzVector recoilMomentum
G4CascadeRecoilMaker(G4double tolerance=0.001 *CLHEP::MeV)
G4double bindingEnergy(G4int A, G4int Z)
G4Fragment * makeRecoilFragment()
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector