ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNACPA100ElasticModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNACPA100ElasticModel.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 // CPA100 elastic model class for electrons
27 //
28 // Based on the work of M. Terrissol and M. C. Bordage
29 //
30 // Users are requested to cite the following papers:
31 // - M. Terrissol, A. Baudre, Radiat. Prot. Dosim. 31 (1990) 175-177
32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terrissol, X. Franceries,
33 // M. Bardies, N. Lampe, S. Incerti, Phys. Med. 32 (2016) 1833-1840
34 //
35 // Authors of this class:
36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bordes, S. Incerti
37 //
38 // 15.01.2014: creation
39 //
40 
41 #include "G4DNACPA100ElasticModel.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4DNAMolecularMaterial.hh"
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 
48 using namespace std;
49 
50 // #define CPA100_VERBOSE
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 
54 G4DNACPA100ElasticModel::G4DNACPA100ElasticModel(const G4ParticleDefinition*,
55  const G4String& nam)
56 :G4VEmModel(nam),isInitialised(false)
57 {
58 
59  SetLowEnergyLimit(11*eV);
60  SetHighEnergyLimit(255955*eV);
61 
62  verboseLevel= 0;
63  // Verbosity scale:
64  // 0 = nothing
65  // 1 = warning for energy non-conservation
66  // 2 = details of energy budget
67  // 3 = calculation of cross sections, file openings, sampling of atoms
68  // 4 = entering in methods
69 
70 #ifdef UEHARA_VERBOSE
71  if( verboseLevel>0 )
72  {
73  G4cout << "CPA100 Elastic model is constructed " << G4endl
74  << "Energy range: "
75  << LowEnergyLimit()/eV << " eV - "
76  << HighEnergyLimit()/ keV << " keV"
77  << G4endl;
78  }
79 #endif
80 
81  fParticleChangeForGamma = 0;
82  fpMolWaterDensity = 0;
83 
84  // Selection of stationary mode
85 
86  statCode = false;
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
91 G4DNACPA100ElasticModel::~G4DNACPA100ElasticModel()
92 {
93  // For total cross section
94 
95  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
96  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
97  {
98  G4DNACrossSectionDataSet* table = pos->second;
99  delete table;
100  }
101 
102  // For final state
103  eVecm.clear();
104 
105 }
106 
107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108 
109 void G4DNACPA100ElasticModel::Initialise(const G4ParticleDefinition*
110  particle,
111  const G4DataVector& /*cuts*/)
112 {
113 
114 #ifdef UEHARA_VERBOSE
115  if (verboseLevel > 3)
116  G4cout << "Calling G4DNACPA100ElasticModel::Initialise()" << G4endl;
117 #endif
118 
119  if(particle->GetParticleName() != "e-")
120  {
121  G4Exception("*** WARNING: the G4DNACPA100ElasticModel is "
122  "not intented to be used with another particle than the electron",
123  "",FatalException,"") ;
124  }
125 
126  // Energy limits
127 
128  if (LowEnergyLimit() < 11.*eV)
129  {
130  G4cout << "G4DNACPA100ElasticModel: low energy limit increased from " <<
131  LowEnergyLimit()/eV << " eV to " << 11 << " eV" << G4endl;
132  SetLowEnergyLimit(11.*eV);
133  }
134 
135  if (HighEnergyLimit() > 255955.*eV)
136  {
137  G4cout << "G4DNACPA100ElasticModel: high energy limit decreased from " <<
138  HighEnergyLimit()/keV << " keV to " << 255.955 << " keV"
139  << G4endl;
140  SetHighEnergyLimit(255955.*eV);
141  }
142 
143  // Reading of data files
144 
145  G4double scaleFactor = 1e-20*m*m;
146 
147  G4String fileElectron("dna/sigma_elastic_e_cpa100");
148 
149  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
150  G4String electron;
151 
152  // *** ELECTRON
153 
154  // For total cross section
155 
156  electron = electronDef->GetParticleName();
157 
158  tableFile[electron] = fileElectron;
159 
160  G4DNACrossSectionDataSet* tableE =
161  new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
162  eV,scaleFactor );
163 
164  /*
165  G4DNACrossSectionDataSet* tableE =
166  new G4DNACrossSectionDataSet(new G4DNACPA100LogLogInterpolation,
167  eV,scaleFactor );
168  */
169 
170  tableE->LoadData(fileElectron);
171 
172  tableData[electron] = tableE;
173 
174  // For final state
175 
176  char *path = getenv("G4LEDATA");
177 
178  if (!path)
179  {
180  G4Exception("G4DNACPA100ElasticModel::Initialise","em0006",
181  FatalException,"G4LEDATA environment variable not set.");
182  return;
183  }
184 
185  std::ostringstream eFullFileName;
186 
187  eFullFileName << path
188  << "/dna/sigmadiff_cumulated_elastic_e_cpa100.dat";
189 
190  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
191 
192  if (!eDiffCrossSection)
193  G4Exception("G4DNACPA100ElasticModel::Initialise","em0003",
194  FatalException,
195  "Missing data file:/dna/sigmadiff_cumulated_elastic_e_cpa100.dat");
196 
197  // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
198  // Added clear for MT
199 
200  eTdummyVec.clear();
201  eVecm.clear();
202  eDiffCrossSectionData.clear();
203 
204  //
205 
206  eTdummyVec.push_back(0.);
207 
208  while(!eDiffCrossSection.eof())
209  {
210  G4double tDummy;
211  G4double eDummy;
212  eDiffCrossSection>>tDummy>>eDummy;
213 
214  // SI : mandatory eVecm initialization
215 
216  if (tDummy != eTdummyVec.back())
217  {
218  eTdummyVec.push_back(tDummy);
219  eVecm[tDummy].push_back(0.);
220  }
221 
222  eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
223 
224  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
225  }
226 
227  // End final state
228 
229 #ifdef UEHARA_VERBOSE
230  if (verboseLevel > 2)
231  G4cout << "Loaded cross section files for CPA100 Elastic model" << G4endl;
232 #endif
233 
234 #ifdef UEHARA_VERBOSE
235  if( verboseLevel>0 )
236  {
237  G4cout << "CPA100 Elastic model is initialized " << G4endl
238  << "Energy range: "
239  << LowEnergyLimit() / eV << " eV - "
240  << HighEnergyLimit() / keV << " keV"
241  << G4endl;
242  }
243 #endif
244 
245  // Initialize water density pointer
246  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()
247  ->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
248 
249  if (isInitialised) { return; }
250  fParticleChangeForGamma = GetParticleChangeForGamma();
251  isInitialised = true;
252 
253 }
254 
255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
256 
257 G4double G4DNACPA100ElasticModel::CrossSectionPerVolume
258 (const G4Material* material,
259  const G4ParticleDefinition* p,
260  G4double ekin,
261  G4double,
262  G4double)
263 {
264 
265 #ifdef UEHARA_VERBOSE
266  if (verboseLevel > 3)
267  G4cout <<
268  "Calling CrossSectionPerVolume() of G4DNACPA100ElasticModel" << G4endl;
269 #endif
270 
271  // Calculate total cross section for model
272 
273  G4double sigma=0;
274 
275  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
276 
277  const G4String& particleName = p->GetParticleName();
278 
279  if (ekin <= HighEnergyLimit() && ekin >= LowEnergyLimit())
280  {
281  //SI : XS must not be zero otherwise sampling of secondaries
282  // method ignored
283 
284  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
285  pos = tableData.find(particleName);
286 
287  if (pos != tableData.end())
288  {
289  G4DNACrossSectionDataSet* table = pos->second;
290  if (table != 0)
291  {
292  sigma = table->FindValue(ekin);
293  //
294  //Dump in non-MT mode
295  //
296  /*
297  G4double minEnergy = 10.481 * eV;
298  G4double maxEnergy = 255955. * eV;
299  G4int nEnergySteps = 1000;
300  G4double energy(minEnergy);
301  G4double
302  stpEnergy(std::pow(maxEnergy/energy,
303  1./static_cast<G4double>(nEnergySteps-1)));
304  G4int step(nEnergySteps);
305  system ("rm -rf elastic-cpa100.out");
306  FILE* myFile=fopen("elastic-cpa100.out","a");
307  while (step>0)
308  {
309  step--;
310  fprintf (myFile,"%16.9le %16.9le\n",
311  energy/eV,
312  table->FindValue(energy)/(1e-20*m*m));
313  energy*=stpEnergy;
314  }
315  fclose (myFile);
316  abort();
317  */
318  //
319  // end of dump
320  //
321  }
322  }
323  else
324  {
325  G4Exception("G4DNACPA100ElasticModel::ComputeCrossSectionPerVolume",
326  "em0002",
327  FatalException,"Model not applicable to particle type.");
328  }
329  }
330 
331 #ifdef UEHARA_VERBOSE
332  if (verboseLevel > 2)
333  {
334  G4cout << "__________________________________" << G4endl;
335  G4cout << "G4DNACPA100ElasticModel - XS INFO START" << G4endl;
336  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
337  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
338  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
339  // G4cout << " - Cross section per water molecule (cm^-1)="
340  // << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
341  G4cout << "G4DNACPA100ElasticModel - XS INFO END" << G4endl;
342  }
343 #endif
344 
345  return sigma*waterDensity;
346 }
347 
348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
349 
350 void G4DNACPA100ElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
351  const G4MaterialCutsCouple* /*couple*/,
352  const G4DynamicParticle* aDynamicElectron,
353  G4double,
354  G4double)
355 {
356 #ifdef UEHARA_VERBOSE
357  if (verboseLevel > 3)
358  G4cout << "Calling SampleSecondaries() of G4DNACPA100ElasticModel" << G4endl;
359 #endif
360 
361  G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
362 
363  G4double cosTheta = RandomizeCosTheta(electronEnergy0);
364  G4double phi = 2. * pi * G4UniformRand();
365 
366  G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
367 
368  //G4ThreeVector xVers = zVers.orthogonal();
369  //G4ThreeVector yVers = zVers.cross(xVers);
370  //G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
371  //G4double yDir = xDir;
372  //xDir *= std::cos(phi);
373  //yDir *= std::sin(phi);
374 
375  // Computation of scattering angles (from Subroutine DIRAN in CPA100)
376 
377  G4double CT1, ST1, CF1, SF1, CT2, ST2, CF2, SF2;
378  G4double sinTheta = std::sqrt (1-cosTheta*cosTheta);
379 
380  CT1=0;
381  ST1=0;
382  CF1=0;
383  SF1=0;
384  CT2=0;
385  ST2=0;
386  CF2=0;
387  SF2=0;
388 
389  CT1 = zVers.z();
390  ST1=std::sqrt(1.-CT1*CT1);
391 
392  if (ST1!=0) CF1 = zVers.x()/ST1; else CF1 = std::cos(2. * pi * G4UniformRand());
393  if (ST1!=0) SF1 = zVers.y()/ST1; else SF1 = std::sqrt(1.-CF1*CF1);
394 
395  G4double A3, A4, A5, A2, A1;
396  A3=0;
397  A4=0;
398  A5=0;
399  A2=0;
400  A1=0;
401 
402  A3 = sinTheta*std::cos(phi);
403  A4 = A3*CT1 + ST1*cosTheta;
404  A5 = sinTheta * std::sin(phi);
405  A2 = A4 * SF1 + A5 * CF1;
406  A1 = A4 * CF1 - A5 * SF1;
407 
408  CT2 = CT1*cosTheta - ST1*A3;
409  ST2 = std::sqrt(1.-CT2*CT2);
410 
411  if (ST2==0) ST2=1E-6;
412  CF2 = A1/ST2;
413  SF2 = A2/ST2;
414 
415  /*
416  G4cout << "CT1=" << CT1 << G4endl;
417  G4cout << "ST1=" << ST1 << G4endl;
418  G4cout << "CF1=" << CF1 << G4endl;
419  G4cout << "SF1=" << SF1 << G4endl;
420  G4cout << "cosTheta=" << cosTheta << G4endl;
421  G4cout << "sinTheta=" << sinTheta << G4endl;
422  G4cout << "cosPhi=" << std::cos(phi) << G4endl;
423  G4cout << "sinPhi=" << std::sin(phi) << G4endl;
424  G4cout << "CT2=" << CT2 << G4endl;
425  G4cout << "ST2=" << ST2 << G4endl;
426  G4cout << "CF2=" << CF2 << G4endl;
427  G4cout << "SF2=" << SF2 << G4endl;
428  */
429 
430  G4ThreeVector zPrimeVers(ST2*CF2,ST2*SF2,CT2);
431 
432  //
433 
434  fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
435 
436  if (!statCode)
437 
438  fParticleChangeForGamma->SetProposedKineticEnergy
439  (electronEnergy0-1.214E-4*(1.-cosTheta)*electronEnergy0);
440 
441  else fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
442 
443  //
444 
445  fParticleChangeForGamma->ProposeLocalEnergyDeposit(1.214E-4*(1.-cosTheta)*electronEnergy0);
446 
447 }
448 
449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
450 
451 G4double G4DNACPA100ElasticModel::Theta
452  (G4ParticleDefinition *, G4double k, G4double integrDiff)
453 {
454 
455  G4double theta = 0.;
456  G4double valueT1 = 0;
457  G4double valueT2 = 0;
458  G4double valueE21 = 0;
459  G4double valueE22 = 0;
460  G4double valueE12 = 0;
461  G4double valueE11 = 0;
462  G4double xs11 = 0;
463  G4double xs12 = 0;
464  G4double xs21 = 0;
465  G4double xs22 = 0;
466 
467  // Protection against out of boundary access
468  if (k==eTdummyVec.back()) k=k*(1.-1e-12);
469  //
470 
471  std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
472  std::vector<G4double>::iterator t1 = t2-1;
473 
474  std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(),
475  integrDiff);
476  std::vector<G4double>::iterator e11 = e12-1;
477 
478  std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(),
479  integrDiff);
480  std::vector<G4double>::iterator e21 = e22-1;
481 
482  valueT1 =*t1;
483  valueT2 =*t2;
484  valueE21 =*e21;
485  valueE22 =*e22;
486  valueE12 =*e12;
487  valueE11 =*e11;
488 
489  xs11 = eDiffCrossSectionData[valueT1][valueE11];
490  xs12 = eDiffCrossSectionData[valueT1][valueE12];
491  xs21 = eDiffCrossSectionData[valueT2][valueE21];
492  xs22 = eDiffCrossSectionData[valueT2][valueE22];
493 
494  //TEST CPA100
495  //if(k==valueT1) xs22 = eDiffCrossSectionData[valueT1][valueE12];
496 
497  if (xs11==0 && xs12==0 && xs21==0 && xs22==0) return (0.);
498 
499  theta = QuadInterpolator(
500  valueE11, valueE12,
501  valueE21, valueE22,
502  xs11, xs12,
503  xs21, xs22,
504  valueT1, valueT2,
505  k, integrDiff);
506 
507  return theta;
508 
509  //TEST CPA100
510  //return xs22;
511 }
512 
513 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
514 
515 G4double G4DNACPA100ElasticModel::LinLogInterpolate(G4double e1,
516  G4double e2,
517  G4double e,
518  G4double xs1,
519  G4double xs2)
520 {
521  G4double d1 = std::log(xs1);
522  G4double d2 = std::log(xs2);
523  G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
524  return value;
525 }
526 
527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
528 
529 G4double G4DNACPA100ElasticModel::LinLinInterpolate(G4double e1,
530  G4double e2,
531  G4double e,
532  G4double xs1,
533  G4double xs2)
534 {
535  G4double d1 = xs1;
536  G4double d2 = xs2;
537  G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
538  return value;
539 }
540 
541 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
542 
543 G4double G4DNACPA100ElasticModel::LogLogInterpolate(G4double e1,
544  G4double e2,
545  G4double e,
546  G4double xs1,
547  G4double xs2)
548 {
549  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
550  G4double b = std::log10(xs2) - a*std::log10(e2);
551  G4double sigma = a*std::log10(e) + b;
552  G4double value = (std::pow(10.,sigma));
553  return value;
554 }
555 
556 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
557 
558 
559 G4double G4DNACPA100ElasticModel::QuadInterpolator(G4double e11, G4double e12,
560  G4double e21, G4double e22,
561  G4double xs11, G4double xs12,
562  G4double xs21, G4double xs22,
563  G4double t1, G4double t2,
564  G4double t, G4double e)
565 {
566  // Log-Log
567 /*
568  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
569  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
570  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
571 
572 
573  // Lin-Log
574  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
575  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
576  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
577 */
578 
579  // Lin-Lin
580  G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
581  G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
582  G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
583 
584  return value;
585 }
586 
587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
588 
589 G4double G4DNACPA100ElasticModel::RandomizeCosTheta(G4double k)
590 {
591 
592  G4double integrdiff=0; // PROBABILITY between 0 and 1.
593  G4double uniformRand=G4UniformRand();
594  integrdiff = uniformRand;
595 
596  G4double cosTheta=0.;
597 
598  // 1 - COS THETA is read from the data file
599  cosTheta = 1 - Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
600 
601  //
602  //
603  //Dump
604  //
605  //G4cout << "theta=" << theta << G4endl;
606  //G4cout << "cos theta=" << std::cos(theta*pi/180) << G4endl;
607  //G4cout << "sin theta=" << std::sin(theta*pi/180) << G4endl;
608  //G4cout << "acos(cos theta)=" << std::acos(cosTheta) << G4endl;
609  //G4cout << "cos theta="<< cosTheta << G4endl;
610  //G4cout << "1 - cos theta="<< 1. - cosTheta << G4endl;
611  //G4cout << "sin theta=" << std::sqrt(1-cosTheta*cosTheta) << G4endl;
612  //
613  /*
614  G4double minProb = 0; // we scan probability between 0 and one
615  G4double maxProb = 1;
616  G4int nProbSteps = 100;
617  G4double prob(minProb);
618  G4double stepProb((maxProb-minProb)/static_cast<G4double>(nProbSteps));
619  G4int step(nProbSteps);
620  system ("rm -rf elastic-cumul-cpa100-100keV.out");
621  FILE* myFile=fopen("elastic-cumul-cpa100-100keV.out","a");
622  while (step>=0)
623  {
624  step--;
625  fprintf (myFile,"%16.9le %16.9le\n",
626  prob,
627  Theta(G4Electron::ElectronDefinition(),100000,prob)); // SELECT NRJ IN eV !!!
628  prob=prob+stepProb;
629  }
630  fclose (myFile);
631  abort();
632  */
633  //
634  // end of dump
635  //
636 
637  return cosTheta;
638 }