Geant4  10.02.p02
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 // 20140930 Change name from "const char*" to "const G4String"
53 // 20150626 Fix logical condition for report relative or absolute violations.
54 
55 #include "G4CascadeCheckBalance.hh"
56 #include "globals.hh"
57 #include "G4CascadParticle.hh"
59 #include "G4InuclNuclei.hh"
60 #include "G4InuclParticle.hh"
61 #include "G4CollisionOutput.hh"
62 #include "G4LorentzVector.hh"
63 #include <vector>
64 
65 
66 // Constructor sets acceptance limits
67 
68 const G4double G4CascadeCheckBalance::tolerance = 1e-6; // How small is zero?
69 
71  : G4VCascadeCollider(owner), relativeLimit(G4CascadeCheckBalance::tolerance),
72  absoluteLimit(G4CascadeCheckBalance::tolerance), initialBaryon(0),
73  finalBaryon(0), initialCharge(0), finalCharge(0), initialStrange(0),
74  finalStrange(0) {}
75 
77  G4double absolute,
78  const G4String& owner)
79  : G4VCascadeCollider(owner), relativeLimit(relative),
80  absoluteLimit(absolute), initialBaryon(0), finalBaryon(0),
81  initialCharge(0), finalCharge(0), initialStrange(0),
82  finalStrange(0) {}
83 
84 
85 // Pseudo-collision just computes input and output four-vectors
86 
88  G4InuclParticle* target,
89  G4CollisionOutput& output) {
90  if (verboseLevel)
91  G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide"
92  << G4endl;
93 
94  initial *= 0.; // Fast reset; some colliders only have one pointer
95  if (bullet) initial += bullet->getMomentum();
96  if (target) initial += target->getMomentum();
97 
98  // Baryon number, charge and strangeness must be computed "by hand"
99  initialCharge = 0;
100  if (bullet) initialCharge += G4int(bullet->getCharge());
101  if (target) initialCharge += G4int(target->getCharge());
102 
103  G4InuclElementaryParticle* pbullet =
104  dynamic_cast<G4InuclElementaryParticle*>(bullet);
105  G4InuclElementaryParticle* ptarget =
106  dynamic_cast<G4InuclElementaryParticle*>(target);
107 
108  G4InuclNuclei* nbullet = dynamic_cast<G4InuclNuclei*>(bullet);
109  G4InuclNuclei* ntarget = dynamic_cast<G4InuclNuclei*>(target);
110 
111  initialBaryon =
112  ((pbullet ? pbullet->baryon() : nbullet ? nbullet->getA() : 0) +
113  (ptarget ? ptarget->baryon() : ntarget ? ntarget->getA() : 0) );
114 
115  // NOTE: Currently we ignore possibility of hypernucleus target
116  initialStrange = 0;
117  if (pbullet) initialStrange += pbullet->getStrangeness();
118  if (ptarget) initialStrange += ptarget->getStrangeness();
119 
120  // Final state totals are computed for us
121  final = output.getTotalOutputMomentum();
123  finalCharge = output.getTotalCharge();
125 
126  // Report results
127  if (verboseLevel) {
128  G4cout << " initial px " << initial.px() << " py " << initial.py()
129  << " pz " << initial.pz() << " E " << initial.e()
130  << " baryon " << initialBaryon << " charge " << initialCharge
131  << " strange " << initialStrange << G4endl
132  << " final px " << final.px() << " py " << final.py()
133  << " pz " << final.pz() << " E " << final.e()
134  << " baryon " << finalBaryon << " charge " << finalCharge
135  << " strange " << finalStrange << G4endl;
136  }
137 }
138 
139 // For de-excitation, take G4Fragment as initial state
141  G4CollisionOutput& output) {
142  if (verboseLevel)
143  G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<FRAG>)"
144  << G4endl;
145 
146  // Copy initial state directly from fragment (no bullet/target sums)
147  initial = fragment.GetMomentum();
148  initialCharge = fragment.GetZ_asInt();
149  initialBaryon = fragment.GetA_asInt();
150  initialStrange = 0; // No hypernuclei at present
151 
152  // Final state totals are computed for us
153  final = output.getTotalOutputMomentum();
155  finalCharge = output.getTotalCharge();
157 
158  // Report results
159  if (verboseLevel) {
160  G4cout << " initial px " << initial.px() << " py " << initial.py()
161  << " pz " << initial.pz() << " E " << initial.e()
162  << " baryon " << initialBaryon << " charge " << initialCharge
163  << " strange " << initialStrange << G4endl
164  << " final px " << final.px() << " py " << final.py()
165  << " pz " << final.pz() << " E " << final.e()
166  << " baryon " << finalBaryon << " charge " << finalCharge
167  << " strange " << finalStrange << G4endl;
168  }
169 }
170 
171 
172 // Take list of output particles directly (e.g., from G4EPCollider internals)
173 
175  G4InuclParticle* target,
176  const std::vector<G4InuclElementaryParticle>& particles) {
177  if (verboseLevel)
178  G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<vector>)"
179  << G4endl;
180 
181  tempOutput.reset(); // Buffer for processing
182  tempOutput.addOutgoingParticles(particles);
183  collide(bullet, target, tempOutput);
184 }
185 
187  const std::vector<G4InuclElementaryParticle>& particles) {
188  if (verboseLevel)
189  G4cout << " >>> G4CascadeCheckBalance(" << theName
190  << ")::collide(<FRAG>,<vector>)" << G4endl;
191 
192  tempOutput.reset(); // Buffer for processing
193  tempOutput.addOutgoingParticles(particles);
194  collide(target, tempOutput);
195 }
196 
197 
198 // Take list of nuclear fragments directly (e.g., from G4Fissioner internals)
199 
201  const std::vector<G4InuclNuclei>& fragments) {
202  if (verboseLevel)
203  G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<vector>)"
204  << G4endl;
205 
206  tempOutput.reset(); // Buffer for processing
207  tempOutput.addOutgoingNuclei(fragments);
208  collide(target, tempOutput);
209 }
210 
211 
212 // Take list of "cparticles" (e.g., from G4NucleiModel internals)
213 
215  G4InuclParticle* target,
216  const std::vector<G4CascadParticle>& particles) {
217  if (verboseLevel)
218  G4cout << " >>> G4CascadeCheckBalance(" << theName
219  << ")::collide(<cparticles>)" << G4endl;
220 
221  tempOutput.reset(); // Buffer for processing
222  tempOutput.addOutgoingParticles(particles);
223  collide(bullet, target, tempOutput);
224 }
225 
226 
227 // Take lists of both G4InuclEP & G4CP (e.g., from G4IntraNucleiCascader)
228 
230  G4InuclParticle* target,
231  G4CollisionOutput& output,
232  const std::vector<G4CascadParticle>& cparticles) {
233  if (verboseLevel)
234  G4cout << " >>> G4CascadeCheckBalance(" << theName
235  << ")::collide(<EP>,<CP>)" << G4endl;
236 
237  tempOutput.reset(); // Buffer for processing
238  tempOutput.add(output);
239  tempOutput.addOutgoingParticles(cparticles);
240  collide(bullet, target, tempOutput);
241 }
242 
243 
244 // Compare relative and absolute violations to limits, and report
245 
247  G4bool relokay = (std::abs(relativeE()) < relativeLimit);
248  G4bool absokay = (std::abs(deltaE()) < absoluteLimit);
249 
250  if (verboseLevel && (!relokay || !absokay)) {
251  G4cerr << theName << ": Energy conservation: relative " << relativeE()
252  << (relokay ? " conserved" : " VIOLATED")
253  << " absolute " << deltaE()
254  << (absokay ? " conserved" : " VIOLATED") << G4endl;
255  } else if (verboseLevel > 1) {
256  G4cout << theName << ": Energy conservation: relative " << relativeE()
257  << " conserved absolute " << deltaE() << " conserved" << G4endl;
258  }
259 
260  return (relokay && absokay);
261 }
262 
264  G4bool relokay = (std::abs(relativeKE()) < relativeLimit);
265  G4bool absokay = (std::abs(deltaKE()) < absoluteLimit);
266 
267  if (verboseLevel && (!relokay || !absokay)) {
268  G4cerr << theName << ": Kinetic energy balance: relative "
269  << relativeKE() << (relokay ? " conserved" : " VIOLATED")
270  << " absolute " << deltaKE()
271  << (absokay ? " conserved" : " VIOLATED") << G4endl;
272  } else if (verboseLevel > 1) {
273  G4cout << theName << ": Kinetic energy balance: relative "
274  << relativeKE() << " conserved absolute " << deltaKE()
275  << " conserved" << G4endl;
276  }
277 
278  return (relokay && absokay);
279 }
280 
282  G4bool relokay = (std::abs(relativeP()) < relativeLimit);
283  G4bool absokay = (std::abs(deltaP()) < absoluteLimit);
284 
285  if (verboseLevel && (!relokay || !absokay)) {
286  G4cerr << theName << ": Momentum conservation: relative " << relativeP()
287  << (relokay ? " conserved" : " VIOLATED")
288  << " absolute " << deltaP()
289  << (absokay ? " conserved" : " VIOLATED") << G4endl;
290  } else if (verboseLevel > 1) {
291  G4cout << theName << ": Momentum conservation: relative " << relativeP()
292  << " conserved absolute " << deltaP() << " conserved" << G4endl;
293  }
294 
295  return (relokay && absokay);
296 }
297 
299  G4bool bokay = (deltaB() == 0); // Must be perfect!
300 
301  if (verboseLevel && !bokay)
302  G4cerr << theName << ": Baryon number VIOLATED " << deltaB() << G4endl;
303 
304  return bokay;
305 }
306 
308  G4bool qokay = (deltaQ() == 0); // Must be perfect!
309 
310  if (verboseLevel && !qokay)
311  G4cerr << theName << ": Charge conservation VIOLATED " << deltaQ()
312  << G4endl;
313 
314  return qokay;
315 }
316 
318  G4bool sokay = (deltaS() == 0); // Must be perfect!
319 
320  if (verboseLevel && !sokay)
321  G4cerr << theName << ": Strangeness conservation VIOLATED " << deltaS()
322  << G4endl;
323 
324  return sokay;
325 }
void addOutgoingNuclei(const std::vector< G4InuclNuclei > &nuclea)
G4LorentzVector getMomentum() const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
static const G4float tolerance
static const G4double tolerance
int G4int
Definition: G4Types.hh:78
G4LorentzVector getTotalOutputMomentum() const
void add(const G4CollisionOutput &right)
G4GLOB_DLL std::ostream G4cout
G4int GetA_asInt() const
Definition: G4Fragment.hh:256
G4int getA() const
bool G4bool
Definition: G4Types.hh:79
G4int getTotalStrangeness() const
G4CascadeCheckBalance(const G4String &owner="G4CascadeCheckBalance")
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:289
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
G4int getTotalBaryonNumber() const
G4int GetZ_asInt() const
Definition: G4Fragment.hh:261
G4int getTotalCharge() const
G4double getCharge() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4GLOB_DLL std::ostream G4cerr