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