Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNASancheExcitationModel.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 // $Id: G4DNASancheExcitationModel.cc 98733 2016-08-09 10:51:58Z gcosmo $
27 //
28 
29 // Created by Z. Francis
30 
32 #include "G4SystemOfUnits.hh"
34 
35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
36 
37 using namespace std;
38 
39 //#define SANCHE_VERBOSE
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 
44  const G4String& nam) :
45  G4VEmModel(nam), isInitialised(false)
46 {
47  fpWaterDensity = 0;
48 
50  SetHighEnergyLimit(100 * eV);
51  nLevels = 9;
52 
53  verboseLevel = 0;
54  // Verbosity scale:
55  // 0 = nothing
56  // 1 = warning for energy non-conservation
57  // 2 = details of energy budget
58  // 3 = calculation of cross sections, file openings, sampling of atoms
59  // 4 = entering in methods
60 
61 #ifdef SANCHE_VERBOSE
62  if (verboseLevel > 0)
63  {
64  G4cout << "Sanche Excitation model is constructed "
65  << G4endl
66  << "Energy range: "
67  << LowEnergyLimit() / eV << " eV - "
68  << HighEnergyLimit() / eV << " eV"
69  << G4endl;
70  }
71 #endif
72 
74  fpWaterDensity = 0;
75 
76  // Selection of stationary mode
77 
78  statCode = false;
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
82 
84 {
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 
89 void
91 Initialise(const G4ParticleDefinition* /*particle*/,
92  const G4DataVector& /*cuts*/)
93 {
94 
95 #ifdef SANCHE_VERBOSE
96  if (verboseLevel > 3)
97  {
98  G4cout << "Calling G4DNASancheExcitationModel::Initialise()"
99  << G4endl;
100  }
101 #endif
102 
103  // Energy limits
104 
105  if (LowEnergyLimit() < 2.*eV)
106  {
107  G4Exception("*** WARNING : the G4DNASancheExcitationModel class is not "
108  "validated below 2 eV !",
109  "", JustWarning, "");
110  }
111 
112  if (HighEnergyLimit() > 100.*eV)
113  {
114  G4cout << "G4DNASancheExcitationModel: high energy limit decreased from " <<
115  HighEnergyLimit()/eV << " eV to " << 100. << " eV" << G4endl;
116  SetHighEnergyLimit(100.*eV);
117  }
118 
119  //
120 #ifdef SANCHE_VERBOSE
121  if (verboseLevel > 0)
122  {
123  G4cout << "Sanche Excitation model is initialized " << G4endl
124  << "Energy range: "
125  << LowEnergyLimit() / eV << " eV - "
126  << HighEnergyLimit() / eV << " eV"
127  << G4endl;
128  }
129 #endif
130 
131  // Initialize water density pointer
132  fpWaterDensity = G4DNAMolecularMaterial::Instance()->
133  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
134 
135  if (isInitialised) {return;} // RETURNS HERE
136 
138  isInitialised = true;
139 
140  char *path = getenv("G4LEDATA");
141  std::ostringstream eFullFileName;
142  eFullFileName << path << "/dna/sigma_excitationvib_e_sanche.dat";
143  std::ifstream input(eFullFileName.str().c_str());
144 
145  if (!input)
146  {
147  G4Exception("G4DNASancheExcitationModel::Initialise","em0003",
148  FatalException,"Missing data file:/dna/sigma_excitationvib_e_sanche.dat");
149  }
150 
151  // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
152  // Added clear for MT
153  tdummyVec.clear();
154  //
155 
156  double t;
157  double xs;
158 
159  while(!input.eof())
160  {
161  input>>t;
162  tdummyVec.push_back(t);
163 
164  fEnergyLevelXS.push_back(std::vector<double>());
165  fEnergyTotalXS.push_back(0);
166  std::vector<double>& levelXS = fEnergyLevelXS.back();
167  levelXS.reserve(9);
168 
169 // G4cout<<t;
170 
171  for(size_t i = 0 ; i < 9 ;++i)
172  {
173  input>>xs;
174  levelXS.push_back(xs);
175  fEnergyTotalXS.back() += xs;
176 // G4cout <<" " << levelXS[i];
177  }
178 
179 // G4cout << G4endl;
180  }
181 }
182 
183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
184 
186  const G4ParticleDefinition*
187 #ifdef SANCHE_VERBOSE
188  particleDefinition
189 #endif
190  ,
191  G4double ekin,
192  G4double,
193  G4double)
194 {
195 #ifdef SANCHE_VERBOSE
196  if (verboseLevel > 3)
197  {
198  G4cout << "Calling CrossSectionPerVolume() of G4DNASancheExcitationModel"
199  << G4endl;
200  }
201 #endif
202 
203  // Calculate total cross section for model
204 
205  G4double sigma = 0.;
206 
207  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
208 
209  if(waterDensity!= 0.0)
210  {
211  if (ekin >= LowEnergyLimit() && ekin < HighEnergyLimit())
212  // for now this is necessary
213  {
214  //sigma = Sum(ekin);
215  sigma = TotalCrossSection(ekin);
216  }
217 
218 #ifdef SANCHE_VERBOSE
219  if (verboseLevel > 2)
220  {
221  G4cout << "__________________________________" << G4endl;
222  G4cout << "=== G4DNASancheExcitationModel - 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 << "=== G4DNASancheExcitationModel - XS INFO END" << G4endl;
227  }
228 #endif
229  } // if water
230 
231  return sigma*2.*waterDensity;
232  // see papers for factor 2 description
233 
234 }
235 
236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
237 
240  const G4MaterialCutsCouple*,
241  const G4DynamicParticle* aDynamicElectron,
242  G4double,
243  G4double)
244 {
245 #ifdef SANCHE_VERBOSE
246  if (verboseLevel > 3)
247  {
248  G4cout << "Calling SampleSecondaries() of G4DNASancheExcitationModel"
249  << G4endl;
250  }
251 #endif
252 
253  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
254  G4int level = RandomSelect(electronEnergy0);
255  G4double excitationEnergy = VibrationEnergy(level); // levels go from 0 to 8
256  G4double newEnergy = electronEnergy0 - excitationEnergy;
257 
258  /*
259  if (electronEnergy0 < highEnergyLimit)
260  {
261  if (newEnergy >= lowEnergyLimit)
262  {
263  fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
264  fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
265  fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
266  }
267 
268  else
269  {
270  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
271  fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
272  }
273  }
274  */
275 
276  if (electronEnergy0 < HighEnergyLimit() && newEnergy>0.)
277  {
278 
279  if (!statCode)
280  {
284  }
285 
286  else
287  {
291  }
292 
293  }
294 
295  //
296 }
297 
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
299 
301  G4int level)
302 {
303  std::vector<double>::iterator t2 = std::upper_bound(tdummyVec.begin(),
304  tdummyVec.end(), t / eV);
305  std::vector<double>::iterator t1 = t2 - 1;
306 
307  size_t i1 = t1 - tdummyVec.begin();
308  size_t i2 = t2 - tdummyVec.begin();
309 
310  double sigma = LinInterpolate((*t1), (*t2),
311  t / eV,
312  fEnergyLevelXS[i1][level],
313  fEnergyLevelXS[i2][level]);
314 
315  static const double conv_factor = 1e-16 * cm * cm;
316 
317  sigma *= conv_factor;
318  if (sigma == 0.) sigma = 1e-30;
319  return (sigma);
320 }
321 
322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
323 
325 {
326  std::vector<double>::iterator t2 = std::upper_bound(tdummyVec.begin(),
327  tdummyVec.end(), t / eV);
328  std::vector<double>::iterator t1 = t2 - 1;
329 
330  size_t i1 = t1 - tdummyVec.begin();
331  size_t i2 = t2 - tdummyVec.begin();
332 
333  double sigma = LinInterpolate((*t1), (*t2),
334  t / eV,
335  fEnergyTotalXS[i1],
336  fEnergyTotalXS[i2]);
337 
338  static const double conv_factor = 1e-16 * cm * cm;
339 
340  sigma *= conv_factor;
341  if (sigma == 0.) sigma = 1e-30;
342  return (sigma);
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
346 
347 G4double G4DNASancheExcitationModel::VibrationEnergy(G4int level)
348 {
349  static G4double energies[9] = { 0.01, 0.024, 0.061, 0.092, 0.204, 0.417, 0.460,
350  0.500, 0.835 };
351  return (energies[level] * eV);
352 }
353 
354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
355 
356 G4int G4DNASancheExcitationModel::RandomSelect(G4double k)
357 {
358 
359  // Level Selection Counting can be done here !
360 
361  G4int i = nLevels;
362  G4double value = 0.;
363  std::deque<double> values;
364 
365  while (i > 0)
366  {
367  i--;
368  G4double partial = PartialCrossSection(k, i);
369  values.push_front(partial);
370  value += partial;
371  }
372 
373  value *= G4UniformRand();
374 
375  i = nLevels;
376 
377  while (i > 0)
378  {
379  i--;
380  if (values[i] > value)
381  {
382  //outcount<<i<<" "<<VibrationEnergy(i)<<G4endl;
383  return i;
384  }
385  value -= values[i];
386  }
387 
388  //outcount<<0<<" "<<VibrationEnergy(0)<<G4endl;
389 
390  return 0;
391 }
392 
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
394 
395 G4double G4DNASancheExcitationModel::Sum(G4double k)
396 {
397  G4double totalCrossSection = 0.;
398 
399  for (G4int i = 0; i < nLevels; i++)
400  {
401  totalCrossSection += PartialCrossSection(k, i);
402  }
403 
404  return totalCrossSection;
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408 
409 G4double G4DNASancheExcitationModel::LinInterpolate(G4double e1,
410  G4double e2,
411  G4double e,
412  G4double xs1,
413  G4double xs2)
414 {
415  G4double a = (xs2 - xs1) / (e2 - e1);
416  G4double b = xs2 - a * e2;
417  G4double value = a * e + b;
418  // G4cout<<"interP >> "<<e1<<" "<<e2<<" "<<e<<" "<<xs1<<" "<<xs2<<" "<<a<<" "<<b<<" "<<value<<G4endl;
419 
420  return value;
421 }
422 
G4DNASancheExcitationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNASancheExcitationModel")
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:602
size_t GetIndex() const
Definition: G4Material.hh:262
G4double PartialCrossSection(G4double energy, G4int level)
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
string material
Definition: eplot.py:19
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4ThreeVector & GetMomentumDirection() const
static constexpr double cm
Definition: G4SIunits.hh:119
G4ParticleChangeForGamma * fParticleChangeForGamma
static constexpr double eV
Definition: G4SIunits.hh:215
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static G4DNAMolecularMaterial * Instance()
tuple t1
Definition: plottest35.py:33
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:739
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132