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