ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNABornExcitationModel1.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNABornExcitationModel1.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 
29 #include "G4SystemOfUnits.hh"
30 #include "G4DNAChemistryManager.hh"
32 #include <map>
33 
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 
36 using namespace std;
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 
41  const G4String& nam) :
42 G4VEmModel(nam), isInitialised(false), fTableData(0)
43 {
45  fHighEnergy = 0;
46  fLowEnergy = 0;
48 
49  verboseLevel = 0;
50  // Verbosity scale:
51  // 0 = nothing
52  // 1 = warning for energy non-conservation
53  // 2 = details of energy budget
54  // 3 = calculation of cross sections, file openings, sampling of atoms
55  // 4 = entering in methods
56 
57  if (verboseLevel > 0)
58  {
59  G4cout << "Born excitation model is constructed " << G4endl;
60  }
62 
63  // Selection of stationary mode
64 
65  statCode = false;
66 }
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 
71 {
72  // Cross section
73  if (fTableData)
74  delete fTableData;
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 
80  const G4DataVector& /*cuts*/)
81 {
82 
83  if (verboseLevel > 3)
84  {
85  G4cout << "Calling G4DNABornExcitationModel1::Initialise()" << G4endl;
86  }
87 
88  if(fParticleDefinition != 0 && fParticleDefinition != particle)
89  {
90  G4Exception("G4DNABornExcitationModel1::Initialise","em0001",
91  FatalException,"Model already initialized for another particle type.");
92  }
93 
95 
96  if(particle->GetParticleName() == "e-")
97  {
98  fTableFile = "dna/sigma_excitation_e_born";
99  fLowEnergy = 9*eV;
100  fHighEnergy = 1*MeV;
101  }
102  else if(particle->GetParticleName() == "proton")
103  {
104  fTableFile = "dna/sigma_excitation_p_born";
105  fLowEnergy = 500. * keV;
106  fHighEnergy = 100. * MeV;
107  }
108 
111 
112  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
115 
116  if( verboseLevel>0 )
117  {
118  G4cout << "Born excitation model is initialized " << G4endl
119  << "Energy range: "
120  << LowEnergyLimit() / eV << " eV - "
121  << HighEnergyLimit() / keV << " keV for "
122  << particle->GetParticleName()
123  << G4endl;
124  }
125 
126  // Initialize water density pointer
128 
129  if (isInitialised)
130  { return;}
132  isInitialised = true;
133 }
134 
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
136 
138  const G4ParticleDefinition* particleDefinition,
139  G4double ekin,
140  G4double,
141  G4double)
142 {
143  if (verboseLevel > 3)
144  {
145  G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel1"
146  << G4endl;
147  }
148 
149  if(particleDefinition != fParticleDefinition) return 0;
150 
151  // Calculate total cross section for model
152 
153  G4double sigma=0;
154 
155  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
156 
157  if (ekin >= fLowEnergy && ekin <= fHighEnergy)
158  {
159  sigma = fTableData->FindValue(ekin);
160  }
161 
162  if (verboseLevel > 2)
163  {
164  G4cout << "__________________________________" << G4endl;
165  G4cout << "G4DNABornExcitationModel1 - XS INFO START" << G4endl;
166  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
167  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
168  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
169  G4cout << "G4DNABornExcitationModel1 - XS INFO END" << G4endl;
170  }
171 
172  return sigma*waterDensity;
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
176 
177 void G4DNABornExcitationModel1::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
178  const G4MaterialCutsCouple* /*couple*/,
179  const G4DynamicParticle* aDynamicParticle,
180  G4double,
181  G4double)
182 {
183 
184  if (verboseLevel > 3)
185  {
186  G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel1"
187  << G4endl;
188  }
189 
190  G4double k = aDynamicParticle->GetKineticEnergy();
191 
192  G4int level = RandomSelect(k);
193  G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
194  G4double newEnergy = k - excitationEnergy;
195 
196  if (newEnergy > 0)
197  {
199 
202 
204  }
205 
206  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
208  level,
209  theIncomingTrack);
210 }
211 
212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
213 
215  G4int level,
217  G4double kineticEnergy)
218 {
219  if (fParticleDefinition != particle)
220  {
221  G4Exception("G4DNABornExcitationModel1::GetPartialCrossSection",
222  "bornParticleType",
224  "Model initialized for another particle type.");
225  }
226 
227  return fTableData->GetComponent(level)->FindValue(kineticEnergy);
228 }
229 
230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
231 
233 {
234  G4int level = 0;
235 
236  G4double* valuesBuffer = new G4double[fTableData->NumberOfComponents()];
237  const size_t n(fTableData->NumberOfComponents());
238  size_t i(n);
239  G4double value = 0.;
240 
241  while (i > 0)
242  {
243  i--;
244  valuesBuffer[i] = fTableData->GetComponent(i)->FindValue(k);
245  value += valuesBuffer[i];
246  }
247 
248  value *= G4UniformRand();
249  i = n;
250 
251  while (i > 0)
252  {
253  i--;
254 
255  if (valuesBuffer[i] > value)
256  {
257  delete[] valuesBuffer;
258  return i;
259  }
260  value -= valuesBuffer[i];
261  }
262 
263  if (valuesBuffer)
264  delete[] valuesBuffer;
265 
266  return level;
267 }
268