Geant4  10.00.p01
G4CascadeCheckBalance.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: G4CascadeCheckBalance.cc 71942 2013-06-28 19:08:11Z mkelsey $
27 //
28 // Verify and report four-momentum conservation for collision output; uses
29 // same interface as collision generators.
30 //
31 // 20100624 M. Kelsey -- Add baryon number, charge, and kinetic energy
32 // 20100624 M. Kelsey -- Bug fix: All checks should be |delta| !
33 // 20100627 M. Kelsey -- Always report violations on cerr, regardless of
34 // verbosity.
35 // 20100628 M. Kelsey -- Add interface to take list of particles directly,
36 // bug fix reporting of charge conservation error.
37 // 20100630 M. Kelsey -- for nuclei, include excitation energies in total.
38 // 20100701 M. Kelsey -- Undo previous change, handled by G4InuclNuclei.
39 // 20100711 M. Kelsey -- Use name of parent collider for reporting messages
40 // 20100713 M. kelsey -- Hide conservation errors behind verbosity
41 // 20100715 M. Kelsey -- Use new G4CollisionOutput totals instead of loops,
42 // move temporary buffer to be data member
43 // 20100719 M. Kelsey -- Change zero tolerance to 10 keV instead of 1 keV.
44 // 20100909 M. Kelsey -- Add interface for both kinds of particle lists
45 // 20101019 M. Kelsey -- CoVerity report: unitialized constructor
46 // 20110328 M. Kelsey -- Add default ctor and explicit limit setting
47 // 20110722 M. Kelsey -- For IntraNucleiCascader, take G4CollOut as argument
48 // 20120525 M. Kelsey -- Follow process-level checking: allow _either_ rel.
49 // or abs. to pass (instead of requiring both)
50 // 20121002 M. Kelsey -- Add strangeness check (useful for Omega- beam)
51 // 20130621 Add interface to take G4Fragment input instead of G4InuclNuclei.
52 
53 #include "G4CascadeCheckBalance.hh"
54 #include "globals.hh"
55 #include "G4CascadParticle.hh"
57 #include "G4InuclNuclei.hh"
58 #include "G4InuclParticle.hh"
59 #include "G4CollisionOutput.hh"
60 #include "G4LorentzVector.hh"
61 #include <vector>
62 
63 
64 // Constructor sets acceptance limits
65 
66 const G4double G4CascadeCheckBalance::tolerance = 1e-6; // How small is zero?
67 
69  : G4VCascadeCollider(owner), relativeLimit(G4CascadeCheckBalance::tolerance),
70  absoluteLimit(G4CascadeCheckBalance::tolerance), initialBaryon(0),
71  finalBaryon(0), initialCharge(0), finalCharge(0), initialStrange(0),
72  finalStrange(0) {}
73 
75  G4double absolute,
76  const char* owner)
77  : G4VCascadeCollider(owner), relativeLimit(relative),
78  absoluteLimit(absolute), initialBaryon(0), finalBaryon(0),
79  initialCharge(0), finalCharge(0), initialStrange(0),
80  finalStrange(0) {}
81 
82 
83 // Pseudo-collision just computes input and output four-vectors
84 
86  G4InuclParticle* target,
87  G4CollisionOutput& output) {
88  if (verboseLevel)
89  G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide"
90  << G4endl;
91 
92  initial *= 0.; // Fast reset; some colliders only have one pointer
93  if (bullet) initial += bullet->getMomentum();
94  if (target) initial += target->getMomentum();
95 
96  // Baryon number, charge and strangeness must be computed "by hand"
97  initialCharge = 0;
98  if (bullet) initialCharge += G4int(bullet->getCharge());
99  if (target) initialCharge += G4int(target->getCharge());
100 
101  G4InuclElementaryParticle* pbullet =
102  dynamic_cast<G4InuclElementaryParticle*>(bullet);
103  G4InuclElementaryParticle* ptarget =
104  dynamic_cast<G4InuclElementaryParticle*>(target);
105 
106  G4InuclNuclei* nbullet = dynamic_cast<G4InuclNuclei*>(bullet);
107  G4InuclNuclei* ntarget = dynamic_cast<G4InuclNuclei*>(target);
108 
109  initialBaryon =
110  ((pbullet ? pbullet->baryon() : nbullet ? nbullet->getA() : 0) +
111  (ptarget ? ptarget->baryon() : ntarget ? ntarget->getA() : 0) );
112 
113  // NOTE: Currently we ignore possibility of hypernucleus target
114  initialStrange = 0;
115  if (pbullet) initialStrange += pbullet->getStrangeness();
116  if (ptarget) initialStrange += ptarget->getStrangeness();
117 
118  // Final state totals are computed for us
119  final = output.getTotalOutputMomentum();
121  finalCharge = output.getTotalCharge();
123 
124  // Report results
125  if (verboseLevel) {
126  G4cout << " initial px " << initial.px() << " py " << initial.py()
127  << " pz " << initial.pz() << " E " << initial.e()
128  << " baryon " << initialBaryon << " charge " << initialCharge
129  << " strange " << initialStrange << G4endl
130  << " final px " << final.px() << " py " << final.py()
131  << " pz " << final.pz() << " E " << final.e()
132  << " baryon " << finalBaryon << " charge " << finalCharge
133  << " strange " << finalStrange << G4endl;
134  }
135 }
136 
137 // For de-excitation, take G4Fragment as initial state
139  G4CollisionOutput& output) {
140  if (verboseLevel)
141  G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<FRAG>)"
142  << G4endl;
143 
144  // Copy initial state directly from fragment (no bullet/target sums)
145  initial = fragment.GetMomentum();
146  initialCharge = G4int(fragment.GetZ());
147  initialBaryon = G4int(fragment.GetA());
148  initialStrange = 0; // No hypernuclei at present
149 
150  // Final state totals are computed for us
151  final = output.getTotalOutputMomentum();
153  finalCharge = output.getTotalCharge();
155 
156  // Report results
157  if (verboseLevel) {
158  G4cout << " initial px " << initial.px() << " py " << initial.py()
159  << " pz " << initial.pz() << " E " << initial.e()
160  << " baryon " << initialBaryon << " charge " << initialCharge
161  << " strange " << initialStrange << G4endl
162  << " final px " << final.px() << " py " << final.py()
163  << " pz " << final.pz() << " E " << final.e()
164  << " baryon " << finalBaryon << " charge " << finalCharge
165  << " strange " << finalStrange << G4endl;
166  }
167 }
168 
169 
170 // Take list of output particles directly (e.g., from G4EPCollider internals)
171 
173  G4InuclParticle* target,
174  const std::vector<G4InuclElementaryParticle>& particles) {
175  if (verboseLevel)
176  G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<vector>)"
177  << G4endl;
178 
179  tempOutput.reset(); // Buffer for processing
180  tempOutput.addOutgoingParticles(particles);
181  collide(bullet, target, tempOutput);
182 }
183 
185  const std::vector<G4InuclElementaryParticle>& particles) {
186  if (verboseLevel)
187  G4cout << " >>> G4CascadeCheckBalance(" << theName
188  << ")::collide(<FRAG>,<vector>)" << G4endl;
189 
190  tempOutput.reset(); // Buffer for processing
191  tempOutput.addOutgoingParticles(particles);
192  collide(target, tempOutput);
193 }
194 
195 
196 // Take list of nuclear fragments directly (e.g., from G4Fissioner internals)
197 
199  const std::vector<G4InuclNuclei>& fragments) {
200  if (verboseLevel)
201  G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<vector>)"
202  << G4endl;
203 
204  tempOutput.reset(); // Buffer for processing
205  tempOutput.addOutgoingNuclei(fragments);
206  collide(target, tempOutput);
207 }
208 
209 
210 // Take list of "cparticles" (e.g., from G4NucleiModel internals)
211 
213  G4InuclParticle* target,
214  const std::vector<G4CascadParticle>& particles) {
215  if (verboseLevel)
216  G4cout << " >>> G4CascadeCheckBalance(" << theName
217  << ")::collide(<cparticles>)" << G4endl;
218 
219  tempOutput.reset(); // Buffer for processing
220  tempOutput.addOutgoingParticles(particles);
221  collide(bullet, target, tempOutput);
222 }
223 
224 
225 // Take lists of both G4InuclEP & G4CP (e.g., from G4IntraNucleiCascader)
226 
228  G4InuclParticle* target,
229  G4CollisionOutput& output,
230  const std::vector<G4CascadParticle>& cparticles) {
231  if (verboseLevel)
232  G4cout << " >>> G4CascadeCheckBalance(" << theName
233  << ")::collide(<EP>,<CP>)" << G4endl;
234 
235  tempOutput.reset(); // Buffer for processing
236  tempOutput.add(output);
237  tempOutput.addOutgoingParticles(cparticles);
238  collide(bullet, target, tempOutput);
239 }
240 
241 
242 // Compare relative and absolute violations to limits, and report
243 
245  G4bool relokay = (std::abs(relativeE()) < relativeLimit);
246  G4bool absokay = (std::abs(deltaE()) < absoluteLimit);
247 
248  if (verboseLevel && !(relokay || absokay)) {
249  G4cerr << theName << ": Energy conservation: relative " << relativeE()
250  << (relokay ? " conserved" : " VIOLATED")
251  << " absolute " << deltaE()
252  << (absokay ? " conserved" : " VIOLATED") << G4endl;
253  } else if (verboseLevel > 1) {
254  G4cout << theName << ": Energy conservation: relative " << relativeE()
255  << " conserved absolute " << deltaE() << " conserved" << G4endl;
256  }
257 
258  return (relokay && absokay);
259 }
260 
262  G4bool relokay = (std::abs(relativeKE()) < relativeLimit);
263  G4bool absokay = (std::abs(deltaKE()) < absoluteLimit);
264 
265  if (verboseLevel && !(relokay || absokay)) {
266  G4cerr << theName << ": Kinetic energy balance: relative "
267  << relativeKE() << (relokay ? " conserved" : " VIOLATED")
268  << " absolute " << deltaKE()
269  << (absokay ? " conserved" : " VIOLATED") << G4endl;
270  } else if (verboseLevel > 1) {
271  G4cout << theName << ": Kinetic energy balance: relative "
272  << relativeKE() << " conserved absolute " << deltaKE()
273  << " conserved" << G4endl;
274  }
275 
276  return (relokay && absokay);
277 }
278 
280  G4bool relokay = (std::abs(relativeP()) < relativeLimit);
281  G4bool absokay = (std::abs(deltaP()) < absoluteLimit);
282 
283  if (verboseLevel && !(relokay || absokay)) {
284  G4cerr << theName << ": Momentum conservation: relative " << relativeP()
285  << (relokay ? " conserved" : " VIOLATED")
286  << " absolute " << deltaP()
287  << (absokay ? " conserved" : " VIOLATED") << G4endl;
288  } else if (verboseLevel > 1) {
289  G4cout << theName << ": Momentum conservation: relative " << relativeP()
290  << " conserved absolute " << deltaP() << " conserved" << G4endl;
291  }
292 
293  return (relokay && absokay);
294 }
295 
297  G4bool bokay = (deltaB() == 0); // Must be perfect!
298 
299  if (verboseLevel && !bokay)
300  G4cerr << theName << ": Baryon number VIOLATED " << deltaB() << G4endl;
301 
302  return bokay;
303 }
304 
306  G4bool qokay = (deltaQ() == 0); // Must be perfect!
307 
308  if (verboseLevel && !qokay)
309  G4cerr << theName << ": Charge conservation VIOLATED " << deltaQ()
310  << G4endl;
311 
312  return qokay;
313 }
314 
316  G4bool sokay = (deltaS() == 0); // Must be perfect!
317 
318  if (verboseLevel && !sokay)
319  G4cerr << theName << ": Strangeness conservation VIOLATED " << deltaS()
320  << G4endl;
321 
322  return sokay;
323 }
void addOutgoingNuclei(const std::vector< G4InuclNuclei > &nuclea)
G4LorentzVector getMomentum() const
G4double GetA() const
Definition: G4Fragment.hh:303
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
static const G4double tolerance
int G4int
Definition: G4Types.hh:78
G4LorentzVector getTotalOutputMomentum() const
void add(const G4CollisionOutput &right)
G4CascadeCheckBalance(const char *owner="G4CascadeCheckBalance")
G4GLOB_DLL std::ostream G4cout
G4int getA() const
bool G4bool
Definition: G4Types.hh:79
G4int getTotalStrangeness() const
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:271
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
G4int getTotalBaryonNumber() const
G4int getTotalCharge() const
G4double getCharge() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double GetZ() const
Definition: G4Fragment.hh:298
G4GLOB_DLL std::ostream G4cerr