ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DNAPTBExcitationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DNAPTBExcitationModel.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 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // Models come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459-479 (2017)
29 //
30 
31 #include "G4DNAPTBExcitationModel.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMolecularMaterial.hh"
35 
36 G4DNAPTBExcitationModel::G4DNAPTBExcitationModel(const G4String& applyToMaterial, const G4ParticleDefinition*,
37  const G4String& nam)
38  : G4VDNAModel(nam, applyToMaterial)
39 {
40  verboseLevel= 0;
41  // Verbosity scale:
42  // 0 = nothing
43  // 1 = warning for energy non-conservation
44  // 2 = details of energy budget
45  // 3 = calculation of cross sections, file openings, sampling of atoms
46  // 4 = entering in methods
47 
48  // initialisation of mean energy loss for each material
49  tableMeanEnergyPTB["THF"] = 8.01*eV;
50  tableMeanEnergyPTB["PY"] = 7.61*eV;
51  tableMeanEnergyPTB["PU"] = 7.61*eV;
52  tableMeanEnergyPTB["TMP"] = 8.01*eV;
53 
54  if( verboseLevel>0 )
55  {
56  G4cout << "PTB excitation model is constructed " << G4endl;
57  }
58 }
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 
62 G4DNAPTBExcitationModel::~G4DNAPTBExcitationModel()
63 {
64 
65 }
66 
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 
69 void G4DNAPTBExcitationModel::Initialise(const G4ParticleDefinition* particle,
70  const G4DataVector& /*cuts*/, G4ParticleChangeForGamma*)
71 {
72  if (verboseLevel > 3)
73  G4cout << "Calling G4DNAPTBExcitationModel::Initialise()" << G4endl;
74 
75  G4double scaleFactor = 1e-16*cm*cm;
76  G4double scaleFactorBorn = (1.e-22 / 3.343) * m*m;
77 
78  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
79 
80  //*******************************************************
81  // Cross section data
82  //*******************************************************
83 
84  if(particle == electronDef)
85  {
86  G4String particleName = particle->GetParticleName();
87 
88  AddCrossSectionData("THF",
89  particleName,
90  "dna/sigma_excitation_e-_PTB_THF",
91  scaleFactor);
92  SetLowELimit("THF", particleName, 9.*eV);
93  SetHighELimit("THF", particleName, 1.*keV);
94 
95  AddCrossSectionData("PY",
96  particleName,
97  "dna/sigma_excitation_e-_PTB_PY",
98  scaleFactor);
99  SetLowELimit("PY", particleName, 9.*eV);
100  SetHighELimit("PY", particleName, 1.*keV);
101 
102  AddCrossSectionData("PU",
103  particleName,
104  "dna/sigma_excitation_e-_PTB_PU",
105  scaleFactor);
106  SetLowELimit("PU", particleName, 9.*eV);
107  SetHighELimit("PU", particleName, 1.*keV);
108 
109  AddCrossSectionData("TMP",
110  particleName,
111  "dna/sigma_excitation_e-_PTB_TMP",
112  scaleFactor);
113  SetLowELimit("TMP", particleName, 9.*eV);
114  SetHighELimit("TMP", particleName, 1.*keV);
115 
116  AddCrossSectionData("G4_WATER",
117  particleName,
118  "dna/sigma_excitation_e_born",
119  scaleFactorBorn);
120  SetLowELimit("G4_WATER", particleName, 9.*eV);
121  SetHighELimit("G4_WATER", particleName, 1.*keV);
122 
123  // DNA materials
124  //
125  AddCrossSectionData("backbone_THF",
126  particleName,
127  "dna/sigma_excitation_e-_PTB_THF",
128  scaleFactor*33./30);
129  SetLowELimit("backbone_THF", particleName, 9.*eV);
130  SetHighELimit("backbone_THF", particleName, 1.*keV);
131 
132  AddCrossSectionData("cytosine_PY",
133  particleName,
134  "dna/sigma_excitation_e-_PTB_PY",
135  scaleFactor*42./30);
136  SetLowELimit("cytosine_PY", particleName, 9.*eV);
137  SetHighELimit("cytosine_PY", particleName, 1.*keV);
138 
139  AddCrossSectionData("thymine_PY",
140  particleName,
141  "dna/sigma_excitation_e-_PTB_PY",
142  scaleFactor*48./30);
143  SetLowELimit("thymine_PY", particleName, 9.*eV);
144  SetHighELimit("thymine_PY", particleName, 1.*keV);
145 
146  AddCrossSectionData("adenine_PU",
147  particleName,
148  "dna/sigma_excitation_e-_PTB_PU",
149  scaleFactor*50./44);
150  SetLowELimit("adenine_PU", particleName, 9.*eV);
151  SetHighELimit("adenine_PU", particleName, 1.*keV);
152 
153  AddCrossSectionData("guanine_PU",
154  particleName,
155  "dna/sigma_excitation_e-_PTB_PU",
156  scaleFactor*56./44);
157  SetLowELimit("guanine_PU", particleName, 9.*eV);
158  SetHighELimit("guanine_PU", particleName, 1.*keV);
159 
160  AddCrossSectionData("backbone_TMP",
161  particleName,
162  "dna/sigma_excitation_e-_PTB_TMP",
163  scaleFactor*33./50);
164  SetLowELimit("backbone_TMP", particleName, 9.*eV);
165  SetHighELimit("backbone_TMP", particleName, 1.*keV);
166  }
167 
168  //*******************************************************
169  // Load data
170  //*******************************************************
171 
172  LoadCrossSectionData(particle->GetParticleName() );
173 
174  //*******************************************************
175  // Verbose
176  //*******************************************************
177 
178  if( verboseLevel>0 )
179  {
180  G4cout << "PTB excitation model is initialized " << G4endl;
181  }
182 }
183 
184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 
186 G4double G4DNAPTBExcitationModel::CrossSectionPerVolume(const G4Material* /*material*/,
187  const G4String& materialName,
188  const G4ParticleDefinition* particleDefinition,
189  G4double ekin,
190  G4double /*emin*/,
191  G4double /*emax*/)
192 {
193  if (verboseLevel > 3)
194  G4cout << "Calling CrossSectionPerVolume() of G4DNAPTBExcitationModel" << G4endl;
195 
196  // Get the name of the current particle
197  G4String particleName = particleDefinition->GetParticleName();
198 
199  // initialise variables
200  G4double lowLim = 0;
201  G4double highLim = 0;
202  G4double sigma=0;
203 
204  // Get the low energy limit for the current particle
205  lowLim = GetLowELimit(materialName, particleName);
206 
207  // Get the high energy limit for the current particle
208  highLim = GetHighELimit(materialName, particleName);
209 
210  // Check that we are in the correct energy range
211  if (ekin >= lowLim && ekin < highLim)
212  {
213  // Get the map with all the data tables
214  TableMapData* tableData = GetTableData();
215 
216  // Retrieve the cross section value
217  sigma = (*tableData)[materialName][particleName]->FindValue(ekin);
218 
219  if (verboseLevel > 2)
220  {
221  G4cout << "__________________________________" << G4endl;
222  G4cout << "°°° G4DNAPTBExcitationModel - XS INFO START" << G4endl;
223  G4cout << "°°° Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
224  G4cout << "°°° Cross section per "<< materialName <<" molecule (cm^2)=" << sigma/cm/cm << G4endl;
225  G4cout << "°°° G4DNAPTBExcitationModel - XS INFO END" << G4endl;
226  }
227 
228  }
229 
230  // Return the cross section value
231  return sigma;
232 }
233 
234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
235 
236 void G4DNAPTBExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
237  const G4MaterialCutsCouple* /*couple*/,
238  const G4String& materialName,
239  const G4DynamicParticle* aDynamicParticle,
240  G4ParticleChangeForGamma* particleChangeForGamma,
241  G4double /*tmin*/,
242  G4double /*tmax*/)
243 {
244  if (verboseLevel > 3)
245  G4cout << "Calling SampleSecondaries() of G4DNAPTBExcitationModel" << G4endl;
246 
247  // Get the incident particle kinetic energy
248  G4double k = aDynamicParticle->GetKineticEnergy();
249 
250  if(materialName!="G4_WATER")
251  {
252  // Retrieve the excitation energy for the current material
253  G4double excitationEnergy = tableMeanEnergyPTB[materialName];
254 
255  // Calculate the new energy of the particle
256  G4double newEnergy = k - excitationEnergy;
257 
258  // Check that the new energy is above zero before applying it the particle.
259  // Otherwise, do nothing.
260  if (newEnergy > 0)
261  {
262  particleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
263  particleChangeForGamma->SetProposedKineticEnergy(newEnergy);
264  particleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
265  }
266  }
267  else
268  {
269  const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
270 
271  G4int level = RandomSelectShell(k,particleName, materialName);
272  G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
273  G4double newEnergy = k - excitationEnergy;
274 
275  if (newEnergy > 0)
276  {
277  particleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
278  particleChangeForGamma->SetProposedKineticEnergy(newEnergy);
279  particleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
280  }
281 
282  const G4Track * theIncomingTrack = particleChangeForGamma->GetCurrentTrack();
283  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
284  level,
285  theIncomingTrack);
286  }
287 }