ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNABornExcitationModel2.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNABornExcitationModel2.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 "G4PhysicsTable.hh"
33 #include "G4PhysicsVector.hh"
34 #include "G4UnitsTable.hh"
35 #include <map>
36 
37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38 
39 using namespace std;
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
44  const G4String& nam) :
45 G4VEmModel(nam), isInitialised(false), fTableData(0)
46 {
48  fHighEnergy = 0;
49  fLowEnergy = 0;
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  if (verboseLevel > 0)
61  {
62  G4cout << "Born excitation model is constructed " << G4endl;
63  }
66  fTotalXS = 0;
67  fTableData = 0;
68 
69  // Selection of stationary mode
70 
71  statCode = false;
72 }
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
77 {
78  // Cross section
79  if (fTableData)
80  delete fTableData;
81 }
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
86  const G4DataVector& /*cuts*/)
87 {
88 
89  if (verboseLevel > 3)
90  {
91  G4cout << "Calling G4DNABornExcitationModel2::Initialise()" << G4endl;
92  }
93 
94  if(fParticleDefinition != 0 && fParticleDefinition != particle)
95  {
96  G4Exception("G4DNABornExcitationModel2::Initialise","em0001",
97  FatalException,"Model already initialized for another particle type.");
98  }
99 
101 
102  std::ostringstream fullFileName;
103  char *path = getenv("G4LEDATA");
104 
105  if(G4String(path) == "")
106  {
107  G4Exception("G4DNABornExcitationModel2::Initialise","G4LEDATA-CHECK",
108  FatalException, "G4LEDATA not defined in environment variables");
109  }
110 
111  fullFileName << path;
112 
113  if(particle->GetParticleName() == "e-")
114  {
115  fullFileName << "/dna/bornExcitation-e.dat";
116  fLowEnergy = 9*eV;
117  fHighEnergy = 1*MeV;
118  }
119  else if(particle->GetParticleName() == "proton")
120  {
121  fullFileName << "/dna/bornExcitation-p.dat";
122  fLowEnergy = 500. * keV;
123  fHighEnergy = 100. * MeV;
124  }
125 
128 
129  //G4double scaleFactor = (1.e-22 / 3.343) * m*m;
130 
131  fTableData = new G4PhysicsTable();
132  fTableData->RetrievePhysicsTable(fullFileName.str().c_str(), true);
133  for(size_t level = 0; level<fTableData->size(); ++level)
134  {
135  //(*fTableData)(level)->ScaleVector(1,scaleFactor);
136  (*fTableData)(level)->SetSpline(true);
137  }
138 
139  size_t finalBin_i = 2000;
140  G4double E_min = fLowEnergy;
141  G4double E_max = fHighEnergy;
142  fTotalXS = new G4PhysicsLogVector(E_min, E_max, finalBin_i);
143  fTotalXS->SetSpline(true);
145  G4double finalXS;
146 
147  for(size_t energy_i = 0; energy_i < finalBin_i; ++energy_i)
148  {
149  energy = fTotalXS->Energy(energy_i);
150  finalXS = 0;
151 
152  for(size_t level = 0; level<fTableData->size(); ++level)
153  {
154  finalXS += (*fTableData)(level)->Value(energy);
155  }
156  fTotalXS->PutValue(energy_i, finalXS);
157  //G4cout << "energy = " << energy << " " << fTotalXS->Value(energy)
158  // << " " << energy_i << " " << finalXS << G4endl;
159  }
160 
161  // for(energy = LowEnergyLimit() ; energy < HighEnergyLimit() ; energy += 1*pow(10,log10(energy)))
162  // {
163  // G4cout << "energy = " << energy << " " << fTotalXS->Value(energy) << G4endl;
164  // }
165 
166  if( verboseLevel>0 )
167  {
168  G4cout << "Born excitation model is initialized " << G4endl
169  << "Energy range: "
170  << LowEnergyLimit() / eV << " eV - "
171  << HighEnergyLimit() / keV << " keV for "
172  << particle->GetParticleName()
173  << G4endl;
174  }
175 
176  // Initialize water density pointer
178 
179  if (isInitialised)
180  { return;}
182  isInitialised = true;
183 }
184 
185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
186 
188  const G4ParticleDefinition* particleDefinition,
189  G4double ekin,
190  G4double,
191  G4double)
192 {
193  if (verboseLevel > 3)
194  {
195  G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel2"
196  << G4endl;
197  }
198 
199  if(particleDefinition != fParticleDefinition) return 0;
200 
201  // Calculate total cross section for model
202 
203  G4double sigma=0;
204 
205  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
206 
207  if (ekin >= fLowEnergy && ekin <= fHighEnergy)
208  {
209  sigma = fTotalXS->Value(ekin, fLastBinCallForFinalXS);
210 
211  // for(size_t i = 0; i < 5; ++i)
212  // sigma += (*fTableData)[i]->Value(ekin);
213 
214  if(sigma == 0)
215  {
216  G4cerr << "PROBLEM SIGMA = 0 at " << G4BestUnit(ekin, "Energy")<< G4endl;
217  }
218  }
219 
220  if (verboseLevel > 2)
221  {
222  G4cout << "__________________________________" << G4endl;
223  G4cout << "G4DNABornExcitationModel2 - XS INFO START" << G4endl;
224  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
225  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
226  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
227  G4cout << "G4DNABornExcitationModel2 - XS INFO END" << G4endl;
228  }
229 
230  return sigma*waterDensity;
231 }
232 
233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
234 
235 void G4DNABornExcitationModel2::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
236  const G4MaterialCutsCouple* /*couple*/,
237  const G4DynamicParticle* aDynamicParticle,
238  G4double,
239  G4double)
240 {
241 
242  if (verboseLevel > 3)
243  {
244  G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel2"
245  << G4endl;
246  }
247 
248  G4double k = aDynamicParticle->GetKineticEnergy();
249 
250  G4int level = RandomSelect(k);
251  G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
252  G4double newEnergy = k - excitationEnergy;
253 
254  if (newEnergy > 0)
255  {
257 
260 
262  }
263 
264  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
266  level,
267  theIncomingTrack);
268 }
269 
270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
271 
273  G4int level,
275  G4double kineticEnergy)
276 {
277  if (fParticleDefinition != particle)
278  {
279  G4Exception("G4DNABornExcitationModel2::GetPartialCrossSection",
280  "bornParticleType",
282  "Model initialized for another particle type.");
283  }
284 
285  return (*fTableData)(level)->Value(kineticEnergy);
286 }
287 
288 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
289 
291 {
292  const size_t n(fTableData->size());
293  size_t i(n);
294 
296 
297  value *= G4UniformRand();
298  i = n;
299 
300  G4double partialXS;
301 
302  while (i > 0)
303  {
304  i--;
305 
306  partialXS = (*fTableData)(i)->Value(k);
307  if (partialXS > value)
308  {
309  return i;
310  }
311  value -= partialXS;
312  }
313 
314  return 0;
315 }
316