Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAEmfietzoglouExcitationModel.cc
Go to the documentation of this file.
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 // Based on the work described in
27 // Rad Res 163, 98-111 (2005)
28 // D. Emfietzoglou, H. Nikjoo
29 //
30 // Authors of the class (2014):
31 // I. Kyriakou (kyriak@cc.uoi.gr)
32 // D. Emfietzoglou (demfietz@cc.uoi.gr)
33 // S. Incerti (incerti@cenbg.in2p3.fr)
34 //
35 
37 #include "G4SystemOfUnits.hh"
38 #include "G4DNAChemistryManager.hh"
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
43 using namespace std;
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 
48  const G4String& nam)
49  :G4VEmModel(nam),isInitialised(false)
50 {
51 
52  // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
53  fpMolWaterDensity = 0;
54 
55  verboseLevel= 0;
56  // Verbosity scale:
57  // 0 = nothing
58  // 1 = warning for energy non-conservation
59  // 2 = details of energy budget
60  // 3 = calculation of cross sections, file openings, sampling of atoms
61  // 4 = entering in methods
62 
63  if( verboseLevel>0 )
64  {
65  G4cout << "Emfietzoglou excitation model is constructed " << G4endl;
66  }
68  SetLowEnergyLimit(8. * eV);
69  SetHighEnergyLimit(10. * keV);
70 
71  // Selection of stationary mode
72 
73  statCode = false;
74 }
75 
76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77 
79 {
80  // Cross section
81 
82  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
83  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
84  {
85  G4DNACrossSectionDataSet* table = pos->second;
86  delete table;
87  }
88 
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
94  const G4DataVector& /*cuts*/)
95 {
96 
97  if (verboseLevel > 3)
98  G4cout << "Calling G4DNAEmfietzoglouExcitationModel::Initialise()" << G4endl;
99 
100  G4String fileElectron("dna/sigma_excitation_e_emfietzoglou");
101 
103 
105 
106  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
107 
108  // *** ELECTRON
109 
110  electron = electronDef->GetParticleName();
111 
112  tableFile[electron] = fileElectron;
113 
114  lowEnergyLimit[electron] = 8. * eV;
115  highEnergyLimit[electron] = 10. * keV;
116 
117  // Cross section
118 
120  tableE->LoadData(fileElectron);
121 
122  tableData[electron] = tableE;
123 
124  //
125 
126  if (particle==electronDef)
127  {
128  SetLowEnergyLimit(lowEnergyLimit[electron]);
129  SetHighEnergyLimit(highEnergyLimit[electron]);
130  }
131 
132  if( verboseLevel>0 )
133  {
134  G4cout << "Emfietzoglou excitation model is initialized " << G4endl
135  << "Energy range: "
136  << LowEnergyLimit() / eV << " eV - "
137  << HighEnergyLimit() / keV << " keV for "
138  << particle->GetParticleName()
139  << G4endl;
140  }
141 
142  // Initialize water density pointer
144 
145  if (isInitialised) { return; }
147  isInitialised = true;
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 
153  const G4ParticleDefinition* particleDefinition,
154  G4double ekin,
155  G4double,
156  G4double)
157 {
158  if (verboseLevel > 3)
159  G4cout << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouExcitationModel" << G4endl;
160 
161  if (particleDefinition != G4Electron::ElectronDefinition()) return 0;
162 
163  // Calculate total cross section for model
164 
165  G4double lowLim = 0;
166  G4double highLim = 0;
167  G4double sigma=0;
168 
169  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
170 
171  if(waterDensity!= 0.0)
172  {
173  const G4String& particleName = particleDefinition->GetParticleName();
174 
175  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
176  pos1 = lowEnergyLimit.find(particleName);
177  if (pos1 != lowEnergyLimit.end())
178  {
179  lowLim = pos1->second;
180  }
181 
182  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
183  pos2 = highEnergyLimit.find(particleName);
184  if (pos2 != highEnergyLimit.end())
185  {
186  highLim = pos2->second;
187  }
188 
189  if (ekin >= lowLim && ekin < highLim)
190  {
191  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
192  pos = tableData.find(particleName);
193 
194  if (pos != tableData.end())
195  {
196  G4DNACrossSectionDataSet* table = pos->second;
197  if (table != 0)
198  {
199  sigma = table->FindValue(ekin);
200  }
201  }
202  else
203  {
204  G4Exception("G4DNAEmfietzoglouExcitationModel::CrossSectionPerVolume","em0002",
205  FatalException,"Model not applicable to particle type.");
206  }
207  }
208 
209  if (verboseLevel > 2)
210  {
211  G4cout << "__________________________________" << G4endl;
212  G4cout << "G4DNAEmfietzoglouExcitationModel - XS INFO START" << G4endl;
213  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
214  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
215  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
216  // G4cout << " Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
217  G4cout << "G4DNAEmfietzoglouExcitationModel - XS INFO END" << G4endl;
218  }
219 
220  } // if (waterMaterial)
221 
222  return sigma*waterDensity;
223 }
224 
225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
226 
227 void G4DNAEmfietzoglouExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
228  const G4MaterialCutsCouple* /*couple*/,
229  const G4DynamicParticle* aDynamicParticle,
230  G4double,
231  G4double)
232 {
233 
234  if (verboseLevel > 3)
235  G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouExcitationModel" << G4endl;
236 
237  G4double k = aDynamicParticle->GetKineticEnergy();
238 
239  const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
240 
241  G4int level = RandomSelect(k,particleName);
242  G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
243  G4double newEnergy = k - excitationEnergy;
244 
245  if (newEnergy > 0)
246  {
248 
249  if (!statCode) fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
251 
253  }
254 
255  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
257  level,
258  theIncomingTrack);
259 }
260 
261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
262 
263 G4int G4DNAEmfietzoglouExcitationModel::RandomSelect(G4double k, const G4String& particle)
264 {
265 
266  G4int level = 0;
267 
268  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
269  pos = tableData.find(particle);
270 
271  if (pos != tableData.end())
272  {
273  G4DNACrossSectionDataSet* table = pos->second;
274 
275  if (table != 0)
276  {
277  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
278  const size_t n(table->NumberOfComponents());
279  size_t i(n);
280  G4double value = 0.;
281 
282  //Check reading of initial xs file
283  //G4cout << table->GetComponent(0)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
284  //G4cout << table->GetComponent(1)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
285  //G4cout << table->GetComponent(2)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
286  //G4cout << table->GetComponent(3)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
287  //G4cout << table->GetComponent(4)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
288  //G4cout << table->GetComponent(5)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
289  //G4cout << table->GetComponent(6)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
290  //abort();
291 
292  while (i>0)
293  {
294  i--;
295  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
296  value += valuesBuffer[i];
297  }
298 
299  value *= G4UniformRand();
300 
301  i = n;
302 
303  while (i > 0)
304  {
305  i--;
306 
307  if (valuesBuffer[i] > value)
308  {
309  delete[] valuesBuffer;
310  return i;
311  }
312  value -= valuesBuffer[i];
313  }
314 
315  if (valuesBuffer) delete[] valuesBuffer;
316 
317  }
318  }
319  else
320  {
321  G4Exception("G4DNAEmfietzoglouExcitationModel::RandomSelect","em0002",
322  FatalException,"Model not applicable to particle type.");
323  }
324  return level;
325 }
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:602
size_t GetIndex() const
Definition: G4Material.hh:262
virtual G4bool LoadData(const G4String &argFileName)
G4ParticleDefinition * GetDefinition() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &=*(new G4DataVector()))
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
string material
Definition: eplot.py:19
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
static constexpr double m
Definition: G4SIunits.hh:129
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4ThreeVector & GetMomentumDirection() const
static constexpr double cm
Definition: G4SIunits.hh:119
static constexpr double eV
Definition: G4SIunits.hh:215
const G4int n
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual size_t NumberOfComponents(void) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4DNAChemistryManager * Instance()
static G4DNAMolecularMaterial * Instance()
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
const G4Track * GetCurrentTrack() const
G4DNAEmfietzoglouExcitationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAEmfietzoglouExcitationModel")
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:739
static constexpr double keV
Definition: G4SIunits.hh:216
static const G4double pos
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132