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