| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Modified by Z. Francis, S. Incerti to handl 27 // Modified by Z. Francis, S. Incerti to handle HZE
28 // && inverse rudd function sampling 26-10-201 28 // && inverse rudd function sampling 26-10-2010
29 // <<
30 // Rewitten by V.Ivanchenko 21.05.2023 <<
31 // <<
32 29
33 #include "G4EmCorrections.hh" <<
34 #include "G4DNARuddIonisationExtendedModel.hh" 30 #include "G4DNARuddIonisationExtendedModel.hh"
35 #include "G4PhysicalConstants.hh" 31 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh"
37 #include "G4UAtomicDeexcitation.hh" 33 #include "G4UAtomicDeexcitation.hh"
38 #include "G4LossTableManager.hh" 34 #include "G4LossTableManager.hh"
39 #include "G4NistManager.hh" <<
40 #include "G4DNAChemistryManager.hh" 35 #include "G4DNAChemistryManager.hh"
41 #include "G4DNAMolecularMaterial.hh" 36 #include "G4DNAMolecularMaterial.hh"
42 37
43 #include "G4IonTable.hh" 38 #include "G4IonTable.hh"
44 #include "G4DNARuddAngle.hh" 39 #include "G4DNARuddAngle.hh"
45 #include "G4DeltaAngle.hh" 40 #include "G4DeltaAngle.hh"
46 #include "G4Exp.hh" 41 #include "G4Exp.hh"
47 #include "G4Log.hh" <<
48 #include "G4Pow.hh" <<
49 #include "G4Alpha.hh" <<
50 #include "G4Proton.hh" <<
51 42
52 //....oooOO0OOooo........oooOO0OOooo........oo 43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 44
54 G4DNACrossSectionDataSet* G4DNARuddIonisationE << 45 using namespace std;
55 G4DNACrossSectionDataSet* G4DNARuddIonisationE <<
56 G4DNACrossSectionDataSet* G4DNARuddIonisationE <<
57 const std::vector<G4double>* G4DNARuddIonisati <<
58 <<
59 namespace <<
60 { <<
61 const G4double scaleFactor = CLHEP::m*CLHEP: <<
62 const G4double tolerance = 1*CLHEP::eV; <<
63 const G4double Ry = 13.6*CLHEP::eV; <<
64 <<
65 // Following values provided by M. Dingfelde <<
66 const G4double Bj[5] = {12.60*CLHEP::eV, 14. <<
67 32.20*CLHEP::eV, 539 <<
68 } <<
69 46
70 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 48
72 G4DNARuddIonisationExtendedModel::G4DNARuddIon 49 G4DNARuddIonisationExtendedModel::G4DNARuddIonisationExtendedModel(const G4ParticleDefinition*,
73 50 const G4String& nam)
74 : G4VEmModel(nam) << 51 :G4VEmModel(nam),isInitialised(false)
75 { 52 {
76 fEmCorrections = G4LossTableManager::Instanc << 53 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
77 fGpow = G4Pow::GetInstance(); << 54 fpWaterDensity = 0;
78 fLowestEnergy = 100*CLHEP::eV; << 55
79 fLimitEnergy = 1*CLHEP::keV; << 56 slaterEffectiveCharge[0]=0.;
>> 57 slaterEffectiveCharge[1]=0.;
>> 58 slaterEffectiveCharge[2]=0.;
>> 59 sCoefficient[0]=0.;
>> 60 sCoefficient[1]=0.;
>> 61 sCoefficient[2]=0.;
>> 62
>> 63 lowEnergyLimitForA[1] = 0 * eV;
>> 64 lowEnergyLimitForA[2] = 0 * eV;
>> 65 lowEnergyLimitForA[3] = 0 * eV;
>> 66 lowEnergyLimitOfModelForA[1] = 100 * eV;
>> 67 lowEnergyLimitOfModelForA[4] = 1 * keV;
>> 68 lowEnergyLimitOfModelForA[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma
>> 69 killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1];
>> 70 killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4];
>> 71 killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5];
>> 72
>> 73 verboseLevel= 0;
>> 74 // Verbosity scale:
>> 75 // 0 = nothing
>> 76 // 1 = warning for energy non-conservation
>> 77 // 2 = details of energy budget
>> 78 // 3 = calculation of cross sections, file openings, sampling of atoms
>> 79 // 4 = entering in methods
>> 80
>> 81 if( verboseLevel>0 )
>> 82 {
>> 83 G4cout << "Rudd ionisation model is constructed " << G4endl;
>> 84 }
80 85
81 // Mark this model as "applicable" for atomi << 86 // Define default angular generator
82 SetDeexcitationFlag(true); << 87 SetAngularDistribution(new G4DNARuddAngle());
83 88
84 // Define default angular generator << 89 // Mark this model as "applicable" for atomic deexcitation
85 SetAngularDistribution(new G4DNARuddAngle()) << 90 SetDeexcitationFlag(true);
>> 91 fAtomDeexcitation = 0;
>> 92 fParticleChangeForGamma = 0;
86 93
87 if (nullptr == xshelium) { LoadData(); } << 94 // Selection of stationary mode
>> 95
>> 96 statCode = false;
88 } 97 }
89 98
90 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 100
92 G4DNARuddIonisationExtendedModel::~G4DNARuddIo 101 G4DNARuddIonisationExtendedModel::~G4DNARuddIonisationExtendedModel()
93 { 102 {
94 if(isFirst) { << 103 // Cross section
95 for(auto & i : xsdata) { delete i; } << 104
96 } << 105 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 106 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
>> 107 {
>> 108 G4DNACrossSectionDataSet* table = pos->second;
>> 109 delete table;
>> 110 }
>> 111
>> 112 // The following removal is forbidden G4VEnergyLossModel takes care of deletion
>> 113 // however coverity will signal this as an error
>> 114 // if (fAtomDeexcitation) {delete fAtomDeexcitation;}
>> 115
97 } 116 }
98 117
99 //....oooOO0OOooo........oooOO0OOooo........oo 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 119
101 void G4DNARuddIonisationExtendedModel::LoadDat << 120 void G4DNARuddIonisationExtendedModel::Initialise(const G4ParticleDefinition* particle,
>> 121 const G4DataVector& /*cuts*/)
102 { 122 {
103 // initialisation of static data once << 123 if (verboseLevel > 3)
104 isFirst = true; << 124 G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl;
105 G4String filename("dna/sigma_ionisation_h_ru << 125
106 xsdata[0] = new G4DNACrossSectionDataSet(new << 126 // Energy limits
107 xsdata[0]->LoadData(filename); << 127
108 << 128 G4String fileProton("dna/sigma_ionisation_p_rudd");
109 filename = "dna/sigma_ionisation_p_rudd"; << 129 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
110 xsdata[1] = new G4DNACrossSectionDataSet(new << 130 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
111 xsdata[1]->LoadData(filename); << 131 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
112 << 132 G4String fileHelium("dna/sigma_ionisation_he_rudd");
113 filename = "dna/sigma_ionisation_alphapluspl << 133 G4String fileLithium("dna/sigma_ionisation_li_rudd");
114 xsdata[2] = new G4DNACrossSectionDataSet(new << 134 G4String fileBeryllium("dna/sigma_ionisation_be_rudd");
115 xsdata[2]->LoadData(filename); << 135 G4String fileBoron("dna/sigma_ionisation_b_rudd");
116 << 136 G4String fileCarbon("dna/sigma_ionisation_c_rudd");
117 filename = "dna/sigma_ionisation_li_rudd"; << 137 G4String fileNitrogen("dna/sigma_ionisation_n_rudd");
118 xsdata[3] = new G4DNACrossSectionDataSet(new << 138 G4String fileOxygen("dna/sigma_ionisation_o_rudd");
119 xsdata[3]->LoadData(filename); << 139 G4String fileSilicon("dna/sigma_ionisation_si_rudd");
120 << 140 G4String fileIron("dna/sigma_ionisation_fe_rudd");
121 filename = "dna/sigma_ionisation_be_rudd"; << 141
122 xsdata[4] = new G4DNACrossSectionDataSet(new << 142 G4DNAGenericIonsManager *instance;
123 xsdata[4]->LoadData(filename); << 143 instance = G4DNAGenericIonsManager::Instance();
124 << 144 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
125 filename = "dna/sigma_ionisation_b_rudd"; << 145 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
126 xsdata[5] = new G4DNACrossSectionDataSet(new << 146 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
127 xsdata[5]->LoadData(filename); << 147 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
128 << 148 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
129 filename = "dna/sigma_ionisation_c_rudd"; << 149
130 xsdata[6] = new G4DNACrossSectionDataSet(new << 150 //G4ParticleDefinition* carbonDef = instance->GetIon("carbon");
131 xsdata[6]->LoadData(filename); << 151 //G4ParticleDefinition* nitrogenDef = instance->GetIon("nitrogen");
132 << 152 //G4ParticleDefinition* oxygenDef = instance->GetIon("oxygen");
133 filename = "dna/sigma_ionisation_n_rudd"; << 153 //G4ParticleDefinition* siliconDef = instance->GetIon("silicon");
134 xsdata[7] = new G4DNACrossSectionDataSet(new << 154 //G4ParticleDefinition* ironDef = instance->GetIon("iron");
135 xsdata[7]->LoadData(filename); << 155 G4ParticleDefinition* lithiumDef = G4IonTable::GetIonTable()->GetIon(3,7);
136 << 156 G4ParticleDefinition* berylliumDef = G4IonTable::GetIonTable()->GetIon(4,9);
137 filename = "dna/sigma_ionisation_o_rudd"; << 157 G4ParticleDefinition* boronDef = G4IonTable::GetIonTable()->GetIon(5,11);
138 xsdata[8] = new G4DNACrossSectionDataSet(new << 158 G4ParticleDefinition* carbonDef = G4IonTable::GetIonTable()->GetIon(6,12);
139 xsdata[8]->LoadData(filename); << 159 G4ParticleDefinition* nitrogenDef = G4IonTable::GetIonTable()->GetIon(7,14);
140 << 160 G4ParticleDefinition* oxygenDef = G4IonTable::GetIonTable()->GetIon(8,16);
141 filename = "dna/sigma_ionisation_si_rudd"; << 161 G4ParticleDefinition* siliconDef = G4IonTable::GetIonTable()->GetIon(14,28);
142 xsdata[14] = new G4DNACrossSectionDataSet(ne << 162 G4ParticleDefinition* ironDef = G4IonTable::GetIonTable()->GetIon(26,56);
143 xsdata[14]->LoadData(filename); << 163 //
144 << 164
145 filename = "dna/sigma_ionisation_fe_rudd"; << 165 G4String proton;
146 xsdata[26] = new G4DNACrossSectionDataSet(ne << 166 G4String hydrogen;
147 xsdata[26]->LoadData(filename); << 167 G4String alphaPlusPlus;
148 << 168 G4String alphaPlus;
149 filename = "dna/sigma_ionisation_alphaplus_r << 169 G4String helium;
150 xsalphaplus = new G4DNACrossSectionDataSet(n << 170 G4String lithium;
151 xsalphaplus->LoadData(filename); << 171 G4String beryllium;
152 << 172 G4String boron;
153 filename = "dna/sigma_ionisation_he_rudd"; << 173 G4String carbon;
154 xshelium = new G4DNACrossSectionDataSet(new << 174 G4String nitrogen;
155 xshelium->LoadData(filename); << 175 G4String oxygen;
156 << 176 G4String silicon;
157 // to avoid possible threading problem fill << 177 G4String iron;
158 auto water = G4NistManager::Instance()->Find << 178
159 fpWaterDensity = << 179 G4double scaleFactor = 1 * m*m;
160 G4DNAMolecularMaterial::Instance()->GetNum << 180
>> 181 // LIMITS AND DATA
>> 182
>> 183 // **********************************************************************************************
>> 184
>> 185 proton = protonDef->GetParticleName();
>> 186 tableFile[proton] = fileProton;
>> 187 lowEnergyLimit[proton] = lowEnergyLimitForA[1];
>> 188 highEnergyLimit[proton] = 500. * keV;
>> 189
>> 190 // Cross section
>> 191
>> 192 G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 193 eV,
>> 194 scaleFactor );
>> 195 tableProton->LoadData(fileProton);
>> 196 tableData[proton] = tableProton;
>> 197
>> 198 // **********************************************************************************************
>> 199
>> 200 hydrogen = hydrogenDef->GetParticleName();
>> 201 tableFile[hydrogen] = fileHydrogen;
>> 202
>> 203 lowEnergyLimit[hydrogen] = lowEnergyLimitForA[1];
>> 204 highEnergyLimit[hydrogen] = 100. * MeV;
>> 205
>> 206 // Cross section
>> 207
>> 208 G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 209 eV,
>> 210 scaleFactor );
>> 211 tableHydrogen->LoadData(fileHydrogen);
>> 212
>> 213 tableData[hydrogen] = tableHydrogen;
>> 214
>> 215 // **********************************************************************************************
>> 216
>> 217 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
>> 218 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
>> 219
>> 220 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForA[4];
>> 221 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
>> 222
>> 223 // Cross section
>> 224
>> 225 G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 226 eV,
>> 227 scaleFactor );
>> 228 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
>> 229
>> 230 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
>> 231
>> 232 // **********************************************************************************************
>> 233
>> 234 alphaPlus = alphaPlusDef->GetParticleName();
>> 235 tableFile[alphaPlus] = fileAlphaPlus;
>> 236
>> 237 lowEnergyLimit[alphaPlus] = lowEnergyLimitForA[4];
>> 238 highEnergyLimit[alphaPlus] = 400. * MeV;
>> 239
>> 240 // Cross section
>> 241
>> 242 G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 243 eV,
>> 244 scaleFactor );
>> 245 tableAlphaPlus->LoadData(fileAlphaPlus);
>> 246 tableData[alphaPlus] = tableAlphaPlus;
>> 247
>> 248 // **********************************************************************************************
>> 249
>> 250 helium = heliumDef->GetParticleName();
>> 251 tableFile[helium] = fileHelium;
>> 252
>> 253 lowEnergyLimit[helium] = lowEnergyLimitForA[4];
>> 254 highEnergyLimit[helium] = 400. * MeV;
>> 255
>> 256 // Cross section
>> 257
>> 258 G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 259 eV,
>> 260 scaleFactor );
>> 261 tableHelium->LoadData(fileHelium);
>> 262 tableData[helium] = tableHelium;
>> 263
>> 264 // **********************************************************************************************
>> 265
>> 266 lithium = lithiumDef->GetParticleName();
>> 267 tableFile[lithium] = fileLithium;
>> 268
>> 269 //SI
>> 270 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
>> 271 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
>> 272 lowEnergyLimit[lithium] = 0.5*7*MeV;
>> 273 highEnergyLimit[lithium] = 1e6*7*MeV;
>> 274 //
>> 275
>> 276 // Cross section
>> 277
>> 278 G4DNACrossSectionDataSet* tableLithium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 279 eV,
>> 280 scaleFactor );
>> 281 tableLithium->LoadData(fileLithium);
>> 282 tableData[lithium] = tableLithium;
>> 283
>> 284 // **********************************************************************************************
>> 285
>> 286 beryllium = berylliumDef->GetParticleName();
>> 287 tableFile[beryllium] = fileBeryllium;
>> 288
>> 289 //SI
>> 290 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
>> 291 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
>> 292 lowEnergyLimit[beryllium] = 0.5*9*MeV;
>> 293 highEnergyLimit[beryllium] = 1e6*9*MeV;
>> 294 //
>> 295
>> 296 // Cross section
>> 297
>> 298 G4DNACrossSectionDataSet* tableBeryllium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 299 eV,
>> 300 scaleFactor );
>> 301 tableBeryllium->LoadData(fileBeryllium);
>> 302 tableData[beryllium] = tableBeryllium;
>> 303
>> 304 // **********************************************************************************************
>> 305
>> 306 boron = boronDef->GetParticleName();
>> 307 tableFile[boron] = fileBoron;
>> 308
>> 309 //SI
>> 310 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
>> 311 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
>> 312 lowEnergyLimit[boron] = 0.5*11*MeV;
>> 313 highEnergyLimit[boron] = 1e6*11*MeV;
>> 314 //
>> 315
>> 316 // Cross section
>> 317
>> 318 G4DNACrossSectionDataSet* tableBoron = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 319 eV,
>> 320 scaleFactor );
>> 321 tableBoron->LoadData(fileBoron);
>> 322 tableData[boron] = tableBoron;
>> 323
>> 324 // **********************************************************************************************
>> 325
>> 326 carbon = carbonDef->GetParticleName();
>> 327 tableFile[carbon] = fileCarbon;
>> 328
>> 329 //SI
>> 330 //lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
>> 331 //highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
>> 332 lowEnergyLimit[carbon] = 0.5*12*MeV;
>> 333 highEnergyLimit[carbon] = 1e6*12*MeV;
>> 334 //
>> 335
>> 336 // Cross section
>> 337
>> 338 G4DNACrossSectionDataSet* tableCarbon = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 339 eV,
>> 340 scaleFactor );
>> 341 tableCarbon->LoadData(fileCarbon);
>> 342 tableData[carbon] = tableCarbon;
>> 343
>> 344 // **********************************************************************************************
>> 345
>> 346 oxygen = oxygenDef->GetParticleName();
>> 347 tableFile[oxygen] = fileOxygen;
>> 348
>> 349 //SI
>> 350 //lowEnergyLimit[oxygen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
>> 351 //highEnergyLimit[oxygen] = 1e6* particle->GetAtomicMass()* MeV;
>> 352 lowEnergyLimit[oxygen] = 0.5*16*MeV;
>> 353 highEnergyLimit[oxygen] = 1e6*16*MeV;
>> 354 //
>> 355
>> 356 // Cross section
>> 357
>> 358 G4DNACrossSectionDataSet* tableOxygen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 359 eV,
>> 360 scaleFactor );
>> 361 tableOxygen->LoadData(fileOxygen);
>> 362 tableData[oxygen] = tableOxygen;
>> 363
>> 364 // **********************************************************************************************
>> 365
>> 366 nitrogen = nitrogenDef->GetParticleName();
>> 367 tableFile[nitrogen] = fileNitrogen;
>> 368
>> 369 //SI
>> 370 //lowEnergyLimit[nitrogen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
>> 371 //highEnergyLimit[nitrogen] = 1e6* particle->GetAtomicMass()* MeV;
>> 372 lowEnergyLimit[nitrogen] = 0.5*14*MeV;
>> 373 highEnergyLimit[nitrogen] = 1e6*14*MeV;
>> 374 //
>> 375
>> 376 // Cross section
>> 377
>> 378 G4DNACrossSectionDataSet* tableNitrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 379 eV,
>> 380 scaleFactor );
>> 381 tableNitrogen->LoadData(fileNitrogen);
>> 382 tableData[nitrogen] = tableNitrogen;
>> 383
>> 384 // **********************************************************************************************
>> 385
>> 386 silicon = siliconDef->GetParticleName();
>> 387 tableFile[silicon] = fileSilicon;
>> 388
>> 389 //lowEnergyLimit[silicon] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
>> 390 //highEnergyLimit[silicon] = 1e6* particle->GetAtomicMass()* MeV;
>> 391 lowEnergyLimit[silicon] = 0.5*28*MeV;
>> 392 highEnergyLimit[silicon] = 1e6*28*MeV;
>> 393 //
>> 394
>> 395 // Cross section
>> 396
>> 397 G4DNACrossSectionDataSet* tableSilicon = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 398 eV,
>> 399 scaleFactor );
>> 400 tableSilicon->LoadData(fileSilicon);
>> 401 tableData[silicon] = tableSilicon;
>> 402
>> 403 // **********************************************************************************************
>> 404
>> 405 iron = ironDef->GetParticleName();
>> 406 tableFile[iron] = fileIron;
>> 407
>> 408 //SI
>> 409 //lowEnergyLimit[iron] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
>> 410 //highEnergyLimit[iron] = 1e6* particle->GetAtomicMass()* MeV;
>> 411 lowEnergyLimit[iron] = 0.5*56*MeV;
>> 412 highEnergyLimit[iron] = 1e6*56*MeV;
>> 413 //
>> 414
>> 415 // Cross section
>> 416
>> 417 G4DNACrossSectionDataSet* tableIron = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 418 eV,
>> 419 scaleFactor );
>> 420 tableIron->LoadData(fileIron);
>> 421 tableData[iron] = tableIron;
>> 422
>> 423 // **********************************************************************************************
>> 424
>> 425 // SI: not anymore
>> 426 // ZF Following lines can be replaced by:
>> 427 // SetLowEnergyLimit(lowEnergyLimit[particle->GetParticleName()]);
>> 428 // SetHighEnergyLimit(highEnergyLimit[particle->GetParticleName()]);
>> 429 // at least for HZE
>> 430
>> 431 if (particle==protonDef)
>> 432 {
>> 433 SetLowEnergyLimit(lowEnergyLimit[proton]);
>> 434 SetHighEnergyLimit(highEnergyLimit[proton]);
>> 435 }
>> 436
>> 437 if (particle==hydrogenDef)
>> 438 {
>> 439 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
>> 440 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
>> 441 }
>> 442
>> 443 if (particle==heliumDef)
>> 444 {
>> 445 SetLowEnergyLimit(lowEnergyLimit[helium]);
>> 446 SetHighEnergyLimit(highEnergyLimit[helium]);
>> 447 }
>> 448
>> 449 if (particle==alphaPlusDef)
>> 450 {
>> 451 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
>> 452 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
>> 453 }
>> 454
>> 455 if (particle==alphaPlusPlusDef)
>> 456 {
>> 457 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
>> 458 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
>> 459 }
>> 460
>> 461 if (particle==lithiumDef)
>> 462 {
>> 463 SetLowEnergyLimit(lowEnergyLimit[lithium]);
>> 464 SetHighEnergyLimit(highEnergyLimit[lithium]);
>> 465 }
>> 466
>> 467 if (particle==berylliumDef)
>> 468 {
>> 469 SetLowEnergyLimit(lowEnergyLimit[beryllium]);
>> 470 SetHighEnergyLimit(highEnergyLimit[beryllium]);
>> 471 }
>> 472
>> 473 if (particle==boronDef)
>> 474 {
>> 475 SetLowEnergyLimit(lowEnergyLimit[boron]);
>> 476 SetHighEnergyLimit(highEnergyLimit[boron]);
>> 477 }
>> 478
>> 479 if (particle==carbonDef)
>> 480 {
>> 481 SetLowEnergyLimit(lowEnergyLimit[carbon]);
>> 482 SetHighEnergyLimit(highEnergyLimit[carbon]);
>> 483 }
>> 484
>> 485 if (particle==nitrogenDef)
>> 486 {
>> 487 SetLowEnergyLimit(lowEnergyLimit[nitrogen]);
>> 488 SetHighEnergyLimit(highEnergyLimit[nitrogen]);
>> 489 }
>> 490
>> 491 if (particle==oxygenDef)
>> 492 {
>> 493 SetLowEnergyLimit(lowEnergyLimit[oxygen]);
>> 494 SetHighEnergyLimit(highEnergyLimit[oxygen]);
>> 495 }
>> 496
>> 497 if (particle==siliconDef)
>> 498 {
>> 499 SetLowEnergyLimit(lowEnergyLimit[silicon]);
>> 500 SetHighEnergyLimit(highEnergyLimit[silicon]);
>> 501 }
>> 502
>> 503 if (particle==ironDef)
>> 504 {
>> 505 SetLowEnergyLimit(lowEnergyLimit[iron]);
>> 506 SetHighEnergyLimit(highEnergyLimit[iron]);
>> 507 }
>> 508
>> 509 //----------------------------------------------------------------------
>> 510
>> 511 if( verboseLevel>0 )
>> 512 {
>> 513 G4cout << "Rudd ionisation model is initialized " << G4endl
>> 514 << "Energy range: "
>> 515 << LowEnergyLimit() / eV << " eV - "
>> 516 << HighEnergyLimit() / keV << " keV for "
>> 517 << particle->GetParticleName()
>> 518 << G4endl;
>> 519 }
>> 520
>> 521 // Initialize water density pointer
>> 522 fpWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
>> 523
>> 524 //
>> 525
>> 526 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
>> 527
>> 528 if (isInitialised) { return; }
>> 529 fParticleChangeForGamma = GetParticleChangeForGamma();
>> 530 isInitialised = true;
161 } 531 }
162 532
163 //....oooOO0OOooo........oooOO0OOooo........oo 533 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164 534
165 void G4DNARuddIonisationExtendedModel::Initial << 535 G4double G4DNARuddIonisationExtendedModel::CrossSectionPerVolume(const G4Material* material,
166 << 536 const G4ParticleDefinition* particleDefinition,
>> 537 G4double k,
>> 538 G4double,
>> 539 G4double)
167 { 540 {
168 if (p != fParticle) { SetParticle(p); } << 541 //SI: particleDefinition->GetParticleName() is for eg. Fe56
>> 542 // particleDefinition->GetPDGMass() is correct
>> 543 // particleDefinition->GetAtomicNumber() is correct
>> 544
>> 545 if (verboseLevel > 3)
>> 546 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" << G4endl;
>> 547
>> 548 // Calculate total cross section for model
>> 549
>> 550 G4DNAGenericIonsManager *instance;
>> 551 instance = G4DNAGenericIonsManager::Instance();
>> 552
>> 553 if (
>> 554 particleDefinition != G4Proton::ProtonDefinition()
>> 555 &&
>> 556 particleDefinition != instance->GetIon("hydrogen")
>> 557 &&
>> 558 particleDefinition != instance->GetIon("alpha++")
>> 559 &&
>> 560 particleDefinition != instance->GetIon("alpha+")
>> 561 &&
>> 562 particleDefinition != instance->GetIon("helium")
>> 563 &&
>> 564 // SI
>> 565 //particleDefinition != instance->GetIon("carbon")
>> 566 //&&
>> 567 //particleDefinition != instance->GetIon("nitrogen")
>> 568 //&&
>> 569 //particleDefinition != instance->GetIon("oxygen")
>> 570 //&&
>> 571 //particleDefinition != instance->GetIon("iron")
>> 572 particleDefinition != G4IonTable::GetIonTable()->GetIon(3,7)
>> 573 &&
>> 574 particleDefinition != G4IonTable::GetIonTable()->GetIon(4,9)
>> 575 &&
>> 576 particleDefinition != G4IonTable::GetIonTable()->GetIon(5,11)
>> 577 &&
>> 578 particleDefinition != G4IonTable::GetIonTable()->GetIon(6,12)
>> 579 &&
>> 580 particleDefinition != G4IonTable::GetIonTable()->GetIon(7,14)
>> 581 &&
>> 582 particleDefinition != G4IonTable::GetIonTable()->GetIon(8,16)
>> 583 &&
>> 584 particleDefinition != G4IonTable::GetIonTable()->GetIon(14,28)
>> 585 &&
>> 586 particleDefinition != G4IonTable::GetIonTable()->GetIon(26,56)
>> 587 //
>> 588 )
>> 589
>> 590 return 0;
>> 591
>> 592 G4double lowLim = 0;
>> 593
>> 594 if ( particleDefinition == G4Proton::ProtonDefinition()
>> 595 || particleDefinition == instance->GetIon("hydrogen")
>> 596 )
>> 597
>> 598 lowLim = lowEnergyLimitOfModelForA[1];
>> 599
>> 600 else if ( particleDefinition == instance->GetIon("alpha++")
>> 601 || particleDefinition == instance->GetIon("alpha+")
>> 602 || particleDefinition == instance->GetIon("helium")
>> 603 )
>> 604
>> 605 lowLim = lowEnergyLimitOfModelForA[4];
>> 606
>> 607 else lowLim = lowEnergyLimitOfModelForA[5];
>> 608
>> 609 G4double highLim = 0;
>> 610 G4double sigma=0;
>> 611
>> 612
>> 613 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
>> 614
>> 615 const G4String& particleName = particleDefinition->GetParticleName();
>> 616
>> 617 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
>> 618 pos2 = highEnergyLimit.find(particleName);
>> 619
>> 620 if (pos2 != highEnergyLimit.end())
>> 621 {
>> 622 highLim = pos2->second;
>> 623 }
169 624
170 // particle change object may be externally << 625 if (k <= highLim)
171 if (nullptr == fParticleChangeForGamma) { << 626 {
172 fParticleChangeForGamma = GetParticleChang << 627
173 } << 628 //SI : XS must not be zero otherwise sampling of secondaries method ignored
>> 629
>> 630 if (k < lowLim) k = lowLim;
>> 631
>> 632 //
>> 633
>> 634 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 635 pos = tableData.find(particleName);
>> 636
>> 637 if (pos != tableData.end())
>> 638 {
>> 639 G4DNACrossSectionDataSet* table = pos->second;
>> 640 if (table != 0)
>> 641 {
>> 642 sigma = table->FindValue(k);
>> 643 }
>> 644 }
>> 645 else
>> 646 {
>> 647 G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume","em0002",
>> 648 FatalException,"Model not applicable to particle type.");
>> 649 }
>> 650
>> 651 } // if (k >= lowLim && k < highLim)
>> 652
>> 653 if (verboseLevel > 2)
>> 654 {
>> 655 G4cout << "__________________________________" << G4endl;
>> 656 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl;
>> 657 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
>> 658 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
>> 659 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
>> 660 //G4cout << " - Cross section per water molecule (cm^-1)="
>> 661 //<< sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
>> 662 G4cout << "G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl;
174 663
175 // initialisation once in each thread <<
176 if (!isInitialised) { <<
177 isInitialised = true; <<
178 const G4String& pname = fParticle->GetPart <<
179 if (pname == "proton") { <<
180 idx = 1; <<
181 xscurrent = xsdata[1]; <<
182 fElow = fLowestEnergy; <<
183 } else if (pname == "hydrogen") { <<
184 idx = 0; <<
185 xscurrent = xsdata[0]; <<
186 fElow = fLowestEnergy; <<
187 } else if (pname == "alpha") { <<
188 idx = 1; <<
189 xscurrent = xsdata[2]; <<
190 isHelium = true; <<
191 fElow = fLimitEnergy; <<
192 } else if (pname == "alpha+") { <<
193 idx = 1; <<
194 isHelium = true; <<
195 xscurrent = xsalphaplus; <<
196 fElow = fLimitEnergy; <<
197 // The following values are provided by <<
198 slaterEffectiveCharge[0]=2.0; <<
199 slaterEffectiveCharge[1]=2.0; <<
200 slaterEffectiveCharge[2]=2.0; <<
201 sCoefficient[0]=0.7; <<
202 sCoefficient[1]=0.15; <<
203 sCoefficient[2]=0.15; <<
204 } else if (pname == "helium") { <<
205 idx = 0; <<
206 isHelium = true; <<
207 fElow = fLimitEnergy; <<
208 xscurrent = xshelium; <<
209 slaterEffectiveCharge[0]=1.7; <<
210 slaterEffectiveCharge[1]=1.15; <<
211 slaterEffectiveCharge[2]=1.15; <<
212 sCoefficient[0]=0.5; <<
213 sCoefficient[1]=0.25; <<
214 sCoefficient[2]=0.25; <<
215 } else { <<
216 isIon = true; <<
217 idx = -1; <<
218 xscurrent = xsdata[1]; <<
219 fElow = fLowestEnergy; <<
220 } <<
221 // defined stationary mode <<
222 statCode = G4EmParameters::Instance()->DNA <<
223 <<
224 // initialise atomic de-excitation <<
225 fAtomDeexcitation = G4LossTableManager::In <<
226 <<
227 if (verbose > 0) { <<
228 G4cout << "### G4DNARuddIonisationExtend <<
229 << "/n idx=" << idx << " isIon=" << <<
230 << " isHelium=" << isHelium << G4endl; <<
231 } 664 }
232 } << 665
>> 666 return sigma*waterDensity;
>> 667
233 } 668 }
234 669
235 //....oooOO0OOooo........oooOO0OOooo........oo 670 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236 671
237 void G4DNARuddIonisationExtendedModel::SetPart << 672 void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
>> 673 const G4MaterialCutsCouple* couple,
>> 674 const G4DynamicParticle* particle,
>> 675 G4double,
>> 676 G4double)
238 { 677 {
239 fParticle = p; << 678 //SI: particle->GetDefinition()->GetParticleName() is for eg. Fe56
240 fMass = p->GetPDGMass(); << 679 // particle->GetDefinition()->GetPDGMass() is correct
241 fMassRate = (isIon) ? CLHEP::proton_mass_c2/ << 680 // particle->GetDefinition()->GetAtomicNumber() is correct
>> 681 // particle->GetDefinition()->GetAtomicMass() is correct
>> 682
>> 683 if (verboseLevel > 3)
>> 684 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl;
>> 685
>> 686 G4double lowLim = 0;
>> 687 G4double highLim = 0;
>> 688
>> 689 // ZF: the following line summarizes the commented part
>> 690
>> 691 if(particle->GetDefinition()->GetAtomicMass() <= 4) lowLim = killBelowEnergyForA[particle->GetDefinition()->GetAtomicMass()];
>> 692
>> 693 else lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
>> 694
>> 695 /*
>> 696
>> 697 if(particle->GetDefinition()->GetAtomicMass() >= 5) lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
>> 698
>> 699 if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
>> 700 || particle->GetDefinition() == instance->GetIon("hydrogen")
>> 701 )
>> 702
>> 703 lowLim = killBelowEnergyForA[1];
>> 704
>> 705 if ( particle->GetDefinition() == instance->GetIon("alpha++")
>> 706 || particle->GetDefinition() == instance->GetIon("alpha+")
>> 707 || particle->GetDefinition() == instance->GetIon("helium")
>> 708 )
>> 709
>> 710 lowLim = killBelowEnergyForA[4];
>> 711
>> 712 */
>> 713 //
>> 714
>> 715 G4double k = particle->GetKineticEnergy();
>> 716
>> 717 const G4String& particleName = particle->GetDefinition()->GetParticleName();
>> 718
>> 719 // SI - the following is useless since lowLim is already defined
>> 720 /*
>> 721 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
>> 722 pos1 = lowEnergyLimit.find(particleName);
>> 723
>> 724 if (pos1 != lowEnergyLimit.end())
>> 725 {
>> 726 lowLim = pos1->second;
>> 727 }
>> 728 */
>> 729
>> 730 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
>> 731 pos2 = highEnergyLimit.find(particleName);
>> 732
>> 733 if (pos2 != highEnergyLimit.end()) highLim = pos2->second;
>> 734
>> 735 if (k >= lowLim && k <= highLim)
>> 736
>> 737 // SI: no strict limits, like in the non extended version of the model
>> 738 {
>> 739 G4ParticleDefinition* definition = particle->GetDefinition();
>> 740 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
>> 741 /*
>> 742 G4double particleMass = definition->GetPDGMass();
>> 743 G4double totalEnergy = k + particleMass;
>> 744 G4double pSquare = k*(totalEnergy+particleMass);
>> 745 G4double totalMomentum = std::sqrt(pSquare);
>> 746 */
>> 747
>> 748 G4int ionizationShell = RandomSelect(k,particleName);
>> 749
>> 750 // sample deexcitation
>> 751 // here we assume that H_{2}O electronic levels are the same as Oxygen.
>> 752 // this can be considered true with a rough 10% error in energy on K-shell,
>> 753
>> 754 G4double bindingEnergy = 0;
>> 755 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
>> 756
>> 757 //SI: additional protection if tcs interpolation method is modified
>> 758 if (k<bindingEnergy) return;
>> 759 //
>> 760
>> 761 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
>> 762
>> 763 G4int Z = 8;
>> 764
>> 765 G4ThreeVector deltaDirection =
>> 766 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
>> 767 Z, ionizationShell,
>> 768 couple->GetMaterial());
>> 769
>> 770 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
>> 771 fvect->push_back(dp);
>> 772
>> 773 fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
>> 774
>> 775 // SI: the following lines are not needed anymore
>> 776 /*
>> 777 G4double cosTheta = 0.;
>> 778 G4double phi = 0.;
>> 779 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi, ionizationShell);
>> 780
>> 781 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
>> 782 G4double dirX = sinTheta*std::cos(phi);
>> 783 G4double dirY = sinTheta*std::sin(phi);
>> 784 G4double dirZ = cosTheta;
>> 785 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
>> 786 deltaDirection.rotateUz(primaryDirection);
>> 787 */
>> 788
>> 789 // Ignored for ions on electrons
>> 790 /*
>> 791 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
>> 792
>> 793 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
>> 794 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
>> 795 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
>> 796 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
>> 797 finalPx /= finalMomentum;
>> 798 finalPy /= finalMomentum;
>> 799 finalPz /= finalMomentum;
>> 800
>> 801 G4ThreeVector direction;
>> 802 direction.set(finalPx,finalPy,finalPz);
>> 803
>> 804 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
>> 805 */
>> 806
>> 807 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
>> 808 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
>> 809
>> 810 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
>> 811
>> 812 // SI: only atomic deexcitation from K shell is considered
>> 813 if(fAtomDeexcitation && ionizationShell == 4)
>> 814 {
>> 815 const G4AtomicShell* shell
>> 816 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
>> 817 secNumberInit = fvect->size();
>> 818 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
>> 819 secNumberFinal = fvect->size();
>> 820
>> 821 if(secNumberFinal > secNumberInit)
>> 822 {
>> 823 for (size_t i=secNumberInit; i<secNumberFinal; ++i)
>> 824 {
>> 825 //Check if there is enough residual energy
>> 826 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
>> 827 {
>> 828 //Ok, this is a valid secondary: keep it
>> 829 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
>> 830 }
>> 831 else
>> 832 {
>> 833 //Invalid secondary: not enough energy to create it!
>> 834 //Keep its energy in the local deposit
>> 835 delete (*fvect)[i];
>> 836 (*fvect)[i]=0;
>> 837 }
>> 838 }
>> 839 }
>> 840
>> 841 }
>> 842
>> 843 //This should never happen
>> 844 if(bindingEnergy < 0.0)
>> 845 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
>> 846 "em2050",FatalException,"Negative local energy deposit");
>> 847
>> 848 //bindingEnergy has been decreased
>> 849 //by the amount of energy taken away by deexc. products
>> 850 if (!statCode)
>> 851 {
>> 852 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
>> 853 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
>> 854 }
>> 855 else
>> 856 {
>> 857 fParticleChangeForGamma->SetProposedKineticEnergy(k);
>> 858 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
>> 859 }
>> 860
>> 861 // TEST //////////////////////////
>> 862 // if (secondaryKinetic<0) abort();
>> 863 // if (scatteredEnergy<0) abort();
>> 864 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
>> 865 // if (k-scatteredEnergy<0) abort();
>> 866 /////////////////////////////////
>> 867
>> 868 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
>> 869 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
>> 870 ionizationShell,
>> 871 theIncomingTrack);
>> 872 }
>> 873
>> 874 // SI - not useful since low energy of model is 0 eV
>> 875
>> 876 if (k < lowLim)
>> 877 {
>> 878 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
>> 879 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
>> 880 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k);
>> 881 }
>> 882
242 } 883 }
243 884
244 //....oooOO0OOooo........oooOO0OOooo........oo << 885 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 886
>> 887 G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
>> 888 G4double k,
>> 889 G4int shell)
>> 890 {
>> 891 //-- Fast sampling method -----
>> 892 G4double proposed_energy;
>> 893 G4double random1;
>> 894 G4double value_sampling;
>> 895 G4double max1;
>> 896
>> 897 do
>> 898 {
>> 899 proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling
>> 900
>> 901 max1=0.;
>> 902
>> 903 for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1)
>> 904 max1=RejectionFunction(particleDefinition, k, en, shell);
245 905
246 G4double << 906 random1 = G4UniformRand()*max1;
247 G4DNARuddIonisationExtendedModel::CrossSection << 907
248 << 908 value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
249 << 909
250 << 910 } while(random1 > value_sampling);
251 { << 911
252 // check if model is applicable for given ma << 912 return(proposed_energy);
253 G4double density = (material->GetIndex() < f <<
254 ? (*fpWaterDensity)[material->GetIndex()] <<
255 if (0.0 == density) { return 0.0; } <<
256 <<
257 // ion may be different <<
258 if (fParticle != part) { SetParticle(part); <<
259 <<
260 // ion shoud be stopped - check on kinetic e <<
261 if (kinE < fLowestEnergy) { return DBL_MAX; <<
262 <<
263 G4double e = kinE*fMassRate; <<
264 <<
265 G4double sigma = (e > fElow) ? xscurrent->Fi <<
266 : xscurrent->FindValue(fElow) * e / fElow; <<
267 <<
268 if (idx == -1) { <<
269 sigma *= fEmCorrections->EffectiveChargeSq <<
270 } <<
271 <<
272 sigma *= density; <<
273 <<
274 if (verbose > 1) { <<
275 G4cout << "G4DNARuddIonisationExtendedMode <<
276 << " Ekin(keV)=" << kinE/CLHEP::keV <<
277 << " sigma(cm^2)=" << sigma/CLHEP:: <<
278 } <<
279 return sigma; <<
280 } 913 }
281 914
282 //....oooOO0OOooo........oooOO0OOooo........oo << 915 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
283 916
284 void << 917 // The following section is not used anymore but is kept for memory
285 G4DNARuddIonisationExtendedModel::SampleSecond << 918 // GetAngularDistribution()->SampleDirectionForShell is used instead
286 <<
287 <<
288 <<
289 { <<
290 const G4ParticleDefinition* pd = dpart->GetD <<
291 if (fParticle != pd) { SetParticle(pd); } <<
292 <<
293 // stop ion with energy below low energy lim <<
294 G4double kinE = dpart->GetKineticEnergy(); <<
295 // ion shoud be stopped - check on kinetic e <<
296 if (kinE <= fLowestEnergy) { <<
297 fParticleChangeForGamma->SetProposedKineti <<
298 fParticleChangeForGamma->ProposeTrackStatu <<
299 fParticleChangeForGamma->ProposeLocalEnerg <<
300 return; <<
301 } <<
302 <<
303 G4int shell = SelectShell(kinE*fMassRate); <<
304 G4double bindingEnergy = (useDNAWaterStructu <<
305 ? waterStructure.IonisationEnergy(shell) : <<
306 919
307 //Si: additional protection if tcs interpola << 920 /*
308 if (kinE < bindingEnergy) { return; } << 921 void G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
309 << 922 G4double k,
310 G4double esec = SampleElectronEnergy(kinE, s << 923 G4double secKinetic,
311 G4double esum = 0.0; << 924 G4double & cosTheta,
>> 925 G4double & phi,
>> 926 G4int shell )
>> 927 {
>> 928 G4double maxSecKinetic = 0.;
>> 929 G4double maximumEnergyTransfer = 0.;
312 930
313 // sample deexcitation << 931 // ZF. generalized & relativistic version
314 // here we assume that H2O electronic levels << 932
315 // this can be considered true with a rough << 933 if( (k/MeV)/(particleDefinition->GetPDGMass()/MeV) <= 0.1 )
316 G4int Z = 8; << 934 {
317 G4ThreeVector deltaDir = << 935 maximumEnergyTransfer= 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
318 GetAngularDistribution()->SampleDirectionF << 936 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
319 << 937 }
320 // SI: only atomic deexcitation from K shell << 938 else
321 if(fAtomDeexcitation != nullptr && shell == << 939 {
322 auto as = G4AtomicShellEnumerator(0); << 940 G4double approx_nuc_number = particleDefinition->GetPDGMass() / proton_mass_c2;
323 auto ashell = fAtomDeexcitation->GetAtomic << 941 G4double en_per_nucleon = k/approx_nuc_number;
324 fAtomDeexcitation->GenerateParticles(fvect << 942 G4double beta2 = 1. - 1./pow( (1.+(en_per_nucleon/electron_mass_c2)*(electron_mass_c2/proton_mass_c2)), 2.);
325 << 943 G4double gamma = 1./sqrt(1.-beta2);
326 // compute energy sum from de-excitation << 944 maximumEnergyTransfer = 2.*electron_mass_c2*(gamma*gamma-1.)/(1.+2.*gamma*(electron_mass_c2/particleDefinition->GetPDGMass())+pow(electron_mass_c2/particleDefinition->GetPDGMass(), 2.) );
327 for (auto const & ptr : *fvect) { << 945 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
328 esum += ptr->GetKineticEnergy(); << 946 }
329 } << 947
330 } << 948 maxSecKinetic = maximumEnergyTransfer-waterStructure.IonisationEnergy(shell);
331 // check energy balance << 949
332 // remaining excitation energy of water mole << 950 phi = twopi * G4UniformRand();
333 G4double exc = bindingEnergy - esum; << 951
334 << 952 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
335 // remaining projectile energy << 953 else cosTheta = (2.*G4UniformRand())-1.;
336 G4double scatteredEnergy = kinE - bindingEne << 954
337 if(scatteredEnergy < -tolerance || exc < -to << 955 }
338 G4cout << "G4DNARuddIonisationExtendedMode << 956 */
339 << "negative final E(keV)=" << scat << 957
340 << kinE/CLHEP::keV << " " << pd->G << 958 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
341 << " Edelta(keV)=" << esec/CLHEP::k << 959
342 << G4endl; << 960 G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition,
343 } << 961 G4double k,
344 << 962 G4double proposed_ws,
345 // projectile << 963 G4int ionizationLevelIndex)
346 if (!statCode) { <<
347 fParticleChangeForGamma->SetProposedKineti <<
348 fParticleChangeForGamma->ProposeLocalEnerg <<
349 } else { <<
350 fParticleChangeForGamma->SetProposedKineti <<
351 fParticleChangeForGamma->ProposeLocalEnerg <<
352 } <<
353 <<
354 // delta-electron <<
355 auto dp = new G4DynamicParticle(G4Electron: <<
356 fvect->push_back(dp); <<
357 <<
358 // create radical <<
359 const G4Track* theIncomingTrack = fParticleC <<
360 G4DNAChemistryManager::Instance()->CreateWat <<
361 theIncomingTrack); <<
362 } <<
363 <<
364 //....oooOO0OOooo........oooOO0OOooo........oo <<
365 <<
366 G4int G4DNARuddIonisationExtendedModel::Select <<
367 { <<
368 G4double sum = 0.0; <<
369 G4double xs; <<
370 for (G4int i=0; i<5; ++i) { <<
371 auto ptr = xscurrent->GetComponent(i); <<
372 xs = (e > fElow) ? ptr->FindValue(e) : ptr <<
373 sum += xs; <<
374 fTemp[i] = sum; <<
375 } <<
376 sum *= G4UniformRand(); <<
377 for (G4int i=0; i<5; ++i) { <<
378 if (sum <= fTemp[i]) { return i; } <<
379 } <<
380 return 0; <<
381 } <<
382 <<
383 //....oooOO0OOooo........oooOO0OOooo........oo <<
384 <<
385 G4double G4DNARuddIonisationExtendedModel::Max <<
386 { <<
387 // kinematic limit <<
388 G4double tau = kine/fMass; <<
389 G4double gam = 1.0 + tau; <<
390 G4double emax = 2.0*CLHEP::electron_mass_c2* <<
391 <<
392 // Initialisation of sampling <<
393 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2; <<
394 if (shell == 4) { <<
395 //Data For Liquid Water K SHELL from Dingf <<
396 A1 = 1.25; <<
397 B1 = 0.5; <<
398 C1 = 1.00; <<
399 D1 = 1.00; <<
400 E1 = 3.00; <<
401 A2 = 1.10; <<
402 B2 = 1.30; <<
403 C2 = 1.00; <<
404 D2 = 0.00; <<
405 alphaConst = 0.66; <<
406 } else { <<
407 //Data For Liquid Water from Dingfelder (P <<
408 A1 = 1.02; <<
409 B1 = 82.0; <<
410 C1 = 0.45; <<
411 D1 = -0.80; <<
412 E1 = 0.38; <<
413 A2 = 1.07; <<
414 // Value provided by M. Dingfelder (priv. <<
415 B2 = 11.6; <<
416 C2 = 0.60; <<
417 D2 = 0.04; <<
418 alphaConst = 0.64; <<
419 } <<
420 bEnergy = Bj[shell]; <<
421 G4double v2 = 0.25*emax/(bEnergy*gam*gam); <<
422 v = std::sqrt(v2); <<
423 u = Ry/bEnergy; <<
424 wc = 4.*v2 - 2.*v - 0.25*u; <<
425 <<
426 G4double L1 = (C1 * fGpow->powA(v, D1)) / (1 <<
427 G4double L2 = C2 * fGpow->powA(v, D2); <<
428 G4double H1 = (A1 * G4Log(1. + v2)) / (v2 + <<
429 G4double H2 = (A2 / v2) + (B2 / (v2 * v2)); <<
430 <<
431 F1 = L1 + H1; <<
432 F2 = (L2 * H2) / (L2 + H2); <<
433 return emax; <<
434 } <<
435 <<
436 //....oooOO0OOooo........oooOO0OOooo........oo <<
437 <<
438 G4double G4DNARuddIonisationExtendedModel::Sam <<
439 <<
440 { <<
441 G4double emax = MaxEnergy(kine, shell); <<
442 // compute cumulative probability function <<
443 G4double step = 1*CLHEP::eV; <<
444 auto nn = (G4int)(emax/step); <<
445 nn = std::min(std::max(nn, 10), 100); <<
446 step = emax/(G4double)nn; <<
447 <<
448 // find max probability <<
449 G4double pmax = ProbabilityFunction(kine, 0. <<
450 //G4cout << "## E(keV)=" << kine/keV << " em <<
451 // << " pmax(0)=" << pmax << " shell=" <<
452 <<
453 G4double e0 = 0.0; // energy with max probab <<
454 // 2 areas after point with max probability <<
455 G4double e1 = emax; <<
456 G4double e2 = emax; <<
457 G4double p1 = 0.0; <<
458 G4double p2 = 0.0; <<
459 const G4double f = 0.25; <<
460 <<
461 // find max probability <<
462 G4double e = 0.0; <<
463 G4double p = 0.0; <<
464 for (G4int i=0; i<nn; ++i) { <<
465 e += step; <<
466 p = ProbabilityFunction(kine, e, shell); <<
467 if (p > pmax) { <<
468 pmax = p; <<
469 e0 = e; <<
470 } else { <<
471 break; <<
472 } <<
473 } <<
474 // increase step to be more effective <<
475 step *= 2.0; <<
476 // 2-nd area <<
477 for (G4int i=0; i<nn; ++i) { <<
478 e += step; <<
479 if (std::abs(e - emax) < step) { <<
480 e1 = emax; <<
481 break; <<
482 } <<
483 p = ProbabilityFunction(kine, e, shell); <<
484 if (p < f*pmax) { <<
485 p1 = p; <<
486 e1 = e; <<
487 break; <<
488 } <<
489 } <<
490 // 3-d area <<
491 if (e < emax) { <<
492 for (G4int i=0; i<nn; ++i) { <<
493 e += step; <<
494 if (std::abs(e - emax) < step) { <<
495 e2 = emax; <<
496 break; <<
497 } <<
498 p = ProbabilityFunction(kine, e, shell); <<
499 if (p < f*p1) { <<
500 p2 = p; <<
501 e2 = e; <<
502 break; <<
503 } <<
504 } <<
505 } <<
506 pmax *= 1.05; <<
507 // regression method with 3 regions <<
508 G4double s0 = pmax*e1; <<
509 G4double s1 = s0 + p1 * (e2 - e1); <<
510 G4double s2 = s1 + p2 * (emax - e2); <<
511 s0 = (s0 == s1) ? 1.0 : s0 / s2; <<
512 s1 = (s1 == s2) ? 1.0 : s1 / s2; <<
513 <<
514 //G4cout << "pmax=" << pmax << " e1(keV)=" < <<
515 // << " p2=" << p2 << " s0=" << s0 << " s1 <<
516 <<
517 // sampling <<
518 G4int count = 0; <<
519 G4double ymax, y, deltae; <<
520 for (G4int i = 0; i<100000; ++i) { <<
521 G4double q = G4UniformRand(); <<
522 if (q <= s0) { <<
523 ymax = pmax; <<
524 deltae = e1 * q / s0; <<
525 } else if (q <= s1) { <<
526 ymax = p1; <<
527 deltae = e1 + (e2 - e1) * (q - s0) / (s1 <<
528 } else { <<
529 ymax = p2; <<
530 deltae = e2 + (emax - e2) * (q - s1) / ( <<
531 } <<
532 y = ProbabilityFunction(kine, deltae, shel <<
533 //G4cout << " " << i << ". deltae=" << <<
534 // << " y=" << y << " ymax=" << ymax <<
535 if (y > ymax && count < 10) { <<
536 ++count; <<
537 G4cout << "G4DNARuddIonisationExtendedMo <<
538 << fParticle->GetParticleName() << " E( <<
539 << " Edelta(keV)=" << deltae/CLHEP::keV <<
540 << " y=" << y << " ymax=" << ymax << " <<
541 } <<
542 if (ymax * G4UniformRand() <= y) { <<
543 return deltae; <<
544 } <<
545 } <<
546 deltae = std::min(e0 + step, 0.5*emax); <<
547 return deltae; <<
548 } <<
549 <<
550 //....oooOO0OOooo........oooOO0OOooo........oo <<
551 <<
552 G4double G4DNARuddIonisationExtendedModel::Pro <<
553 <<
554 <<
555 { <<
556 // Shells ids are 0 1 2 3 4 (4 is k shell) <<
557 // !!Attention, "energyTransfer" here is the <<
558 // that the secondary kinetic en <<
559 // <<
560 // ds S F1(nu) + <<
561 // ---- = G(k) * ---- ----------------- <<
562 // dw Bj (1+w)^3 * [1 + e <<
563 // <<
564 // w is the secondary electron kinetic Energ <<
565 // <<
566 // All the other parameters can be found in <<
567 // <<
568 // M.Eugene Rudd, 1988, User-Friendly model <<
569 // electrons from protons or electron collis <<
570 // <<
571 G4double w = deltae/bEnergy; <<
572 G4double x = alphaConst*(w - wc)/v; <<
573 G4double y = (x > -15.) ? 1.0 + G4Exp(x) : 1 <<
574 <<
575 G4double res = CorrectionFactor(kine, shell) <<
576 (fGpow->powN((1. + w)/u, 3) * y); <<
577 <<
578 if (isHelium) { <<
579 G4double energyTransfer = deltae + bEnergy <<
580 G4double Zeff = 2.0 - <<
581 (sCoefficient[0] * S_1s(kine, energyTran <<
582 sCoefficient[1] * S_2s(kine, energyTran <<
583 sCoefficient[2] * S_2p(kine, energyTran <<
584 <<
585 res *= Zeff * Zeff; <<
586 } <<
587 return std::max(res, 0.0); <<
588 } <<
589 <<
590 //....oooOO0OOooo........oooOO0OOooo........oo <<
591 <<
592 G4double G4DNARuddIonisationExtendedModel::Com <<
593 const G4ParticleDefinition* p, G4doub <<
594 { <<
595 if (fParticle != p) { SetParticle(p); } <<
596 MaxEnergy(e, shell); <<
597 return ProbabilityFunction(e, deltae, shell) <<
598 } <<
599 <<
600 //....oooOO0OOooo........oooOO0OOooo........oo <<
601 <<
602 G4double G4DNARuddIonisationExtendedModel::S_1 <<
603 <<
604 <<
605 <<
606 { 964 {
607 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2) << 965 const G4int j=ionizationLevelIndex;
608 // Dingfelder, in Chattanooga 2005 proceedin << 966 G4double Bj_energy, alphaConst;
>> 967 G4double Ry = 13.6*eV;
>> 968 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
>> 969
>> 970 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
>> 971
>> 972 // Following values provided by M. Dingfelder (priv. comm)
>> 973 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
>> 974
>> 975 if (j == 4)
>> 976 {
>> 977 alphaConst = 0.66;
>> 978 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
>> 979 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
>> 980 //---
>> 981 }
>> 982 else
>> 983 {
>> 984 alphaConst = 0.64;
>> 985 Bj_energy = Bj[ionizationLevelIndex];
>> 986 }
>> 987
>> 988 G4double energyTransfer = proposed_ws + Bj_energy;
>> 989 proposed_ws/=Bj_energy;
>> 990 G4DNAGenericIonsManager *instance;
>> 991 instance = G4DNAGenericIonsManager::Instance();
>> 992 G4double tau = 0.;
>> 993 G4double A_ion = 0.;
>> 994 tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
>> 995 A_ion = particleDefinition->GetAtomicMass();
>> 996
>> 997 G4double v2;
>> 998 G4double beta2;
>> 999
>> 1000 if((tau/MeV)<5.447761194e-2)
>> 1001 {
>> 1002 v2 = tau / Bj_energy;
>> 1003 beta2 = 2.*tau / electron_mass_c2;
>> 1004 }
>> 1005 // Relativistic
>> 1006 else
>> 1007 {
>> 1008 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
>> 1009 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
>> 1010 }
>> 1011
>> 1012 G4double v = std::sqrt(v2);
>> 1013 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
>> 1014 G4double rejection_term = 1.+G4Exp(alphaConst*(proposed_ws - wc) / v);
>> 1015 rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
>> 1016 //* (S/Bj_energy) ; Not needed anymore
>> 1017
>> 1018 G4bool isHelium = false;
>> 1019
>> 1020 if ( particleDefinition == G4Proton::ProtonDefinition()
>> 1021 || particleDefinition == instance->GetIon("hydrogen")
>> 1022 )
>> 1023 {
>> 1024 return(rejection_term);
>> 1025 }
>> 1026
>> 1027 else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium
>> 1028 {
>> 1029 G4double Z = particleDefinition->GetAtomicNumber();
>> 1030
>> 1031 G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.));
>> 1032 G4double Zeffion = Z*(1.-G4Exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
>> 1033 rejection_term*=Zeffion*Zeffion;
>> 1034 }
>> 1035
>> 1036 else if (particleDefinition == instance->GetIon("alpha++") )
>> 1037 {
>> 1038 isHelium = true;
>> 1039 slaterEffectiveCharge[0]=0.;
>> 1040 slaterEffectiveCharge[1]=0.;
>> 1041 slaterEffectiveCharge[2]=0.;
>> 1042 sCoefficient[0]=0.;
>> 1043 sCoefficient[1]=0.;
>> 1044 sCoefficient[2]=0.;
>> 1045 }
>> 1046
>> 1047 else if (particleDefinition == instance->GetIon("alpha+") )
>> 1048 {
>> 1049 isHelium = true;
>> 1050 slaterEffectiveCharge[0]=2.0;
>> 1051 // The following values are provided by M. Dingfelder (priv. comm)
>> 1052 slaterEffectiveCharge[1]=2.0;
>> 1053 slaterEffectiveCharge[2]=2.0;
>> 1054 //
>> 1055 sCoefficient[0]=0.7;
>> 1056 sCoefficient[1]=0.15;
>> 1057 sCoefficient[2]=0.15;
>> 1058 }
>> 1059
>> 1060 else if (particleDefinition == instance->GetIon("helium") )
>> 1061 {
>> 1062 isHelium = true;
>> 1063 slaterEffectiveCharge[0]=1.7;
>> 1064 slaterEffectiveCharge[1]=1.15;
>> 1065 slaterEffectiveCharge[2]=1.15;
>> 1066 sCoefficient[0]=0.5;
>> 1067 sCoefficient[1]=0.25;
>> 1068 sCoefficient[2]=0.25;
>> 1069 }
>> 1070
>> 1071 // if ( particleDefinition == instance->GetIon("helium")
>> 1072 // || particleDefinition == instance->GetIon("alpha+")
>> 1073 // || particleDefinition == instance->GetIon("alpha++")
>> 1074 // )
>> 1075
>> 1076 if (isHelium)
>> 1077 {
>> 1078
>> 1079 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
>> 1080
>> 1081 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
>> 1082 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
>> 1083 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
>> 1084
>> 1085 rejection_term*= zEff * zEff;
>> 1086 }
>> 1087
>> 1088 return (rejection_term);
>> 1089 }
>> 1090
>> 1091
>> 1092 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1093
>> 1094
>> 1095 G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle,
>> 1096 G4double k,
>> 1097 G4int ionizationLevelIndex)
>> 1098 {
>> 1099
>> 1100 const G4int j=ionizationLevelIndex;
>> 1101
>> 1102 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2;
>> 1103 //G4double alphaConst ;
>> 1104 G4double Bj_energy;
>> 1105
>> 1106 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
>> 1107 // Following values provided by M. Dingfelder (priv. comm)
>> 1108
>> 1109 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
>> 1110
>> 1111 if (j == 4)
>> 1112 {
>> 1113 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
>> 1114 A1 = 1.25;
>> 1115 B1 = 0.5;
>> 1116 C1 = 1.00;
>> 1117 D1 = 1.00;
>> 1118 E1 = 3.00;
>> 1119 A2 = 1.10;
>> 1120 B2 = 1.30;
>> 1121 C2 = 1.00;
>> 1122 D2 = 0.00;
>> 1123 //alphaConst = 0.66;
>> 1124 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
>> 1125 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
>> 1126 //---
>> 1127 }
>> 1128 else
>> 1129 {
>> 1130 //Data For Liquid Water from Dingfelder (Protons in Water)
>> 1131 A1 = 1.02;
>> 1132 B1 = 82.0;
>> 1133 C1 = 0.45;
>> 1134 D1 = -0.80;
>> 1135 E1 = 0.38;
>> 1136 A2 = 1.07;
>> 1137 //B2 = 14.6; From Ding Paper
>> 1138 // Value provided by M. Dingfelder (priv. comm)
>> 1139 B2 = 11.6;
>> 1140 //
>> 1141 C2 = 0.60;
>> 1142 D2 = 0.04;
>> 1143 //alphaConst = 0.64;
>> 1144
>> 1145 Bj_energy = Bj[ionizationLevelIndex];
>> 1146 }
>> 1147
>> 1148 G4double tau = 0.;
>> 1149 G4double A_ion = 0.;
>> 1150 tau = (electron_mass_c2 / particle->GetPDGMass()) * k;
>> 1151
>> 1152 A_ion = particle->GetAtomicMass();
>> 1153
>> 1154 G4double v2;
>> 1155 G4double beta2;
>> 1156 if((tau/MeV)<5.447761194e-2)
>> 1157 {
>> 1158 v2 = tau / Bj_energy;
>> 1159 beta2 = 2.*tau / electron_mass_c2;
>> 1160 }
>> 1161 // Relativistic
>> 1162 else
>> 1163 {
>> 1164 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
>> 1165 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
>> 1166 }
>> 1167
>> 1168 G4double v = std::sqrt(v2);
>> 1169 //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
>> 1170 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
>> 1171 G4double L2 = C2*std::pow(v,(D2));
>> 1172 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
>> 1173 G4double H2 = (A2/v2) + (B2/(v2*v2));
>> 1174 G4double F1 = L1+H1;
>> 1175 G4double F2 = (L2*H2)/(L2+H2);
>> 1176
>> 1177 // ZF. generalized & relativistic version
>> 1178 G4double maximumEnergy;
>> 1179
>> 1180 //---- maximum kinetic energy , non relativistic ------
>> 1181 if( (k/MeV)/(particle->GetPDGMass()/MeV) <= 0.1 )
>> 1182 {
>> 1183 maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k;
>> 1184 }
>> 1185 //---- relativistic -----------------------------------
>> 1186 else
>> 1187 {
>> 1188 G4double gamma = 1./sqrt(1.-beta2);
>> 1189 maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
>> 1190 (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) );
>> 1191 }
>> 1192
>> 1193 //either it is transfered energy or secondary electron energy ...
>> 1194 //maximumEnergy-=Bj_energy;
>> 1195
>> 1196 //-----------------------------------------------------
>> 1197 G4double wmax = maximumEnergy/Bj_energy;
>> 1198 G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
>> 1199 c=1./c; //!!!!!!!!!!! manual calculus leads to c=1/c
>> 1200 G4double randVal = G4UniformRand();
>> 1201 G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
>> 1202 proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
>> 1203 // proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws);
>> 1204 proposed_ws/= ( F1*c + F2*c - 2.*randVal );
>> 1205 proposed_ws*=Bj_energy;
>> 1206
>> 1207 return(proposed_ws);
609 1208
610 G4double r = Rh(kine, energyTransfer, slater <<
611 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. <<
612 return value; <<
613 } 1209 }
614 1210
615 //....oooOO0OOooo........oooOO0OOooo........oo 1211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
616 1212
617 G4double G4DNARuddIonisationExtendedModel::S_2 << 1213 G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t,
618 << 1214 G4double energyTransferred,
619 << 1215 G4double slaterEffectiveChg,
620 1216 G4double shellNumber)
621 { 1217 {
622 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4) << 1218 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
623 // Dingfelder, in Chattanooga 2005 proceedin << 1219 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
624 1220
625 G4double r = Rh(kine, energyTransfer, slater << 1221 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
626 G4double value = << 1222 G4double value = 1. - G4Exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
627 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * <<
628 1223
629 return value; << 1224 return value;
630 } 1225 }
631 1226
632 //....oooOO0OOooo........oooOO0OOooo........oo 1227 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
633 1228
634 G4double G4DNARuddIonisationExtendedModel::S_2 << 1229 G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t,
635 << 1230 G4double energyTransferred,
636 << 1231 G4double slaterEffectiveChg,
637 1232 G4double shellNumber)
638 { 1233 {
639 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^ << 1234 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
640 // Dingfelder, in Chattanooga 2005 proceedin << 1235 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
>> 1236
>> 1237 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
>> 1238 G4double value = 1. - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
641 1239
642 G4double r = Rh(kine, energyTransfer, slater << 1240 return value;
643 G4double value = <<
644 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3 <<
645 1241
646 return value; <<
647 } 1242 }
648 1243
649 //....oooOO0OOooo........oooOO0OOooo........oo 1244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
650 1245
651 G4double G4DNARuddIonisationExtendedModel::Rh( << 1246 G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t,
652 << 1247 G4double energyTransferred,
>> 1248 G4double slaterEffectiveChg,
>> 1249 G4double shellNumber)
653 { 1250 {
654 // The following values are provided by M. D << 1251 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
655 // Dingfelder, in Chattanooga 2005 proceedin << 1252 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
656 1253
657 G4double escaled = CLHEP::electron_mass_c2/f << 1254 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
658 const G4double H = 13.60569172 * CLHEP::eV; << 1255 G4double value = 1. - G4Exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
659 G4double value = 2.0*std::sqrt(escaled / H)* <<
660 1256
661 return value; << 1257 return value;
662 } 1258 }
663 1259
664 //....oooOO0OOooo........oooOO0OOooo........oo 1260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
665 1261
666 G4double << 1262 G4double G4DNARuddIonisationExtendedModel::R(G4double t,
667 G4DNARuddIonisationExtendedModel::CorrectionFa << 1263 G4double energyTransferred,
>> 1264 G4double slaterEffectiveChg,
>> 1265 G4double shellNumber)
668 { 1266 {
669 // ZF Shortened << 1267 // tElectron = m_electron / m_alpha * t
670 G4double res = 1.0; << 1268 // Dingfelder, in Chattanooga 2005 proceedings, p 4
671 if (shell < 4 && 0 != idx) { << 1269
672 const G4double ln10 = fGpow->logZ(10); << 1270 G4double tElectron = 0.511/3728. * t;
673 G4double x = 2.0*((G4Log(kine/CLHEP::eV)/l <<
674 // The following values are provided by M. 1271 // The following values are provided by M. Dingfelder (priv. comm)
675 res = 0.6/(1.0 + G4Exp(x)) + 0.9; << 1272 G4double H = 2.*13.60569172 * eV;
676 } << 1273 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
677 return res; << 1274
>> 1275 return value;
>> 1276 }
>> 1277
>> 1278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1279
>> 1280 G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
>> 1281 {
>> 1282 // ZF Shortened
>> 1283 G4DNAGenericIonsManager *instance;
>> 1284 instance = G4DNAGenericIonsManager::Instance();
>> 1285
>> 1286 if (particleDefinition == instance->GetIon("hydrogen") && shell < 4)
>> 1287 {
>> 1288 G4double value = (std::log10(k/eV)-4.2)/0.5;
>> 1289 // The following values are provided by M. Dingfelder (priv. comm)
>> 1290 return((0.6/(1+G4Exp(value))) + 0.9);
>> 1291 }
>> 1292 else
>> 1293 {
>> 1294 return(1.);
>> 1295 }
>> 1296 }
>> 1297
>> 1298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1299
>> 1300 G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k, const G4String& particle )
>> 1301 {
>> 1302
>> 1303 G4int level = 0;
>> 1304
>> 1305 // Retrieve data table corresponding to the current particle type
>> 1306
>> 1307 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 1308 pos = tableData.find(particle);
>> 1309
>> 1310 if (pos != tableData.end())
>> 1311 {
>> 1312 G4DNACrossSectionDataSet* table = pos->second;
>> 1313
>> 1314 if (table != 0)
>> 1315 {
>> 1316 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
>> 1317
>> 1318 const size_t n(table->NumberOfComponents());
>> 1319 size_t i(n);
>> 1320 G4double value = 0.;
>> 1321
>> 1322 while (i>0)
>> 1323 {
>> 1324 i--;
>> 1325 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
>> 1326
>> 1327 value += valuesBuffer[i];
>> 1328 }
>> 1329
>> 1330 value *= G4UniformRand();
>> 1331
>> 1332 i = n;
>> 1333
>> 1334 while (i > 0)
>> 1335 {
>> 1336 i--;
>> 1337
>> 1338 if (valuesBuffer[i] > value)
>> 1339 {
>> 1340 delete[] valuesBuffer;
>> 1341 return i;
>> 1342 }
>> 1343 value -= valuesBuffer[i];
>> 1344 }
>> 1345
>> 1346 if (valuesBuffer) delete[] valuesBuffer;
>> 1347
>> 1348 }
>> 1349 }
>> 1350 else
>> 1351 {
>> 1352 G4Exception("G4DNARuddIonisationExtendedModel::RandomSelect","em0002",
>> 1353 FatalException,"Model not applicable to particle type.");
>> 1354 }
>> 1355
>> 1356 return level;
678 } 1357 }
679 1358
680 //....oooOO0OOooo........oooOO0OOooo........oo 1359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1360
>> 1361 G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track )
>> 1362 {
>> 1363 G4double sigma = 0.;
>> 1364
>> 1365 const G4DynamicParticle* particle = track.GetDynamicParticle();
>> 1366 G4double k = particle->GetKineticEnergy();
>> 1367
>> 1368 G4double lowLim = 0;
>> 1369 G4double highLim = 0;
>> 1370
>> 1371 const G4String& particleName = particle->GetDefinition()->GetParticleName();
>> 1372
>> 1373 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
>> 1374 pos1 = lowEnergyLimit.find(particleName);
>> 1375
>> 1376 if (pos1 != lowEnergyLimit.end())
>> 1377 {
>> 1378 lowLim = pos1->second;
>> 1379 }
>> 1380
>> 1381 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
>> 1382 pos2 = highEnergyLimit.find(particleName);
>> 1383
>> 1384 if (pos2 != highEnergyLimit.end())
>> 1385 {
>> 1386 highLim = pos2->second;
>> 1387 }
>> 1388
>> 1389 if (k >= lowLim && k <= highLim)
>> 1390 {
>> 1391 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 1392 pos = tableData.find(particleName);
>> 1393
>> 1394 if (pos != tableData.end())
>> 1395 {
>> 1396 G4DNACrossSectionDataSet* table = pos->second;
>> 1397 if (table != 0)
>> 1398 {
>> 1399 sigma = table->FindValue(k);
>> 1400 }
>> 1401 }
>> 1402 else
>> 1403 {
>> 1404 G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002",
>> 1405 FatalException,"Model not applicable to particle type.");
>> 1406 }
>> 1407 }
>> 1408
>> 1409 return sigma;
>> 1410 }
>> 1411
>> 1412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1413
>> 1414 G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
>> 1415 {
>> 1416 return 0;
>> 1417 }
>> 1418
681 1419