Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 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........oooOO0OOooo........oooOO0OOooo.... 35 36 using namespace std; 37 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 39 40 G4DNABornExcitationModel1::G4DNABornExcitationModel1(const G4ParticleDefinition*, 41 const G4String& nam) : 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 openings, sampling of atoms 55 // 4 = entering in methods 56 57 if (verboseLevel > 0) 58 { 59 G4cout << "Born excitation model is constructed " << G4endl; 60 } 61 fParticleChangeForGamma = nullptr; 62 63 // Selection of stationary mode 64 65 statCode = false; 66 } 67 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 69 70 G4DNABornExcitationModel1::~G4DNABornExcitationModel1() 71 { 72 // Cross section 73 74 delete fTableData; 75 } 76 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 78 79 void G4DNABornExcitationModel1::Initialise(const G4ParticleDefinition* particle, 80 const G4DataVector& /*cuts*/) 81 { 82 83 if (verboseLevel > 3) 84 { 85 G4cout << "Calling G4DNABornExcitationModel1::Initialise()" << G4endl; 86 } 87 88 if(fParticleDefinition != nullptr && fParticleDefinition != particle) 89 { 90 G4Exception("G4DNABornExcitationModel1::Initialise","em0001", 91 FatalException,"Model already initialized for another particle type."); 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() == "proton") 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*m; 113 fTableData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor ); 114 fTableData->LoadData(fTableFile); 115 116 if( verboseLevel>0 ) 117 { 118 G4cout << "Born excitation model is initialized " << G4endl 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::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER")); 128 129 if (isInitialised) 130 { return;} 131 fParticleChangeForGamma = GetParticleChangeForGamma(); 132 isInitialised = true; 133 } 134 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 136 137 G4double G4DNABornExcitationModel1::CrossSectionPerVolume(const G4Material* material, 138 const G4ParticleDefinition* particleDefinition, 139 G4double ekin, 140 G4double, 141 G4double) 142 { 143 if (verboseLevel > 3) 144 { 145 G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel1" 146 << G4endl; 147 } 148 149 if(particleDefinition != fParticleDefinition) return 0; 150 151 // Calculate total cross section for model 152 153 G4double sigma=0; 154 155 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()]; 156 157 if (ekin >= fLowEnergy && ekin <= fHighEnergy) 158 { 159 sigma = fTableData->FindValue(ekin); 160 } 161 162 if (verboseLevel > 2) 163 { 164 G4cout << "__________________________________" << G4endl; 165 G4cout << "G4DNABornExcitationModel1 - XS INFO START" << G4endl; 166 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl; 167 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl; 168 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl; 169 G4cout << "G4DNABornExcitationModel1 - XS INFO END" << G4endl; 170 } 171 172 return sigma*waterDensity; 173 } 174 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 176 177 void G4DNABornExcitationModel1::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/, 178 const G4MaterialCutsCouple* /*couple*/, 179 const G4DynamicParticle* aDynamicParticle, 180 G4double, 181 G4double) 182 { 183 184 if (verboseLevel > 3) 185 { 186 G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel1" 187 << G4endl; 188 } 189 190 G4double k = aDynamicParticle->GetKineticEnergy(); 191 192 G4int level = RandomSelect(k); 193 G4double excitationEnergy = waterStructure.ExcitationEnergy(level); 194 G4double newEnergy = k - excitationEnergy; 195 196 if (newEnergy > 0) 197 { 198 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection()); 199 200 if (!statCode) fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy); 201 else fParticleChangeForGamma->SetProposedKineticEnergy(k); 202 203 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy); 204 } 205 206 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack(); 207 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule, 208 level, 209 theIncomingTrack); 210 } 211 212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 213 214 G4double G4DNABornExcitationModel1::GetPartialCrossSection(const G4Material*, 215 G4int level, 216 const G4ParticleDefinition* particle, 217 G4double kineticEnergy) 218 { 219 if (fParticleDefinition != particle) 220 { 221 G4Exception("G4DNABornExcitationModel1::GetPartialCrossSection", 222 "bornParticleType", 223 FatalException, 224 "Model initialized for another particle type."); 225 } 226 227 return fTableData->GetComponent(level)->FindValue(kineticEnergy); 228 } 229 230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 231 232 G4int G4DNABornExcitationModel1::RandomSelect(G4double k) 233 { 234 G4int level = 0; 235 236 auto valuesBuffer = new G4double[fTableData->NumberOfComponents()]; 237 const auto n = (G4int)fTableData->NumberOfComponents(); 238 G4int i(n); 239 G4double value = 0.; 240 241 while (i > 0) 242 { 243 i--; 244 valuesBuffer[i] = fTableData->GetComponent(i)->FindValue(k); 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