ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAEmfietzoglouExcitationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAEmfietzoglouExcitationModel.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 // 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 {
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  if( verboseLevel>0 )
62  {
63  G4cout << "Emfietzoglou excitation model is constructed " << G4endl;
64  }
66 
69 
70  // Selection of stationary mode
71  statCode = false;
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
77 {
78  // Cross section
79 
80  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
81  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
82  {
83  G4DNACrossSectionDataSet* table = pos->second;
84  delete table;
85  }
86 
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
92  const G4DataVector& /*cuts*/)
93 {
94 
95  if (verboseLevel > 3)
96  G4cout << "Calling G4DNAEmfietzoglouExcitationModel::Initialise()" << G4endl;
97 
98  G4String fileElectron("dna/sigma_excitation_e_emfietzoglou");
99 
101 
103 
104  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
105 
106  // *** ELECTRON
107 
108  electron = electronDef->GetParticleName();
109 
110  tableFile[electron] = fileElectron;
111 
112  // Cross section
113 
115  tableE->LoadData(fileElectron);
116 
117  tableData[electron] = tableE;
118 
119  //
120 
121  if( verboseLevel>0 )
122  {
123  G4cout << "Emfietzoglou excitation model is initialized " << G4endl
124  << "Energy range: "
125  << LowEnergyLimit() / eV << " eV - "
126  << HighEnergyLimit() / keV << " keV for "
127  << particle->GetParticleName()
128  << G4endl;
129  }
130 
131  // Initialize water density pointer
133 
134  if (isInitialised) return;
136  isInitialised = true;
137 }
138 
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140 
142  const G4ParticleDefinition* particleDefinition,
143  G4double ekin,
144  G4double,
145  G4double)
146 {
147  if (verboseLevel > 3)
148  G4cout << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouExcitationModel" << G4endl;
149 
150  if (particleDefinition != G4Electron::ElectronDefinition()) return 0;
151 
152  // Calculate total cross section for model
153 
154  G4double sigma=0;
155 
156  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
157 
158  const G4String& particleName = particleDefinition->GetParticleName();
159 
160  if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit())
161  {
162  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
163  pos = tableData.find(particleName);
164 
165  if (pos != tableData.end())
166  {
167  G4DNACrossSectionDataSet* table = pos->second;
168  if (table != 0) sigma = table->FindValue(ekin);
169  }
170  else
171  {
172  G4Exception("G4DNAEmfietzoglouExcitationModel::CrossSectionPerVolume","em0002",
173  FatalException,"Model not applicable to particle type.");
174  }
175  }
176 
177  if (verboseLevel > 2)
178  {
179  G4cout << "__________________________________" << G4endl;
180  G4cout << "G4DNAEmfietzoglouExcitationModel - XS INFO START" << G4endl;
181  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
182  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
183  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
184  //G4cout << " Cross section per water molecule (cm^-1)=" <<
186  G4cout << "G4DNAEmfietzoglouExcitationModel - XS INFO END" << G4endl;
187  }
188 
189  return sigma*waterDensity;
190 }
191 
192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
193 
194 void G4DNAEmfietzoglouExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
195  const G4MaterialCutsCouple* /*couple*/,
196  const G4DynamicParticle* aDynamicParticle,
197  G4double,
198  G4double)
199 {
200 
201  if (verboseLevel > 3)
202  G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouExcitationModel" << G4endl;
203 
204  G4double k = aDynamicParticle->GetKineticEnergy();
205 
206  const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
207 
208  G4int level = RandomSelect(k,particleName);
209  G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
210  G4double newEnergy = k - excitationEnergy;
211 
212  if (newEnergy > 0)
213  {
215 
218 
220  }
221 
222  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
224  level,
225  theIncomingTrack);
226 }
227 
228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
229 
231 {
232 
233  G4int level = 0;
234 
235  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
236  pos = tableData.find(particle);
237 
238  if (pos != tableData.end())
239  {
240  G4DNACrossSectionDataSet* table = pos->second;
241 
242  if (table != 0)
243  {
244  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
245  const size_t n(table->NumberOfComponents());
246  size_t i(n);
247  G4double value = 0.;
248 
249  //Check reading of initial xs file
250  //G4cout << table->GetComponent(0)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
251  //G4cout << table->GetComponent(1)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
252  //G4cout << table->GetComponent(2)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
253  //G4cout << table->GetComponent(3)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
254  //G4cout << table->GetComponent(4)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
255  //G4cout << table->GetComponent(5)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
256  //G4cout << table->GetComponent(6)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
257  //abort();
258 
259  while (i>0)
260  {
261  i--;
262  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
263  value += valuesBuffer[i];
264  }
265 
266  value *= G4UniformRand();
267 
268  i = n;
269 
270  while (i > 0)
271  {
272  i--;
273 
274  if (valuesBuffer[i] > value)
275  {
276  delete[] valuesBuffer;
277  return i;
278  }
279  value -= valuesBuffer[i];
280  }
281 
282  if (valuesBuffer) delete[] valuesBuffer;
283 
284  }
285  }
286  else
287  {
288  G4Exception("G4DNAEmfietzoglouExcitationModel::RandomSelect","em0002",
289  FatalException,"Model not applicable to particle type.");
290  }
291  return level;
292 }