1 | #ifndef G4CASCADE_CHECK_BALANCE_HH |
---|
2 | #define G4CASCADE_CHECK_BALANCE_HH |
---|
3 | // |
---|
4 | // ******************************************************************** |
---|
5 | // * License and Disclaimer * |
---|
6 | // * * |
---|
7 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
8 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
9 | // * conditions of the Geant4 Software License, included in the file * |
---|
10 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
11 | // * include a list of copyright holders. * |
---|
12 | // * * |
---|
13 | // * Neither the authors of this software system, nor their employing * |
---|
14 | // * institutes,nor the agencies providing financial support for this * |
---|
15 | // * work make any representation or warranty, express or implied, * |
---|
16 | // * regarding this software system or assume any liability for its * |
---|
17 | // * use. Please see the license in the file LICENSE and URL above * |
---|
18 | // * for the full disclaimer and the limitation of liability. * |
---|
19 | // * * |
---|
20 | // * This code implementation is the result of the scientific and * |
---|
21 | // * technical work of the GEANT4 collaboration. * |
---|
22 | // * By using, copying, modifying or distributing the software (or * |
---|
23 | // * any work based on the software) you agree to acknowledge its * |
---|
24 | // * use in resulting scientific publications, and indicate your * |
---|
25 | // * acceptance of all terms of the Geant4 Software license. * |
---|
26 | // ******************************************************************** |
---|
27 | // |
---|
28 | // $Id: G4CascadeCheckBalance.hh,v 1.11 2010/09/24 06:26:06 mkelsey Exp $ |
---|
29 | // Geant4 tag: $Name: hadr-casc-V09-03-85 $ |
---|
30 | // |
---|
31 | // Verify and report four-momentum conservation for collision output; uses |
---|
32 | // same interface as collision generators. |
---|
33 | // |
---|
34 | // 20100624 M. Kelsey -- Add baryon conservation check and kinetic energy |
---|
35 | // 20100628 M. Kelsey -- Add interface to take list of particles directly |
---|
36 | // 20100711 M. Kelsey -- Add name of parent Collider for reporting messages, |
---|
37 | // allow changing parent name, add interface for nuclear fragments |
---|
38 | // 20100713 M. Kelsey -- Add (publicly adjustable) tolerance for zeroes |
---|
39 | // 20100715 M. Kelsey -- FPE! Need to check initial values before computing |
---|
40 | // relative error. |
---|
41 | // 20100715 M. Kelsey -- Add G4CascadParticle interface for G4NucleiModel; |
---|
42 | // do momentum check on direction, not just magnitude. Move |
---|
43 | // temporary G4CollisionOutput buffer here, for thread-safety |
---|
44 | // 20100909 M. Kelsey -- Add interface to get four-vector difference, and |
---|
45 | // to supply both kinds of particle lists (G4IntraNucleiCascader) |
---|
46 | // 20100923 M. Kelsey -- Baryon and charge deltas should have been integer |
---|
47 | |
---|
48 | #include "G4VCascadeCollider.hh" |
---|
49 | #include "globals.hh" |
---|
50 | #include "G4CollisionOutput.hh" |
---|
51 | #include "G4LorentzVector.hh" |
---|
52 | #include <cmath> |
---|
53 | #include <vector> |
---|
54 | |
---|
55 | class G4CascadParticle; |
---|
56 | class G4InuclElementaryParticle; |
---|
57 | class G4InuclNuclei; |
---|
58 | class G4InuclParticle; |
---|
59 | |
---|
60 | class G4CascadeCheckBalance : public G4VCascadeCollider { |
---|
61 | public: |
---|
62 | static const G4double tolerance; // Don't do floating zero! |
---|
63 | |
---|
64 | G4CascadeCheckBalance(G4double relative, G4double absolute, |
---|
65 | const char* owner="G4CascadeCheckBalance"); |
---|
66 | virtual ~G4CascadeCheckBalance() {}; |
---|
67 | |
---|
68 | void setOwner(const char* owner) { setName(owner); } |
---|
69 | |
---|
70 | void collide(G4InuclParticle* bullet, G4InuclParticle* target, |
---|
71 | G4CollisionOutput& output); |
---|
72 | |
---|
73 | // This is for use with G4EPCollider internal checks |
---|
74 | void collide(G4InuclParticle* bullet, G4InuclParticle* target, |
---|
75 | const std::vector<G4InuclElementaryParticle>& particles); |
---|
76 | |
---|
77 | // This is for use with G4Fissioner internal checks |
---|
78 | void collide(G4InuclParticle* bullet, G4InuclParticle* target, |
---|
79 | const std::vector<G4InuclNuclei>& fragments); |
---|
80 | |
---|
81 | // This is for use with G4NucleiModel internal checks |
---|
82 | void collide(G4InuclParticle* bullet, G4InuclParticle* target, |
---|
83 | const std::vector<G4CascadParticle>& particles); |
---|
84 | |
---|
85 | // This is for use with G4IntraNucleiCascader |
---|
86 | void collide(G4InuclParticle* bullet, G4InuclParticle* target, |
---|
87 | const std::vector<G4InuclElementaryParticle>& particles, |
---|
88 | const std::vector<G4CascadParticle>& cparticles); |
---|
89 | |
---|
90 | // Checks on conservation laws (kinematics, baryon number, charge) |
---|
91 | G4bool energyOkay() const; |
---|
92 | G4bool ekinOkay() const; |
---|
93 | G4bool momentumOkay() const; |
---|
94 | G4bool baryonOkay() const; |
---|
95 | G4bool chargeOkay() const; |
---|
96 | |
---|
97 | // Global check, used by G4CascadeInterface validation loop |
---|
98 | G4bool okay() const { return (energyOkay() && momentumOkay() && |
---|
99 | baryonOkay() && chargeOkay()); } |
---|
100 | |
---|
101 | // Calculations of conserved quantities from initial and final state |
---|
102 | // FIXME: Relative comparisons don't work for zero! |
---|
103 | G4double deltaE() const { return (final.e() - initial.e()); } |
---|
104 | G4double relativeE() const { |
---|
105 | return ( (std::abs(deltaE())<tolerance) ? 0. : |
---|
106 | (initial.e()<tolerance) ? 1. : deltaE()/initial.e() ); |
---|
107 | } |
---|
108 | |
---|
109 | G4double deltaKE() const { return (ekin(final) - ekin(initial)); } |
---|
110 | G4double relativeKE() const { |
---|
111 | return ( (std::abs(deltaKE())<tolerance) ? 0. : |
---|
112 | (ekin(initial)<tolerance) ? 1. : deltaKE()/ekin(initial) ); |
---|
113 | } |
---|
114 | |
---|
115 | G4double deltaP() const { return deltaLV().rho(); } |
---|
116 | G4double relativeP() const { |
---|
117 | return ( (std::abs(deltaP())<tolerance) ? 0. : |
---|
118 | (initial.rho()<tolerance) ? 1. : deltaP()/initial.rho() ); |
---|
119 | } |
---|
120 | |
---|
121 | G4LorentzVector deltaLV() const { return final - initial; } |
---|
122 | |
---|
123 | // Baryon number and charge are discrete; no bounds and no "relative" scale |
---|
124 | G4int deltaB() const { return (finalBaryon - initialBaryon); } |
---|
125 | G4int deltaQ() const { return (finalCharge - initialCharge); } |
---|
126 | |
---|
127 | protected: |
---|
128 | // Utility function for kinetic energy |
---|
129 | G4double ekin(const G4LorentzVector& p) const { return (p.e() - p.m()); } |
---|
130 | |
---|
131 | private: |
---|
132 | G4double relativeLimit; |
---|
133 | G4double absoluteLimit; |
---|
134 | |
---|
135 | G4LorentzVector initial; // Four-vectors for computing violations |
---|
136 | G4LorentzVector final; |
---|
137 | |
---|
138 | G4int initialBaryon; // Total baryon number |
---|
139 | G4int finalBaryon; |
---|
140 | |
---|
141 | G4int initialCharge; // Total charge |
---|
142 | G4int finalCharge; |
---|
143 | |
---|
144 | G4CollisionOutput tempOutput; // Buffer for direct-list interfaces |
---|
145 | }; |
---|
146 | |
---|
147 | #endif /* G4CASCADE_CHECK_BALANCE_HH */ |
---|