ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNACPA100ExcitationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNACPA100ExcitationModel.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 excitation 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 "G4DNACPA100ExcitationModel.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4PhysicalConstants.hh"
44 #include "G4DNAChemistryManager.hh"
45 #include "G4DNAMolecularMaterial.hh"
46 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 using namespace std;
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52 
53 G4DNACPA100ExcitationModel::G4DNACPA100ExcitationModel(const G4ParticleDefinition*,
54  const G4String& nam)
55 :G4VEmModel(nam),isInitialised(false)
56 {
57  fpMolWaterDensity = 0;
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  if( verboseLevel>0 )
71  {
72  G4cout << "CPA100 excitation model is constructed " << G4endl;
73  }
74  fParticleChangeForGamma = 0;
75 
76  // Selection of stationary mode
77 
78  statCode = false;
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82 
83 G4DNACPA100ExcitationModel::~G4DNACPA100ExcitationModel()
84 {
85  // Cross section
86 
87  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
88  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
89  {
90  G4DNACrossSectionDataSet* table = pos->second;
91  delete table;
92  }
93 
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 
98 void G4DNACPA100ExcitationModel::Initialise(const G4ParticleDefinition* particle,
99  const G4DataVector& /*cuts*/)
100 {
101 
102  if (verboseLevel > 3)
103  G4cout << "Calling G4DNACPA100ExcitationModel::Initialise()" << G4endl;
104 
105  G4String fileElectron("dna/sigma_excitation_e_cpa100");
106 
107  G4double scaleFactor = 1.e-20 *m*m;
108 
109  // *** ELECTRON
110 
111  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
112  G4String electron;
113  electron = electronDef->GetParticleName();
114 
115  tableFile[electron] = fileElectron;
116 
117  // Cross section
118 
119  G4DNACrossSectionDataSet* tableE
120  = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV, scaleFactor );
121 
122  /*
123  G4DNACrossSectionDataSet* tableE =
124  new G4DNACrossSectionDataSet(new G4DNACPA100LogLogInterpolation, eV, scaleFactor );
125  */
126 
127  tableE->LoadData(fileElectron);
128 
129  tableData[electron] = tableE;
130 
131  //
132 
133  if( verboseLevel>0 )
134  {
135  G4cout << "CPA100 excitation model is initialized " << G4endl
136  << "Energy range: "
137  << LowEnergyLimit() / eV << " eV - "
138  << HighEnergyLimit() / keV << " keV for "
139  << particle->GetParticleName()
140  << G4endl;
141  }
142 
143  // Initialize water density pointer
144  fpMolWaterDensity =
145  G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
146 
147  if (isInitialised) return;
148  fParticleChangeForGamma = GetParticleChangeForGamma();
149  isInitialised = true;
150 }
151 
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
153 
154 G4double G4DNACPA100ExcitationModel::CrossSectionPerVolume(const G4Material* material,
155  const G4ParticleDefinition* particleDefinition,
156  G4double ekin,
157  G4double,
158  G4double)
159 {
160 
161  if (verboseLevel > 3)
162  G4cout << "Calling CrossSectionPerVolume() of G4DNACPA100ExcitationModel" << G4endl;
163 
164  if (particleDefinition != G4Electron::ElectronDefinition()) return 0;
165 
166  // Calculate total cross section for model
167 
168  G4double sigma=0;
169 
170  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
171 
172  const G4String& particleName = particleDefinition->GetParticleName();
173 
174  if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit())
175  {
176  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
177  pos = tableData.find(particleName);
178 
179  if (pos != tableData.end())
180  {
181  G4DNACrossSectionDataSet* table = pos->second;
182  if (table != 0)
183  {
184  sigma = table->FindValue(ekin);
185  }
186  }
187  else
188  {
189  G4Exception("G4DNACPA100ExcitationModel::CrossSectionPerVolume","em0002",
190  FatalException,"Model not applicable to particle type.");
191  }
192  }
193 
194  if (verboseLevel > 2)
195  {
196  G4cout << "__________________________________" << G4endl;
197  G4cout << "G4DNACPA100ExcitationModel - XS INFO START" << G4endl;
198  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
199  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
200  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
201  // G4cout << " - Cross section per water molecule (cm^-1)="
202  // << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
203  G4cout << "G4DNACPA100ExcitationModel - XS INFO END" << G4endl;
204  }
205 
206  return sigma*waterDensity;
207 
208 }
209 
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211 
212 void G4DNACPA100ExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* ,
213  const G4MaterialCutsCouple*,
214  const G4DynamicParticle* aDynamicParticle,
215  G4double,
216  G4double)
217 {
218 
219  if (verboseLevel > 3)
220  G4cout << "Calling SampleSecondaries() of G4DNACPA100ExcitationModel" << G4endl;
221 
222  G4double k = aDynamicParticle->GetKineticEnergy();
223 
224  const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
225 
226  G4int level = RandomSelect(k,particleName);
227  G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
228  G4double newEnergy = k - excitationEnergy;
229 
230  if (newEnergy > 0)
231  {
232  // fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
233 
234  // We take into account direction change as described page 87 (II.92) in thesis by S. Edel
235 
236  G4double cosTheta =
237 
238  (excitationEnergy/k) / (1. + (k/(2*electron_mass_c2))*(1.-excitationEnergy/k) );
239 
240  cosTheta = std::sqrt(1.-cosTheta);
241 
242  G4double phi = 2. * pi * G4UniformRand();
243 
244  G4ThreeVector zVers = aDynamicParticle->GetMomentumDirection();
245 
246  //G4ThreeVector xVers = zVers.orthogonal();
247  //G4ThreeVector yVers = zVers.cross(xVers);
248  //G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
249  //G4double yDir = xDir;
250  //xDir *= std::cos(phi);
251  //yDir *= std::sin(phi);
252  // G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
253 
254  // Computation of scattering angles (from Subroutine DIRAN in CPA100)
255 
256  G4double CT1, ST1, CF1, SF1, CT2, ST2, CF2, SF2;
257  G4double sinTheta = std::sqrt (1-cosTheta*cosTheta);
258 
259  CT1=0;
260  ST1=0;
261  CF1=0;
262  SF1=0;
263  CT2=0;
264  ST2=0;
265  CF2=0;
266  SF2=0;
267 
268  CT1 = zVers.z();
269  ST1=std::sqrt(1.-CT1*CT1);
270 
271  if (ST1!=0) CF1 = zVers.x()/ST1; else CF1 = std::cos(2. * pi * G4UniformRand());
272  if (ST1!=0) SF1 = zVers.y()/ST1; else SF1 = std::sqrt(1.-CF1*CF1);
273 
274  G4double A3, A4, A5, A2, A1;
275  A3=0;
276  A4=0;
277  A5=0;
278  A2=0;
279  A1=0;
280 
281  A3 = sinTheta*std::cos(phi);
282  A4 = A3*CT1 + ST1*cosTheta;
283  A5 = sinTheta * std::sin(phi);
284  A2 = A4 * SF1 + A5 * CF1;
285  A1 = A4 * CF1 - A5 * SF1;
286 
287  CT2 = CT1*cosTheta - ST1*A3;
288  ST2 = std::sqrt(1.-CT2*CT2);
289 
290  if (ST2==0) ST2=1E-6;
291  CF2 = A1/ST2;
292  SF2 = A2/ST2;
293 
294  /*
295  G4cout << "CT1=" << CT1 << G4endl;
296  G4cout << "ST1=" << ST1 << G4endl;
297  G4cout << "CF1=" << CF1 << G4endl;
298  G4cout << "SF1=" << SF1 << G4endl;
299  G4cout << "cosTheta=" << cosTheta << G4endl;
300  G4cout << "sinTheta=" << sinTheta << G4endl;
301  G4cout << "cosPhi=" << std::cos(phi) << G4endl;
302  G4cout << "sinPhi=" << std::sin(phi) << G4endl;
303  G4cout << "CT2=" << CT2 << G4endl;
304  G4cout << "ST2=" << ST2 << G4endl;
305  G4cout << "CF2=" << CF2 << G4endl;
306  G4cout << "SF2=" << SF2 << G4endl;
307  */
308 
309  G4ThreeVector zPrimeVers(ST2*CF2,ST2*SF2,CT2);
310 
311  //
312 
313  fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
314 
315  //
316 
317  if (!statCode) fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
318  else fParticleChangeForGamma->SetProposedKineticEnergy(k);
319 
320  fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
321  }
322 
323  // Chemistry
324 
325  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
326  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
327  level,
328  theIncomingTrack);
329 }
330 
331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
332 
333 G4int G4DNACPA100ExcitationModel::RandomSelect(G4double k, const G4String& particle)
334 {
335  G4int level = 0;
336 
337  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
338  pos = tableData.find(particle);
339 
340  if (pos != tableData.end())
341  {
342  G4DNACrossSectionDataSet* table = pos->second;
343 
344  if (table != 0)
345  {
346  G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
347  const size_t n(table->NumberOfComponents());
348  size_t i(n);
349  G4double value = 0.;
350 
351  //Verification
352  /*
353  G4double tmp=10.481*eV;
354  G4cout << table->GetComponent(0)->FindValue(tmp)/(1e-20*m*m) << G4endl;
355  G4cout << table->GetComponent(1)->FindValue(tmp)/(1e-20*m*m) << G4endl;
356  G4cout << table->GetComponent(2)->FindValue(tmp)/(1e-20*m*m) << G4endl;
357  G4cout << table->GetComponent(3)->FindValue(tmp)/(1e-20*m*m) << G4endl;
358  G4cout << table->GetComponent(4)->FindValue(tmp)/(1e-20*m*m) << G4endl;
359  G4cout <<
360  table->GetComponent(0)->FindValue(tmp)/(1e-20*m*m) +
361  table->GetComponent(1)->FindValue(tmp)/(1e-20*m*m) +
362  table->GetComponent(2)->FindValue(tmp)/(1e-20*m*m) +
363  table->GetComponent(3)->FindValue(tmp)/(1e-20*m*m) +
364  table->GetComponent(4)->FindValue(tmp)/(1e-20*m*m)
365  << G4endl;
366  abort();
367  */
368  //
369  //Dump
370  //
371  /*
372  G4double minEnergy = 10.481 * eV;
373  G4double maxEnergy = 255955. * eV;
374  G4int nEnergySteps = 1000;
375  G4double energy(minEnergy);
376  G4double stpEnergy(std::pow(maxEnergy/energy, 1./static_cast<G4double>(nEnergySteps-1)));
377  G4int step(nEnergySteps);
378  system ("rm -rf excitation-cap100.out");
379  FILE* myFile=fopen("excitation-cpa100.out","a");
380  while (step>0)
381  {
382  step--;
383  fprintf (myFile,"%16.9le %16.9le %16.9le %16.9le %16.9le %16.9le %16.9le \n",
384  energy/eV,
385  table->GetComponent(0)->FindValue(energy)/(1e-20*m*m),
386  table->GetComponent(1)->FindValue(energy)/(1e-20*m*m),
387  table->GetComponent(2)->FindValue(energy)/(1e-20*m*m),
388  table->GetComponent(3)->FindValue(energy)/(1e-20*m*m),
389  table->GetComponent(4)->FindValue(energy)/(1e-20*m*m),
390  table->GetComponent(0)->FindValue(energy)/(1e-20*m*m)+
391  table->GetComponent(1)->FindValue(energy)/(1e-20*m*m)+
392  table->GetComponent(2)->FindValue(energy)/(1e-20*m*m)+
393  table->GetComponent(3)->FindValue(energy)/(1e-20*m*m)+
394  table->GetComponent(4)->FindValue(energy)/(1e-20*m*m)
395  );
396  energy*=stpEnergy;
397  }
398  fclose (myFile);
399  abort();
400  */
401  //
402  // end of dump
403  //
404 
405  while (i>0)
406  {
407  i--;
408  valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
409  value += valuesBuffer[i];
410  }
411 
412  value *= G4UniformRand();
413 
414  i = n;
415 
416  while (i > 0)
417  {
418  i--;
419 
420  if (valuesBuffer[i] > value)
421  {
422  delete[] valuesBuffer;
423  return i;
424  }
425  value -= valuesBuffer[i];
426  }
427 
428  if (valuesBuffer) delete[] valuesBuffer;
429 
430  }
431  }
432  else
433  {
434  G4Exception("G4DNACPA100ExcitationModel::RandomSelect","em0002",
435  FatalException,"Model not applicable to particle type.");
436  }
437  return level;
438 }