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 // 27 28 #include "G4DNABornExcitationModel1.hh" 29 #include "G4SystemOfUnits.hh" 30 #include "G4DNAChemistryManager.hh" 31 #include "G4DNAMolecularMaterial.hh" 32 #include <map> 33 34 //....oooOO0OOooo........oooOO0OOooo........oo 35 36 using namespace std; 37 38 //....oooOO0OOooo........oooOO0OOooo........oo 39 40 G4DNABornExcitationModel1::G4DNABornExcitation 41 42 G4VEmModel(nam) 43 { 44 fpMolWaterDensity = nullptr; 45 fHighEnergy = 0; 46 fLowEnergy = 0; 47 fParticleDefinition = nullptr; 48 49 verboseLevel = 0; 50 // Verbosity scale: 51 // 0 = nothing 52 // 1 = warning for energy non-conservation 53 // 2 = details of energy budget 54 // 3 = calculation of cross sections, file o 55 // 4 = entering in methods 56 57 if (verboseLevel > 0) 58 { 59 G4cout << "Born excitation model is constr 60 } 61 fParticleChangeForGamma = nullptr; 62 63 // Selection of stationary mode 64 65 statCode = false; 66 } 67 68 //....oooOO0OOooo........oooOO0OOooo........oo 69 70 G4DNABornExcitationModel1::~G4DNABornExcitatio 71 { 72 // Cross section 73 74 delete fTableData; 75 } 76 77 //....oooOO0OOooo........oooOO0OOooo........oo 78 79 void G4DNABornExcitationModel1::Initialise(con 80 con 81 { 82 83 if (verboseLevel > 3) 84 { 85 G4cout << "Calling G4DNABornExcitationMode 86 } 87 88 if(fParticleDefinition != nullptr && fPartic 89 { 90 G4Exception("G4DNABornExcitationModel1::In 91 FatalException,"Model already initiali 92 } 93 94 fParticleDefinition = particle; 95 96 if(particle->GetParticleName() == "e-") 97 { 98 fTableFile = "dna/sigma_excitation_e_born" 99 fLowEnergy = 9*eV; 100 fHighEnergy = 1*MeV; 101 } 102 else if(particle->GetParticleName() == "prot 103 { 104 fTableFile = "dna/sigma_excitation_p_born" 105 fLowEnergy = 500. * keV; 106 fHighEnergy = 100. * MeV; 107 } 108 109 SetLowEnergyLimit(fLowEnergy); 110 SetHighEnergyLimit(fHighEnergy); 111 112 G4double scaleFactor = (1.e-22 / 3.343) * m* 113 fTableData = new G4DNACrossSectionDataSet(ne 114 fTableData->LoadData(fTableFile); 115 116 if( verboseLevel>0 ) 117 { 118 G4cout << "Born excitation model is initia 119 << "Energy range: " 120 << LowEnergyLimit() / eV << " eV - " 121 << HighEnergyLimit() / keV << " keV for " 122 << particle->GetParticleName() 123 << G4endl; 124 } 125 126 // Initialize water density pointer 127 fpMolWaterDensity = G4DNAMolecularMaterial:: 128 129 if (isInitialised) 130 { return;} 131 fParticleChangeForGamma = GetParticleChangeF 132 isInitialised = true; 133 } 134 135 //....oooOO0OOooo........oooOO0OOooo........oo 136 137 G4double G4DNABornExcitationModel1::CrossSecti 138 139 140 141 142 { 143 if (verboseLevel > 3) 144 { 145 G4cout << "Calling CrossSectionPerVolume() 146 << G4endl; 147 } 148 149 if(particleDefinition != fParticleDefinition 150 151 // Calculate total cross section for model 152 153 G4double sigma=0; 154 155 G4double waterDensity = (*fpMolWaterDensity) 156 157 if (ekin >= fLowEnergy && ekin <= fHighEnerg 158 { 159 sigma = fTableData->FindValue(ekin); 160 } 161 162 if (verboseLevel > 2) 163 { 164 G4cout << "_______________________________ 165 G4cout << "G4DNABornExcitationModel1 - XS 166 G4cout << "Kinetic energy(eV)=" << ekin/eV 167 G4cout << "Cross section per water molecul 168 G4cout << "Cross section per water molecul 169 G4cout << "G4DNABornExcitationModel1 - XS 170 } 171 172 return sigma*waterDensity; 173 } 174 175 //....oooOO0OOooo........oooOO0OOooo........oo 176 177 void G4DNABornExcitationModel1::SampleSecondar 178 179 180 181 182 { 183 184 if (verboseLevel > 3) 185 { 186 G4cout << "Calling SampleSecondaries() of 187 << G4endl; 188 } 189 190 G4double k = aDynamicParticle->GetKineticEne 191 192 G4int level = RandomSelect(k); 193 G4double excitationEnergy = waterStructure.E 194 G4double newEnergy = k - excitationEnergy; 195 196 if (newEnergy > 0) 197 { 198 fParticleChangeForGamma->ProposeMomentumDi 199 200 if (!statCode) fParticleChangeForGamma->Se 201 else fParticleChangeForGamma->SetProposedK 202 203 fParticleChangeForGamma->ProposeLocalEnerg 204 } 205 206 const G4Track * theIncomingTrack = fParticle 207 G4DNAChemistryManager::Instance()->CreateWat 208 level, 209 theIncomingTrack); 210 } 211 212 //....oooOO0OOooo........oooOO0OOooo........oo 213 214 G4double G4DNABornExcitationModel1::GetPartial 215 216 217 218 { 219 if (fParticleDefinition != particle) 220 { 221 G4Exception("G4DNABornExcitationModel1::Ge 222 "bornParticleType", 223 FatalException, 224 "Model initialized for another 225 } 226 227 return fTableData->GetComponent(level)->Find 228 } 229 230 //....oooOO0OOooo........oooOO0OOooo........oo 231 232 G4int G4DNABornExcitationModel1::RandomSelect( 233 { 234 G4int level = 0; 235 236 auto valuesBuffer = new G4double[fTableData 237 const auto n = (G4int)fTableData->NumberOfC 238 G4int i(n); 239 G4double value = 0.; 240 241 while (i > 0) 242 { 243 i--; 244 valuesBuffer[i] = fTableData->GetComponent 245 value += valuesBuffer[i]; 246 } 247 248 value *= G4UniformRand(); 249 i = n; 250 251 while (i > 0) 252 { 253 i--; 254 255 if (valuesBuffer[i] > value) 256 { 257 delete[] valuesBuffer; 258 return i; 259 } 260 value -= valuesBuffer[i]; 261 } 262 263 264 delete[] valuesBuffer; 265 266 return level; 267 } 268 269