Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/models/src/G4DNABornExcitationModel2.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 "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........oooOO0OOooo........oooOO0OOooo....
 38 
 39 using namespace std;
 40 
 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 42 
 43 G4DNABornExcitationModel2::G4DNABornExcitationModel2(const G4ParticleDefinition*,
 44                                                      const G4String& nam) :
 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 openings, sampling of atoms
 58   // 4 = entering in methods
 59 
 60   if (verboseLevel > 0)
 61   {
 62     G4cout << "Born excitation model is constructed " << G4endl;
 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........oooOO0OOooo........oooOO0OOooo....
 75 
 76 G4DNABornExcitationModel2::~G4DNABornExcitationModel2()
 77 {
 78   // Cross section
 79   
 80     delete fTableData;
 81 }
 82 
 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 84 
 85 void G4DNABornExcitationModel2::Initialise(const G4ParticleDefinition* particle,
 86                                            const G4DataVector& /*cuts*/)
 87 {
 88 
 89   if (verboseLevel > 3)
 90   {
 91     G4cout << "Calling G4DNABornExcitationModel2::Initialise()" << G4endl;
 92   }
 93 
 94   if(fParticleDefinition != nullptr && fParticleDefinition != particle)
 95   {
 96     G4Exception("G4DNABornExcitationModel2::Initialise","em0001",
 97         FatalException,"Model already initialized for another particle type.");
 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::Initialise","G4LEDATA-CHECK",
108         FatalException, "G4LEDATA not defined in environment variables");
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() == "proton")
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) * m*m;
130 
131   fTableData = new G4PhysicsTable();
132   fTableData->RetrievePhysicsTable(fullFileName.str().c_str(), true);
133   /*
134   for(std::size_t level = 0; level<fTableData->size(); ++level)
135   {
136     //(*fTableData)(level)->ScaleVector(1,scaleFactor);
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_max, finalBin_i, true);
143   G4double energy;
144   G4double finalXS;
145 
146   for(std::size_t energy_i = 0; energy_i < finalBin_i; ++energy_i)
147   {
148     energy = fTotalXS->Energy(energy_i);
149     finalXS = 0;
150 
151     for(std::size_t level = 0; level<fTableData->size(); ++level)
152     {
153       finalXS += (*fTableData)(level)->Value(energy);
154     }
155     fTotalXS->PutValue(energy_i, finalXS);
156     //G4cout << "energy = " << energy << " " << fTotalXS->Value(energy)
157     //       << " " << energy_i << " " << finalXS << G4endl;
158   }
159 
160   //    for(energy = LowEnergyLimit() ; energy < HighEnergyLimit() ; energy += 1*pow(10,log10(energy)))
161   //    {
162   //      G4cout << "energy = " << energy << " " << fTotalXS->Value(energy) << G4endl;
163   //    }
164 
165   if( verboseLevel>0 )
166   {
167     G4cout << "Born excitation model is initialized " << G4endl
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::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
177 
178   if (isInitialised)
179   { return;}
180   fParticleChangeForGamma = GetParticleChangeForGamma();
181   isInitialised = true;
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 
186 G4double G4DNABornExcitationModel2::CrossSectionPerVolume(const G4Material* material,
187                                                           const G4ParticleDefinition* particleDefinition,
188                                                           G4double ekin,
189                                                           G4double,
190                                                           G4double)
191 {
192   if (verboseLevel > 3)
193   {
194     G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel2"
195         << G4endl;
196   }
197 
198   if(particleDefinition != fParticleDefinition) return 0;
199 
200   // Calculate total cross section for model
201 
202   G4double sigma=0;
203 
204   G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
205 
206   if (ekin >= fLowEnergy && ekin <= fHighEnergy)
207   {
208     sigma = fTotalXS->Value(ekin, fLastBinCallForFinalXS);
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 " << G4BestUnit(ekin, "Energy")<< G4endl;
216     }
217   }
218 
219   if (verboseLevel > 2)
220   {
221     G4cout << "__________________________________" << G4endl;
222     G4cout << "G4DNABornExcitationModel2 - XS INFO START" << G4endl;
223     G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
224     G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
225     G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
226     G4cout << "G4DNABornExcitationModel2 - XS INFO END" << G4endl;
227   }
228 
229   return sigma*waterDensity;
230 }
231 
232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233 
234 void G4DNABornExcitationModel2::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
235                                                   const G4MaterialCutsCouple* /*couple*/,
236                                                   const G4DynamicParticle* aDynamicParticle,
237                                                   G4double,
238                                                   G4double)
239 {
240 
241   if (verboseLevel > 3)
242   {
243     G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel2"
244            << G4endl;
245   }
246 
247   G4double k = aDynamicParticle->GetKineticEnergy();
248 
249   G4int level = RandomSelect(k);
250   G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
251   G4double newEnergy = k - excitationEnergy;
252 
253   if (newEnergy > 0)
254   {
255     fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
256 
257     if (!statCode) fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
258     else fParticleChangeForGamma->SetProposedKineticEnergy(k);
259     
260     fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
261   }
262 
263   const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
264   G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
265       level,
266       theIncomingTrack);
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 
271 G4double G4DNABornExcitationModel2::GetPartialCrossSection(const G4Material*,
272                                                            G4int level,
273                                                            const G4ParticleDefinition* particle,
274                                                            G4double kineticEnergy)
275 {
276   if (fParticleDefinition != particle)
277   {
278     G4Exception("G4DNABornExcitationModel2::GetPartialCrossSection",
279                 "bornParticleType",
280                 FatalException,
281                 "Model initialized for another particle type.");
282   }
283 
284   return (*fTableData)(level)->Value(kineticEnergy);
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 
289 G4int G4DNABornExcitationModel2::RandomSelect(G4double k)
290 {
291   const std::size_t n(fTableData->size());
292   std::size_t i(n);
293 
294   G4double value = fTotalXS->Value(k, fLastBinCallForFinalXS);
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