Geant4  9.6.p02
 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$
27 //
28 
29 // Created by Z. Francis
30 
31 #include "G4DNASancheExcitationModel.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4DNAMolecularMaterial.hh"
34 
35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
36 
37 using namespace std;
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40 
41 G4DNASancheExcitationModel::G4DNASancheExcitationModel(const G4ParticleDefinition*,
42  const G4String& nam)
43  :G4VEmModel(nam),isInitialised(false)
44 {
45  // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
46  fpWaterDensity = 0;
47 
48  lowEnergyLimit = 2 * eV;
49  highEnergyLimit = 100 * eV;
50  SetLowEnergyLimit(lowEnergyLimit);
51  SetHighEnergyLimit(highEnergyLimit);
52  nLevels = 9;
53 
54  verboseLevel= 0;
55  // Verbosity scale:
56  // 0 = nothing
57  // 1 = warning for energy non-conservation
58  // 2 = details of energy budget
59  // 3 = calculation of cross sections, file openings, sampling of atoms
60  // 4 = entering in methods
61 
62  if (verboseLevel > 0)
63  {
64  G4cout << "Sanche Excitation model is constructed " << G4endl
65  << "Energy range: "
66  << lowEnergyLimit / eV << " eV - "
67  << highEnergyLimit / eV << " eV"
68  << G4endl;
69  }
70  fParticleChangeForGamma = 0;
71  fpWaterDensity = 0;
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
76 G4DNASancheExcitationModel::~G4DNASancheExcitationModel()
77 {}
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
81 void G4DNASancheExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
82  const G4DataVector& /*cuts*/)
83 {
84 
85 
86  if (verboseLevel > 3)
87  G4cout << "Calling G4DNASancheExcitationModel::Initialise()" << G4endl;
88 
89  // Energy limits
90 
91  if (LowEnergyLimit() < lowEnergyLimit)
92  {
93  G4cout << "G4DNASancheExcitationModel: low energy limit increased from " <<
94  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
95  SetLowEnergyLimit(lowEnergyLimit);
96  }
97 
98  if (HighEnergyLimit() > highEnergyLimit)
99  {
100  G4cout << "G4DNASancheExcitationModel: high energy limit decreased from " <<
101  HighEnergyLimit()/eV << " eV to " << highEnergyLimit/eV << " eV" << G4endl;
102  SetHighEnergyLimit(highEnergyLimit);
103  }
104 
105  //
106 
107  if (verboseLevel > 0)
108  {
109  G4cout << "Sanche Excitation model is initialized " << G4endl
110  << "Energy range: "
111  << LowEnergyLimit() / eV << " eV - "
112  << HighEnergyLimit() / eV << " eV"
113  << G4endl;
114  }
115 
116  // Initialize water density pointer
117  fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
118 
119  if (isInitialised) { return; }
120  fParticleChangeForGamma = GetParticleChangeForGamma();
121  isInitialised = true;
122 
123  char *path = getenv("G4LEDATA");
124  std::ostringstream eFullFileName;
125  eFullFileName << path << "/dna/sigma_excitationvib_e_sanche.dat";
126  std::ifstream input(eFullFileName.str().c_str());
127 
128  if (!input)
129  {
130  G4Exception("G4DNASancheExcitationModel::Initialise","em0003",
131  FatalException,"Missing data file:/dna/sigma_excitationvib_e_sanche.dat");
132  }
133 
134  while(!input.eof())
135  {
136  double t;
137  input>>t;
138  tdummyVec.push_back(t);
139  input>>map1[t][0]>>map1[t][1]>>map1[t][2]>>map1[t][3]>>map1[t][4]>>map1[t][5]>>map1[t][6]>>map1[t][7]>>map1[t][8];
140  //G4cout<<t<<" "<<map1[t][0]<<map1[t][1]<<map1[t][2]<<map1[t][3]<<map1[t][4]<<map1[t][5]<<map1[t][6]<<map1[t][7]<<map1[t][8]<<G4endl;
141  }
142 }
143 
144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 
146 G4double G4DNASancheExcitationModel::CrossSectionPerVolume(const G4Material* material,
147  const G4ParticleDefinition* particleDefinition,
148  G4double ekin,
149  G4double,
150  G4double)
151 {
152  if (verboseLevel > 3)
153  G4cout << "Calling CrossSectionPerVolume() of G4DNASancheExcitationModel" << G4endl;
154 
155  // Calculate total cross section for model
156 
157  G4double sigma=0;
158 
159  G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
160 
161  if(waterDensity!= 0.0)
162  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
163  {
164 
165  if (particleDefinition == G4Electron::ElectronDefinition())
166  {
167  if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
168  {
169  sigma = Sum(ekin);
170  }
171  }
172 
173  if (verboseLevel > 2)
174  {
175  G4cout << "__________________________________" << G4endl;
176  G4cout << "°°° G4DNASancheExcitationModel - XS INFO START" << G4endl;
177  G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
178  G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
179  G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
180  // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
181  G4cout << "°°° G4DNASancheExcitationModel - XS INFO END" << G4endl;
182  }
183 
184  } // if water
185 
186 
187  // return sigma*2*material->GetAtomicNumDensityVector()[1];
188  return sigma*2*waterDensity;
189  // see papers for factor 2 description
190 
191 }
192 
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
194 
195 void G4DNASancheExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
196  const G4MaterialCutsCouple*,
197  const G4DynamicParticle* aDynamicElectron,
198  G4double,
199  G4double)
200 {
201 
202  if (verboseLevel > 3)
203  G4cout << "Calling SampleSecondaries() of G4DNASancheExcitationModel" << G4endl;
204 
205  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
206  G4int level = RandomSelect(electronEnergy0);
207  G4double excitationEnergy = VibrationEnergy(level); // levels go from 0 to 8
208  G4double newEnergy = electronEnergy0 - excitationEnergy;
209 
210  /*
211  if (electronEnergy0 < highEnergyLimit)
212  {
213  if (newEnergy >= lowEnergyLimit)
214  {
215  fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
216  fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
217  fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
218  }
219 
220  else
221  {
222  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
223  fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
224  }
225  }
226 */
227 
228  if (electronEnergy0 < highEnergyLimit && newEnergy>0.)
229  {
230  fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
231  fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
232  fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
233  }
234 
235  //
236 }
237 
238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
239 
240 G4double G4DNASancheExcitationModel::PartialCrossSection(G4double t, G4int level)
241 {
242  std::vector<double>::iterator t2 = std::upper_bound(tdummyVec.begin(),tdummyVec.end(), t/eV);
243  std::vector<double>::iterator t1 = t2-1;
244 
245  double sigma = LinInterpolate((*t1), (*t2), t/eV, map1[*t1][level], map1[*t2][level]);
246  sigma*=1e-16*cm*cm;
247  if(sigma==0.)sigma=1e-30;
248  return (sigma);
249 }
250 
251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
252 
253 G4double G4DNASancheExcitationModel::VibrationEnergy(G4int level)
254 {
255  G4double energies[9] = {0.01, 0.024, 0.061, 0.092, 0.204, 0.417, 0.460, 0.500, 0.835};
256  return(energies[level]*eV);
257 }
258 
259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
260 
261 G4int G4DNASancheExcitationModel::RandomSelect(G4double k)
262 {
263 
264  // Level Selection Counting can be done here !
265 
266  G4int i = nLevels;
267  G4double value = 0.;
268  std::deque<double> values;
269 
270  while (i > 0)
271  {
272  i--;
273  G4double partial = PartialCrossSection(k,i);
274  values.push_front(partial);
275  value += partial;
276  }
277 
278  value *= G4UniformRand();
279 
280  i = nLevels;
281 
282  while (i > 0)
283  {
284  i--;
285  if (values[i] > value)
286  {
287  //outcount<<i<<" "<<VibrationEnergy(i)<<G4endl;
288  return i;
289  }
290  value -= values[i];
291  }
292 
293  //outcount<<0<<" "<<VibrationEnergy(0)<<G4endl;
294 
295  return 0;
296 }
297 
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
299 
300 G4double G4DNASancheExcitationModel::Sum(G4double k)
301 {
302  G4double totalCrossSection = 0.;
303 
304  for (G4int i=0; i<nLevels; i++)
305  {
306  totalCrossSection += PartialCrossSection(k,i);
307  }
308  return totalCrossSection;
309 }
310 
311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312 
313 G4double G4DNASancheExcitationModel::LinInterpolate(G4double e1,
314  G4double e2,
315  G4double e,
316  G4double xs1,
317  G4double xs2)
318 {
319  G4double a = (xs2 - xs1) / (e2 - e1);
320  G4double b = xs2 - a*e2;
321  G4double value = a*e + b;
322  // G4cout<<"interP >> "<<e1<<" "<<e2<<" "<<e<<" "<<xs1<<" "<<xs2<<" "<<a<<" "<<b<<" "<<value<<G4endl;
323 
324  return value;
325 }
326 
Generated on Sat May 25 2013 14:33:31 for Geant4 by   doxygen 1.8.4