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