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