Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAChampionElasticModel.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 #include "G4DNAChampionElasticModel.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "G4DNAMolecularMaterial.hh"
33 
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35 
36 using namespace std;
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 
40 G4DNAChampionElasticModel::G4DNAChampionElasticModel(const G4ParticleDefinition*,
41  const G4String& nam)
42 :G4VEmModel(nam),isInitialised(false)
43 {
44 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
45 
46  killBelowEnergy = 7.4*eV;
47  lowEnergyLimit = 0 * eV;
48  highEnergyLimit = 1. * MeV;
49  SetLowEnergyLimit(lowEnergyLimit);
50  SetHighEnergyLimit(highEnergyLimit);
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 << "Champion Elastic model is constructed " << G4endl
63  << "Energy range: "
64  << lowEnergyLimit / eV << " eV - "
65  << highEnergyLimit / MeV << " MeV"
66  << G4endl;
67  }
68  fParticleChangeForGamma = 0;
69  fpMolWaterDensity = 0;
70 }
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
74 G4DNAChampionElasticModel::~G4DNAChampionElasticModel()
75 {
76  // For total cross section
77 
78  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
79  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
80  {
81  G4DNACrossSectionDataSet* table = pos->second;
82  delete table;
83  }
84 
85  // For final state
86 
87  eVecm.clear();
88 
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
93 void G4DNAChampionElasticModel::Initialise(const G4ParticleDefinition* /*particle*/,
94  const G4DataVector& /*cuts*/)
95 {
96 
97  if (verboseLevel > 3)
98  G4cout << "Calling G4DNAChampionElasticModel::Initialise()" << G4endl;
99 
100  // Energy limits
101 
102  if (LowEnergyLimit() < lowEnergyLimit)
103  {
104  G4cout << "G4DNAChampionElasticModel: low energy limit increased from " <<
105  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
106  SetLowEnergyLimit(lowEnergyLimit);
107  }
108 
109  if (HighEnergyLimit() > highEnergyLimit)
110  {
111  G4cout << "G4DNAChampionElasticModel: high energy limit decreased from " <<
112  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
113  SetHighEnergyLimit(highEnergyLimit);
114  }
115 
116  // Reading of data files
117 
118  G4double scaleFactor = 1e-16*cm*cm;
119 
120  G4String fileElectron("dna/sigma_elastic_e_champion");
121 
122  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
123  G4String electron;
124 
125  // *** ELECTRON
126 
127  // For total cross section
128 
129  electron = electronDef->GetParticleName();
130 
131  tableFile[electron] = fileElectron;
132 
133  G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
134  tableE->LoadData(fileElectron);
135  tableData[electron] = tableE;
136 
137  // For final state
138 
139  char *path = getenv("G4LEDATA");
140 
141  if (!path)
142  {
143  G4Exception("G4ChampionElasticModel::Initialise","em0006",
144  FatalException,"G4LEDATA environment variable not set.");
145  return;
146  }
147 
148  std::ostringstream eFullFileName;
149  eFullFileName << path << "/dna/sigmadiff_cumulated_elastic_e_champion.dat";
150  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
151 
152  if (!eDiffCrossSection)
153  G4Exception("G4DNAChampionElasticModel::Initialise","em0003",
154  FatalException,"Missing data file:/dna/sigmadiff_cumulated_elastic_e_champion.dat");
155 
156  eTdummyVec.push_back(0.);
157 
158  while(!eDiffCrossSection.eof())
159  {
160  double tDummy;
161  double eDummy;
162  eDiffCrossSection>>tDummy>>eDummy;
163 
164  // SI : mandatory eVecm initialization
165 
166  if (tDummy != eTdummyVec.back())
167  {
168  eTdummyVec.push_back(tDummy);
169  eVecm[tDummy].push_back(0.);
170  }
171 
172  eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
173 
174  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
175 
176  }
177 
178  // End final state
179 
180  if (verboseLevel > 2)
181  G4cout << "Loaded cross section files for Champion Elastic model" << G4endl;
182 
183  if( verboseLevel>0 )
184  {
185  G4cout << "Champion Elastic model is initialized " << G4endl
186  << "Energy range: "
187  << LowEnergyLimit() / eV << " eV - "
188  << HighEnergyLimit() / MeV << " MeV"
189  << G4endl;
190  }
191 
192  // Initialize water density pointer
193  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
194 
195  if (isInitialised) { return; }
196  fParticleChangeForGamma = GetParticleChangeForGamma();
197  isInitialised = true;
198 
199 }
200 
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
202 
203 G4double G4DNAChampionElasticModel::CrossSectionPerVolume(const G4Material* material,
204  const G4ParticleDefinition* p,
205  G4double ekin,
206  G4double,
207  G4double)
208 {
209  if (verboseLevel > 3)
210  G4cout << "Calling CrossSectionPerVolume() of G4DNAChampionElasticModel" << G4endl;
211 
212  // Calculate total cross section for model
213 
214  G4double sigma=0;
215 
216  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
217 
218  if(waterDensity!= 0.0)
219 // if (material == nistwater || material->GetBaseMaterial() == nistwater)
220  {
221  const G4String& particleName = p->GetParticleName();
222 
223  if (ekin < highEnergyLimit)
224  {
225  //SI : XS must not be zero otherwise sampling of secondaries method ignored
226  if (ekin < killBelowEnergy) return DBL_MAX;
227  //
228 
229  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
230  pos = tableData.find(particleName);
231 
232  if (pos != tableData.end())
233  {
234  G4DNACrossSectionDataSet* table = pos->second;
235  if (table != 0)
236  {
237  sigma = table->FindValue(ekin);
238  }
239  }
240  else
241  {
242  G4Exception("G4DNAChampionElasticModel::ComputeCrossSectionPerVolume","em0002",
243  FatalException,"Model not applicable to particle type.");
244  }
245  }
246 
247  if (verboseLevel > 2)
248  {
249  G4cout << "__________________________________" << G4endl;
250  G4cout << "°°° G4DNAChampionElasticModel - XS INFO START" << G4endl;
251  G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
252  G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
253  G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
254  // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
255  G4cout << "°°° G4DNAChampionElasticModel - XS INFO END" << G4endl;
256  }
257 
258  }
259 
260  return sigma*waterDensity;
261 // return sigma*material->GetAtomicNumDensityVector()[1];
262 }
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
265 
266 void G4DNAChampionElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
267  const G4MaterialCutsCouple* /*couple*/,
268  const G4DynamicParticle* aDynamicElectron,
269  G4double,
270  G4double)
271 {
272 
273  if (verboseLevel > 3)
274  G4cout << "Calling SampleSecondaries() of G4DNAChampionElasticModel" << G4endl;
275 
276  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
277 
278  if (electronEnergy0 < killBelowEnergy)
279  {
280  fParticleChangeForGamma->SetProposedKineticEnergy(0.);
281  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
282  fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
283  return ;
284  }
285 
286  if (electronEnergy0>= killBelowEnergy && electronEnergy0 < highEnergyLimit)
287  {
288 
289  G4double cosTheta = RandomizeCosTheta(electronEnergy0);
290 
291  G4double phi = 2. * pi * G4UniformRand();
292 
293  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
294  G4ThreeVector xVers = zVers.orthogonal();
295  G4ThreeVector yVers = zVers.cross(xVers);
296 
297  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
298  G4double yDir = xDir;
299  xDir *= std::cos(phi);
300  yDir *= std::sin(phi);
301 
302  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
303 
304  fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
305 
306  fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
307  }
308 
309 }
310 
311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
312 
313 G4double G4DNAChampionElasticModel::Theta
314  (G4ParticleDefinition * particleDefinition, G4double k, G4double integrDiff)
315 {
316  G4double theta = 0.;
317  G4double valueT1 = 0;
318  G4double valueT2 = 0;
319  G4double valueE21 = 0;
320  G4double valueE22 = 0;
321  G4double valueE12 = 0;
322  G4double valueE11 = 0;
323  G4double xs11 = 0;
324  G4double xs12 = 0;
325  G4double xs21 = 0;
326  G4double xs22 = 0;
327 
328  if (particleDefinition == G4Electron::ElectronDefinition())
329  {
330  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
331  std::vector<double>::iterator t1 = t2-1;
332 
333  std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(), integrDiff);
334  std::vector<double>::iterator e11 = e12-1;
335 
336  std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(), integrDiff);
337  std::vector<double>::iterator e21 = e22-1;
338 
339  valueT1 =*t1;
340  valueT2 =*t2;
341  valueE21 =*e21;
342  valueE22 =*e22;
343  valueE12 =*e12;
344  valueE11 =*e11;
345 
346  xs11 = eDiffCrossSectionData[valueT1][valueE11];
347  xs12 = eDiffCrossSectionData[valueT1][valueE12];
348  xs21 = eDiffCrossSectionData[valueT2][valueE21];
349  xs22 = eDiffCrossSectionData[valueT2][valueE22];
350  }
351 
352  if (xs11==0 && xs12==0 && xs21==0 && xs22==0) return (0.);
353 
354  theta = QuadInterpolator ( valueE11, valueE12,
355  valueE21, valueE22,
356  xs11, xs12,
357  xs21, xs22,
358  valueT1, valueT2,
359  k, integrDiff );
360 
361  return theta;
362 }
363 
364 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
365 
366 G4double G4DNAChampionElasticModel::LinLogInterpolate(G4double e1,
367  G4double e2,
368  G4double e,
369  G4double xs1,
370  G4double xs2)
371 {
372  G4double d1 = std::log(xs1);
373  G4double d2 = std::log(xs2);
374  G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
375  return value;
376 }
377 
378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
379 
380 G4double G4DNAChampionElasticModel::LinLinInterpolate(G4double e1,
381  G4double e2,
382  G4double e,
383  G4double xs1,
384  G4double xs2)
385 {
386  G4double d1 = xs1;
387  G4double d2 = xs2;
388  G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
389  return value;
390 }
391 
392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
393 
394 G4double G4DNAChampionElasticModel::LogLogInterpolate(G4double e1,
395  G4double e2,
396  G4double e,
397  G4double xs1,
398  G4double xs2)
399 {
400  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
401  G4double b = std::log10(xs2) - a*std::log10(e2);
402  G4double sigma = a*std::log10(e) + b;
403  G4double value = (std::pow(10.,sigma));
404  return value;
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408 
409 
410 G4double G4DNAChampionElasticModel::QuadInterpolator(G4double e11, G4double e12,
411  G4double e21, G4double e22,
412  G4double xs11, G4double xs12,
413  G4double xs21, G4double xs22,
414  G4double t1, G4double t2,
415  G4double t, G4double e)
416 {
417  // Log-Log
418 /*
419  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
420  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
421  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
422 
423 
424  // Lin-Log
425  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
426  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
427  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
428 */
429 
430  // Lin-Lin
431  G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
432  G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
433  G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
434 
435  return value;
436 }
437 
438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
439 
440 G4double G4DNAChampionElasticModel::RandomizeCosTheta(G4double k)
441 {
442 
443  G4double integrdiff=0;
444  G4double uniformRand=G4UniformRand();
445  integrdiff = uniformRand;
446 
447  G4double theta=0.;
448  G4double cosTheta=0.;
449  theta = Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
450 
451  cosTheta= std::cos(theta*pi/180);
452 
453  return cosTheta;
454 }
Generated on Sat May 25 2013 14:33:30 for Geant4 by   doxygen 1.8.4