| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // 081024 G4NucleiPropertiesTable:: to G4Nucle
27 //
28 #include <numeric>
29
30 #include "G4QMDNucleus.hh"
31 #include "G4Pow.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4Proton.hh"
34 #include "G4Neutron.hh"
35 #include "G4NucleiProperties.hh"
36 #include "G4HadronicException.hh"
37
38 G4QMDNucleus::G4QMDNucleus()
39 {
40 G4QMDParameters* parameters = G4QMDParamete
41 hbc = parameters->Get_hbc();
42
43 jj = 0; // will be calcualted in CalEnergyA
44 potentialEnergy = 0.0; // will be set throu
45 excitationEnergy = 0.0;
46 }
47
48
49
50 //G4QMDNucleus::~G4QMDNucleus()
51 //{
52 // ;
53 //}
54
55
56 G4LorentzVector G4QMDNucleus::Get4Momentum()
57 {
58 G4LorentzVector p( 0 );
59 std::vector< G4QMDParticipant* >::iterator
60 for ( it = participants.begin() ; it != par
61 p += (*it)->Get4Momentum();
62
63 return p;
64 }
65
66
67
68 G4int G4QMDNucleus::GetMassNumber()
69 {
70
71 G4int A = 0;
72 std::vector< G4QMDParticipant* >::iterator
73 for ( it = participants.begin() ; it != par
74 {
75 if ( (*it)->GetDefinition() == G4Proton:
76 || (*it)->GetDefinition() == G4Neutron
77 A++;
78 }
79
80 if ( A == 0 ) {
81 throw G4HadronicException(__FILE__, __LI
82 }
83
84 return A;
85 }
86
87
88
89 G4int G4QMDNucleus::GetAtomicNumber()
90 {
91 G4int Z = 0;
92 std::vector< G4QMDParticipant* >::iterator
93 for ( it = participants.begin() ; it != par
94 {
95 if ( (*it)->GetDefinition() == G4Proton:
96 Z++;
97 }
98 return Z;
99 }
100
101
102
103 G4double G4QMDNucleus::GetNuclearMass()
104 {
105
106 G4double mass = G4NucleiProperties::GetNucl
107
108 if ( mass == 0.0 )
109 {
110
111 G4int Z = GetAtomicNumber();
112 G4int A = GetMassNumber();
113 G4int N = A - Z;
114
115 // Weizsacker-Bethe
116
117 G4double Av = 16*MeV;
118 G4double As = 17*MeV;
119 G4double Ac = 0.7*MeV;
120 G4double Asym = 23*MeV;
121
122 G4double BE = Av * A
123 - As * G4Pow::GetInstance()-
124 - Ac * Z*Z/G4Pow::GetInstanc
125 - Asym * ( N - Z )* ( N - Z
126
127 mass = Z * G4Proton::Proton()->GetPDGMas
128 + N * G4Neutron::Neutron()->GetPDGM
129 - BE;
130
131 }
132
133 return mass;
134 }
135
136
137
138 void G4QMDNucleus::CalEnergyAndAngularMomentum
139 {
140
141 //G4cout << "CalEnergyAndAngularMomentumInC
142
143 G4double gamma = Get4Momentum().gamma();
144 G4ThreeVector beta = Get4Momentum().v()/ Ge
145
146 G4ThreeVector pcm0( 0.0 ) ;
147
148 G4int n = GetTotalNumberOfParticipant();
149 pcm.resize( n );
150
151 for ( G4int i= 0; i < n ; i++ )
152 {
153 G4ThreeVector p_i = GetParticipant( i )-
154
155 G4double trans = gamma / ( gamma + 1.0 )
156 pcm[i] = p_i - trans*beta;
157
158 pcm0 += pcm[i];
159 }
160
161 pcm0 = pcm0 / double ( n );
162
163 //G4cout << "pcm0 " << pcm0 << G4endl;
164
165 for ( G4int i= 0; i < n ; i++ )
166 {
167 pcm[i] += -pcm0;
168 //G4cout << "pcm " << i << " " << pcm[i]
169 }
170
171
172 G4double tmass = 0;
173 G4ThreeVector rcm0( 0.0 ) ;
174 rcm.resize( n );
175 es.resize( n );
176
177 for ( G4int i= 0; i < n ; i++ )
178 {
179 G4ThreeVector ri = GetParticipant( i )->
180 G4double trans = gamma / ( gamma + 1.0 )
181
182 es[i] = std::sqrt ( G4Pow::GetInstance()
183
184 rcm[i] = ri + trans*beta;
185
186 rcm0 += rcm[i]*es[i];
187
188 tmass += es[i];
189 }
190
191 rcm0 = rcm0/tmass;
192
193 for ( G4int i= 0; i < n ; i++ )
194 {
195 rcm[i] += -rcm0;
196 //G4cout << "rcm " << i << " " << rcm[i]
197 }
198
199 // Angular momentum
200
201 G4ThreeVector rl ( 0.0 );
202 for ( G4int i= 0; i < n ; i++ )
203 {
204 rl += rcm[i].cross ( pcm[i] );
205 }
206
207 // DHW: move hbc outside of sqrt to get correc
208 // jj = int ( std::sqrt ( rl*rl / hbc ) + 0.5
209
210 jj = int (std::sqrt(rl*rl)/hbc + 0.5);
211
212 // kinetic energy per nucleon in CM
213
214 G4double totalMass = 0.0;
215 for ( G4int i= 0; i < n ; i++ )
216 {
217 // following two lines are equivalent
218 //totalMass += GetParticipant( i )->GetD
219 totalMass += GetParticipant( i )->GetMas
220 }
221
222 //G4double kineticEnergyPerNucleon = ( std:
223
224 // Total (not per nucleion ) Binding Energy
225 G4double bindingEnergy = ( std::accumulate
226
227 //G4cout << "KineticEnergyPerNucleon in GeV
228 //G4cout << "KineticEnergySum in GeV " << s
229 //G4cout << "PotentialEnergy in GeV " << po
230 //G4cout << "BindingEnergy in GeV " << bind
231 //G4cout << "G4BindingEnergy in GeV " << G4
232
233 excitationEnergy = bindingEnergy + G4Nuclei
234 //G4cout << "excitationEnergy in GeV " << e
235 if ( excitationEnergy < 0 ) excitationEnerg
236
237 }
238