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> << 29 << 30 #include "G4QMDNucleus.hh" 28 #include "G4QMDNucleus.hh" 31 #include "G4Pow.hh" << 32 #include "G4SystemOfUnits.hh" << 33 #include "G4Proton.hh" 29 #include "G4Proton.hh" 34 #include "G4Neutron.hh" 30 #include "G4Neutron.hh" 35 #include "G4NucleiProperties.hh" 31 #include "G4NucleiProperties.hh" 36 #include "G4HadronicException.hh" << 32 >> 33 #include <numeric> 37 34 38 G4QMDNucleus::G4QMDNucleus() 35 G4QMDNucleus::G4QMDNucleus() 39 { 36 { 40 G4QMDParameters* parameters = G4QMDParamete 37 G4QMDParameters* parameters = G4QMDParameters::GetInstance(); 41 hbc = parameters->Get_hbc(); 38 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 } 39 } 47 40 48 41 49 42 50 //G4QMDNucleus::~G4QMDNucleus() << 43 G4QMDNucleus::~G4QMDNucleus() 51 //{ << 44 { 52 // ; << 45 ; 53 //} << 46 } 54 47 55 48 56 G4LorentzVector G4QMDNucleus::Get4Momentum() 49 G4LorentzVector G4QMDNucleus::Get4Momentum() 57 { 50 { 58 G4LorentzVector p( 0 ); 51 G4LorentzVector p( 0 ); 59 std::vector< G4QMDParticipant* >::iterator 52 std::vector< G4QMDParticipant* >::iterator it; 60 for ( it = participants.begin() ; it != par 53 for ( it = participants.begin() ; it != participants.end() ; it++ ) 61 p += (*it)->Get4Momentum(); 54 p += (*it)->Get4Momentum(); 62 55 63 return p; 56 return p; 64 } 57 } 65 58 66 59 67 60 68 G4int G4QMDNucleus::GetMassNumber() 61 G4int G4QMDNucleus::GetMassNumber() 69 { 62 { 70 63 71 G4int A = 0; 64 G4int A = 0; 72 std::vector< G4QMDParticipant* >::iterator 65 std::vector< G4QMDParticipant* >::iterator it; 73 for ( it = participants.begin() ; it != par 66 for ( it = participants.begin() ; it != participants.end() ; it++ ) 74 { 67 { 75 if ( (*it)->GetDefinition() == G4Proton: 68 if ( (*it)->GetDefinition() == G4Proton::Proton() 76 || (*it)->GetDefinition() == G4Neutron 69 || (*it)->GetDefinition() == G4Neutron::Neutron() ) 77 A++; 70 A++; 78 } 71 } 79 72 80 if ( A == 0 ) { << 81 throw G4HadronicException(__FILE__, __LI << 82 } << 83 << 84 return A; 73 return A; 85 } 74 } 86 75 87 76 88 77 89 G4int G4QMDNucleus::GetAtomicNumber() 78 G4int G4QMDNucleus::GetAtomicNumber() 90 { 79 { 91 G4int Z = 0; 80 G4int Z = 0; 92 std::vector< G4QMDParticipant* >::iterator 81 std::vector< G4QMDParticipant* >::iterator it; 93 for ( it = participants.begin() ; it != par 82 for ( it = participants.begin() ; it != participants.end() ; it++ ) 94 { 83 { 95 if ( (*it)->GetDefinition() == G4Proton: 84 if ( (*it)->GetDefinition() == G4Proton::Proton() ) 96 Z++; 85 Z++; 97 } 86 } 98 return Z; 87 return Z; 99 } 88 } 100 89 101 90 102 91 103 G4double G4QMDNucleus::GetNuclearMass() 92 G4double G4QMDNucleus::GetNuclearMass() 104 { 93 { 105 94 106 G4double mass = G4NucleiProperties::GetNucl 95 G4double mass = G4NucleiProperties::GetNuclearMass( GetMassNumber() , GetAtomicNumber() ); 107 96 108 if ( mass == 0.0 ) 97 if ( mass == 0.0 ) 109 { 98 { 110 99 111 G4int Z = GetAtomicNumber(); 100 G4int Z = GetAtomicNumber(); 112 G4int A = GetMassNumber(); 101 G4int A = GetMassNumber(); 113 G4int N = A - Z; 102 G4int N = A - Z; 114 103 115 // Weizsacker-Bethe 104 // Weizsacker-Bethe 116 105 117 G4double Av = 16*MeV; 106 G4double Av = 16*MeV; 118 G4double As = 17*MeV; 107 G4double As = 17*MeV; 119 G4double Ac = 0.7*MeV; 108 G4double Ac = 0.7*MeV; 120 G4double Asym = 23*MeV; 109 G4double Asym = 23*MeV; 121 110 122 G4double BE = Av * A 111 G4double BE = Av * A 123 - As * G4Pow::GetInstance()- << 112 - As * std::pow ( G4double ( A ) , 2.0/3.0 ) 124 - Ac * Z*Z/G4Pow::GetInstanc << 113 - Ac * Z*Z/std::pow ( G4double ( A ) , 1.0/3.0 ) 125 - Asym * ( N - Z )* ( N - Z 114 - Asym * ( N - Z )* ( N - Z ) / A; 126 115 127 mass = Z * G4Proton::Proton()->GetPDGMas 116 mass = Z * G4Proton::Proton()->GetPDGMass() 128 + N * G4Neutron::Neutron()->GetPDGM 117 + N * G4Neutron::Neutron()->GetPDGMass() 129 - BE; 118 - BE; 130 119 131 } 120 } 132 121 133 return mass; 122 return mass; 134 } 123 } 135 124 136 125 137 126 138 void G4QMDNucleus::CalEnergyAndAngularMomentum 127 void G4QMDNucleus::CalEnergyAndAngularMomentumInCM() 139 { 128 { 140 129 141 //G4cout << "CalEnergyAndAngularMomentumInC 130 //G4cout << "CalEnergyAndAngularMomentumInCM " << this->GetAtomicNumber() << " " << GetMassNumber() << G4endl; 142 131 143 G4double gamma = Get4Momentum().gamma(); 132 G4double gamma = Get4Momentum().gamma(); 144 G4ThreeVector beta = Get4Momentum().v()/ Ge 133 G4ThreeVector beta = Get4Momentum().v()/ Get4Momentum().e(); 145 134 146 G4ThreeVector pcm0( 0.0 ) ; 135 G4ThreeVector pcm0( 0.0 ) ; 147 136 148 G4int n = GetTotalNumberOfParticipant(); 137 G4int n = GetTotalNumberOfParticipant(); 149 pcm.resize( n ); 138 pcm.resize( n ); 150 139 151 for ( G4int i= 0; i < n ; i++ ) 140 for ( G4int i= 0; i < n ; i++ ) 152 { 141 { 153 G4ThreeVector p_i = GetParticipant( i )- 142 G4ThreeVector p_i = GetParticipant( i )->GetMomentum(); 154 143 155 G4double trans = gamma / ( gamma + 1.0 ) 144 G4double trans = gamma / ( gamma + 1.0 ) * p_i * beta; 156 pcm[i] = p_i - trans*beta; 145 pcm[i] = p_i - trans*beta; 157 146 158 pcm0 += pcm[i]; 147 pcm0 += pcm[i]; 159 } 148 } 160 149 161 pcm0 = pcm0 / double ( n ); 150 pcm0 = pcm0 / double ( n ); 162 151 163 //G4cout << "pcm0 " << pcm0 << G4endl; 152 //G4cout << "pcm0 " << pcm0 << G4endl; 164 153 165 for ( G4int i= 0; i < n ; i++ ) 154 for ( G4int i= 0; i < n ; i++ ) 166 { 155 { 167 pcm[i] += -pcm0; 156 pcm[i] += -pcm0; 168 //G4cout << "pcm " << i << " " << pcm[i] 157 //G4cout << "pcm " << i << " " << pcm[i] << G4endl; 169 } 158 } 170 159 171 160 172 G4double tmass = 0; 161 G4double tmass = 0; 173 G4ThreeVector rcm0( 0.0 ) ; 162 G4ThreeVector rcm0( 0.0 ) ; 174 rcm.resize( n ); 163 rcm.resize( n ); 175 es.resize( n ); 164 es.resize( n ); 176 165 177 for ( G4int i= 0; i < n ; i++ ) 166 for ( G4int i= 0; i < n ; i++ ) 178 { 167 { 179 G4ThreeVector ri = GetParticipant( i )-> 168 G4ThreeVector ri = GetParticipant( i )->GetPosition(); 180 G4double trans = gamma / ( gamma + 1.0 ) 169 G4double trans = gamma / ( gamma + 1.0 ) * ri * beta; 181 170 182 es[i] = std::sqrt ( G4Pow::GetInstance() << 171 es[i] = std::sqrt ( std::pow ( GetParticipant( i )->GetMass() , 2 ) + pcm[i]*pcm[i] ); 183 172 184 rcm[i] = ri + trans*beta; 173 rcm[i] = ri + trans*beta; 185 174 186 rcm0 += rcm[i]*es[i]; 175 rcm0 += rcm[i]*es[i]; 187 176 188 tmass += es[i]; 177 tmass += es[i]; 189 } 178 } 190 179 191 rcm0 = rcm0/tmass; 180 rcm0 = rcm0/tmass; 192 181 193 for ( G4int i= 0; i < n ; i++ ) 182 for ( G4int i= 0; i < n ; i++ ) 194 { 183 { 195 rcm[i] += -rcm0; 184 rcm[i] += -rcm0; 196 //G4cout << "rcm " << i << " " << rcm[i] 185 //G4cout << "rcm " << i << " " << rcm[i] << G4endl; 197 } 186 } 198 187 199 // Angular momentum << 188 // Angluar momentum 200 189 201 G4ThreeVector rl ( 0.0 ); 190 G4ThreeVector rl ( 0.0 ); 202 for ( G4int i= 0; i < n ; i++ ) 191 for ( G4int i= 0; i < n ; i++ ) 203 { 192 { 204 rl += rcm[i].cross ( pcm[i] ); 193 rl += rcm[i].cross ( pcm[i] ); 205 } 194 } 206 195 207 // DHW: move hbc outside of sqrt to get correc << 196 jj = int ( std::sqrt ( rl*rl / hbc ) + 0.5 ); 208 // jj = int ( std::sqrt ( rl*rl / hbc ) + 0.5 << 209 197 210 jj = int (std::sqrt(rl*rl)/hbc + 0.5); << 211 198 212 // kinetic energy per nucleon in CM 199 // kinetic energy per nucleon in CM 213 200 214 G4double totalMass = 0.0; 201 G4double totalMass = 0.0; 215 for ( G4int i= 0; i < n ; i++ ) 202 for ( G4int i= 0; i < n ; i++ ) 216 { 203 { 217 // following two lines are equivalent 204 // following two lines are equivalent 218 //totalMass += GetParticipant( i )->GetD 205 //totalMass += GetParticipant( i )->GetDefinition()->GetPDGMass()/GeV; 219 totalMass += GetParticipant( i )->GetMas 206 totalMass += GetParticipant( i )->GetMass(); 220 } 207 } 221 208 222 //G4double kineticEnergyPerNucleon = ( std: << 209 kineticEnergyPerNucleon = ( std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass )/n; 223 210 224 // Total (not per nucleion ) Binding Energy 211 // Total (not per nucleion ) Binding Energy 225 G4double bindingEnergy = ( std::accumulate << 212 bindingEnergy = ( std::accumulate ( es.begin() , es.end() , 0.0 ) -totalMass ) + potentialEnergy; 226 213 227 //G4cout << "KineticEnergyPerNucleon in GeV 214 //G4cout << "KineticEnergyPerNucleon in GeV " << kineticEnergyPerNucleon << G4endl; 228 //G4cout << "KineticEnergySum in GeV " << s 215 //G4cout << "KineticEnergySum in GeV " << std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass << G4endl; 229 //G4cout << "PotentialEnergy in GeV " << po 216 //G4cout << "PotentialEnergy in GeV " << potentialEnergy << G4endl; 230 //G4cout << "BindingEnergy in GeV " << bind 217 //G4cout << "BindingEnergy in GeV " << bindingEnergy << G4endl; 231 //G4cout << "G4BindingEnergy in GeV " << G4 218 //G4cout << "G4BindingEnergy in GeV " << G4NucleiProperties::GetBindingEnergy( GetAtomicNumber() , GetMassNumber() )/GeV << G4endl; 232 219 233 excitationEnergy = bindingEnergy + G4Nuclei 220 excitationEnergy = bindingEnergy + G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() )/GeV; 234 //G4cout << "excitationEnergy in GeV " << e 221 //G4cout << "excitationEnergy in GeV " << excitationEnergy << G4endl; 235 if ( excitationEnergy < 0 ) excitationEnerg 222 if ( excitationEnergy < 0 ) excitationEnergy = 0.0; 236 223 237 } 224 } 238 225