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