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 "G4DNABornExcitationModel2.hh" 29 #include "G4SystemOfUnits.hh" 30 #include "G4DNAChemistryManager.hh" 31 #include "G4DNAMolecularMaterial.hh" 32 #include "G4PhysicsTable.hh" 33 #include "G4PhysicsVector.hh" 34 #include "G4UnitsTable.hh" 35 #include <map> 36 37 //....oooOO0OOooo........oooOO0OOooo........oo 38 39 using namespace std; 40 41 //....oooOO0OOooo........oooOO0OOooo........oo 42 43 G4DNABornExcitationModel2::G4DNABornExcitation 44 45 G4VEmModel(nam) 46 { 47 fpMolWaterDensity = nullptr; 48 fHighEnergy = 0; 49 fLowEnergy = 0; 50 fParticleDefinition = nullptr; 51 52 verboseLevel = 0; 53 // Verbosity scale: 54 // 0 = nothing 55 // 1 = warning for energy non-conservation 56 // 2 = details of energy budget 57 // 3 = calculation of cross sections, file o 58 // 4 = entering in methods 59 60 if (verboseLevel > 0) 61 { 62 G4cout << "Born excitation model is constr 63 } 64 fParticleChangeForGamma = nullptr; 65 fLastBinCallForFinalXS = 0; 66 fTotalXS = nullptr; 67 fTableData = nullptr; 68 69 // Selection of stationary mode 70 71 statCode = false; 72 } 73 74 //....oooOO0OOooo........oooOO0OOooo........oo 75 76 G4DNABornExcitationModel2::~G4DNABornExcitatio 77 { 78 // Cross section 79 80 delete fTableData; 81 } 82 83 //....oooOO0OOooo........oooOO0OOooo........oo 84 85 void G4DNABornExcitationModel2::Initialise(con 86 con 87 { 88 89 if (verboseLevel > 3) 90 { 91 G4cout << "Calling G4DNABornExcitationMode 92 } 93 94 if(fParticleDefinition != nullptr && fPartic 95 { 96 G4Exception("G4DNABornExcitationModel2::In 97 FatalException,"Model already initiali 98 } 99 100 fParticleDefinition = particle; 101 102 std::ostringstream fullFileName; 103 const char* path = G4FindDataDir("G4LEDATA") 104 105 if(G4String(path).empty()) 106 { 107 G4Exception("G4DNABornExcitationModel2::In 108 FatalException, "G4LEDATA not defined 109 } 110 111 fullFileName << path; 112 113 if(particle->GetParticleName() == "e-") 114 { 115 fullFileName << "/dna/bornExcitation-e.dat 116 fLowEnergy = 9*eV; 117 fHighEnergy = 1*MeV; 118 } 119 else if(particle->GetParticleName() == "prot 120 { 121 fullFileName << "/dna/bornExcitation-p.dat 122 fLowEnergy = 500. * keV; 123 fHighEnergy = 100. * MeV; 124 } 125 126 SetLowEnergyLimit(fLowEnergy); 127 SetHighEnergyLimit(fHighEnergy); 128 129 //G4double scaleFactor = (1.e-22 / 3.343) * 130 131 fTableData = new G4PhysicsTable(); 132 fTableData->RetrievePhysicsTable(fullFileNam 133 /* 134 for(std::size_t level = 0; level<fTableData- 135 { 136 //(*fTableData)(level)->ScaleVector(1,scal 137 } 138 */ 139 std::size_t finalBin_i = 2000; 140 G4double E_min = fLowEnergy; 141 G4double E_max = fHighEnergy; 142 fTotalXS = new G4PhysicsLogVector(E_min, E_m 143 G4double energy; 144 G4double finalXS; 145 146 for(std::size_t energy_i = 0; energy_i < fin 147 { 148 energy = fTotalXS->Energy(energy_i); 149 finalXS = 0; 150 151 for(std::size_t level = 0; level<fTableDat 152 { 153 finalXS += (*fTableData)(level)->Value(e 154 } 155 fTotalXS->PutValue(energy_i, finalXS); 156 //G4cout << "energy = " << energy << " " < 157 // << " " << energy_i << " " << fina 158 } 159 160 // for(energy = LowEnergyLimit() ; energy 161 // { 162 // G4cout << "energy = " << energy << " 163 // } 164 165 if( verboseLevel>0 ) 166 { 167 G4cout << "Born excitation model is initia 168 << "Energy range: " 169 << LowEnergyLimit() / eV << " eV - " 170 << HighEnergyLimit() / keV << " keV for " 171 << particle->GetParticleName() 172 << G4endl; 173 } 174 175 // Initialize water density pointer 176 fpMolWaterDensity = G4DNAMolecularMaterial:: 177 178 if (isInitialised) 179 { return;} 180 fParticleChangeForGamma = GetParticleChangeF 181 isInitialised = true; 182 } 183 184 //....oooOO0OOooo........oooOO0OOooo........oo 185 186 G4double G4DNABornExcitationModel2::CrossSecti 187 188 189 190 191 { 192 if (verboseLevel > 3) 193 { 194 G4cout << "Calling CrossSectionPerVolume() 195 << G4endl; 196 } 197 198 if(particleDefinition != fParticleDefinition 199 200 // Calculate total cross section for model 201 202 G4double sigma=0; 203 204 G4double waterDensity = (*fpMolWaterDensity) 205 206 if (ekin >= fLowEnergy && ekin <= fHighEnerg 207 { 208 sigma = fTotalXS->Value(ekin, fLastBinCall 209 210 // for(std::size_t i = 0; i < 5; ++i) 211 // sigma += (*fTableData)[i]->Value(ekin); 212 213 if(sigma == 0) 214 { 215 G4cerr << "PROBLEM SIGMA = 0 at " << G4B 216 } 217 } 218 219 if (verboseLevel > 2) 220 { 221 G4cout << "_______________________________ 222 G4cout << "G4DNABornExcitationModel2 - XS 223 G4cout << "Kinetic energy(eV)=" << ekin/eV 224 G4cout << "Cross section per water molecul 225 G4cout << "Cross section per water molecul 226 G4cout << "G4DNABornExcitationModel2 - XS 227 } 228 229 return sigma*waterDensity; 230 } 231 232 //....oooOO0OOooo........oooOO0OOooo........oo 233 234 void G4DNABornExcitationModel2::SampleSecondar 235 236 237 238 239 { 240 241 if (verboseLevel > 3) 242 { 243 G4cout << "Calling SampleSecondaries() of 244 << G4endl; 245 } 246 247 G4double k = aDynamicParticle->GetKineticEne 248 249 G4int level = RandomSelect(k); 250 G4double excitationEnergy = waterStructure.E 251 G4double newEnergy = k - excitationEnergy; 252 253 if (newEnergy > 0) 254 { 255 fParticleChangeForGamma->ProposeMomentumDi 256 257 if (!statCode) fParticleChangeForGamma->Se 258 else fParticleChangeForGamma->SetProposedK 259 260 fParticleChangeForGamma->ProposeLocalEnerg 261 } 262 263 const G4Track * theIncomingTrack = fParticle 264 G4DNAChemistryManager::Instance()->CreateWat 265 level, 266 theIncomingTrack); 267 } 268 269 //....oooOO0OOooo........oooOO0OOooo........oo 270 271 G4double G4DNABornExcitationModel2::GetPartial 272 273 274 275 { 276 if (fParticleDefinition != particle) 277 { 278 G4Exception("G4DNABornExcitationModel2::Ge 279 "bornParticleType", 280 FatalException, 281 "Model initialized for another 282 } 283 284 return (*fTableData)(level)->Value(kineticEn 285 } 286 287 //....oooOO0OOooo........oooOO0OOooo........oo 288 289 G4int G4DNABornExcitationModel2::RandomSelect( 290 { 291 const std::size_t n(fTableData->size()); 292 std::size_t i(n); 293 294 G4double value = fTotalXS->Value(k, fLastBin 295 296 value *= G4UniformRand(); 297 i = n; 298 299 G4double partialXS; 300 301 while (i > 0) 302 { 303 i--; 304 305 partialXS = (*fTableData)(i)->Value(k); 306 if (partialXS > value) 307 { 308 return (G4int)i; 309 } 310 value -= partialXS; 311 } 312 313 return 0; 314 } 315 316