Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4DNARuddIonisationExtendedModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id$
27// GEANT4 tag $Name: $
28//
29
30
31// Modified by Z. Francis to handle HZE && inverse rudd function sampling 26-10-2010
32
35#include "G4SystemOfUnits.hh"
37#include "G4LossTableManager.hh"
40
41//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42
43using namespace std;
44
45//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46
48 const G4String& nam)
49 :G4VEmModel(nam),isInitialised(false)
50{
51 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
52 fpWaterDensity = 0;
53
54 slaterEffectiveCharge[0]=0.;
55 slaterEffectiveCharge[1]=0.;
56 slaterEffectiveCharge[2]=0.;
57 sCoefficient[0]=0.;
58 sCoefficient[1]=0.;
59 sCoefficient[2]=0.;
60
61 lowEnergyLimitForA[1] = 0 * eV;
62 lowEnergyLimitForA[2] = 0 * eV;
63 lowEnergyLimitForA[3] = 0 * eV;
64 lowEnergyLimitOfModelForA[1] = 100 * eV;
65 lowEnergyLimitOfModelForA[4] = 1 * keV;
66 lowEnergyLimitOfModelForA[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma
67 killBelowEnergyForA[1] = lowEnergyLimitOfModelForA[1];
68 killBelowEnergyForA[4] = lowEnergyLimitOfModelForA[4];
69 killBelowEnergyForA[5] = lowEnergyLimitOfModelForA[5];
70
71
72 verboseLevel= 0;
73 // Verbosity scale:
74 // 0 = nothing
75 // 1 = warning for energy non-conservation
76 // 2 = details of energy budget
77 // 3 = calculation of cross sections, file openings, sampling of atoms
78 // 4 = entering in methods
79
80 if( verboseLevel>0 )
81 {
82 G4cout << "Rudd ionisation model is constructed " << G4endl;
83 }
84
85 //Mark this model as "applicable" for atomic deexcitation
87 fAtomDeexcitation = 0;
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
94{
95 // Cross section
96
97 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
98 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
99 {
100 G4DNACrossSectionDataSet* table = pos->second;
101 delete table;
102 }
103
104 // The following removal is forbidden G4VEnergyLossModel takes care of deletion
105 // however coverity will signal this as an error
106 //if (fAtomDeexcitation) {delete fAtomDeexcitation;}
107
108}
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111
113 const G4DataVector& /*cuts*/)
114{
115
116 if (verboseLevel > 3)
117 G4cout << "Calling G4DNARuddIonisationExtendedModel::Initialise()" << G4endl;
118
119 // Energy limits
120
121 G4String fileProton("dna/sigma_ionisation_p_rudd");
122 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
123 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
124 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
125 G4String fileHelium("dna/sigma_ionisation_he_rudd");
126 G4String fileCarbon("dna/sigma_ionisation_c_rudd");
127 G4String fileNitrogen("dna/sigma_ionisation_n_rudd");
128 G4String fileOxygen("dna/sigma_ionisation_o_rudd");
129 G4String fileIron("dna/sigma_ionisation_fe_rudd");
130
131 G4DNAGenericIonsManager *instance;
134 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
135 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
136 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
137 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
138 G4ParticleDefinition* carbonDef = instance->GetIon("carbon");
139 G4ParticleDefinition* nitrogenDef = instance->GetIon("nitrogen");
140 G4ParticleDefinition* oxygenDef = instance->GetIon("oxygen");
141 G4ParticleDefinition* ironDef = instance->GetIon("iron");
142
143 G4String proton;
144 G4String hydrogen;
145 G4String alphaPlusPlus;
146 G4String alphaPlus;
147 G4String helium;
148 G4String carbon;
149 G4String nitrogen;
150 G4String oxygen;
151 G4String iron;
152
153 G4double scaleFactor = 1 * m*m;
154
155 // LIMITS AND DATA
156
157 proton = protonDef->GetParticleName();
158 tableFile[proton] = fileProton;
159 lowEnergyLimit[proton] = lowEnergyLimitForA[1];
160 highEnergyLimit[proton] = 500. * keV;
161
162 // Cross section
163
165 eV,
166 scaleFactor );
167 tableProton->LoadData(fileProton);
168 tableData[proton] = tableProton;
169
170 // **********************************************************************************************
171
172 hydrogen = hydrogenDef->GetParticleName();
173 tableFile[hydrogen] = fileHydrogen;
174
175 lowEnergyLimit[hydrogen] = lowEnergyLimitForA[1];
176 highEnergyLimit[hydrogen] = 100. * MeV;
177
178 // Cross section
179
181 eV,
182 scaleFactor );
183 tableHydrogen->LoadData(fileHydrogen);
184
185 tableData[hydrogen] = tableHydrogen;
186
187 // **********************************************************************************************
188
189 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
190 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
191
192 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForA[4];
193 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
194
195 // Cross section
196
198 eV,
199 scaleFactor );
200 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
201
202 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
203
204 // **********************************************************************************************
205
206 alphaPlus = alphaPlusDef->GetParticleName();
207 tableFile[alphaPlus] = fileAlphaPlus;
208
209 lowEnergyLimit[alphaPlus] = lowEnergyLimitForA[4];
210 highEnergyLimit[alphaPlus] = 400. * MeV;
211
212 // Cross section
213
215 eV,
216 scaleFactor );
217 tableAlphaPlus->LoadData(fileAlphaPlus);
218 tableData[alphaPlus] = tableAlphaPlus;
219
220 // **********************************************************************************************
221
222 helium = heliumDef->GetParticleName();
223 tableFile[helium] = fileHelium;
224
225 lowEnergyLimit[helium] = lowEnergyLimitForA[4];
226 highEnergyLimit[helium] = 400. * MeV;
227
228 // Cross section
229
231 eV,
232 scaleFactor );
233 tableHelium->LoadData(fileHelium);
234 tableData[helium] = tableHelium;
235
236 // **********************************************************************************************
237
238 carbon = carbonDef->GetParticleName();
239 tableFile[carbon] = fileCarbon;
240
241 lowEnergyLimit[carbon] = lowEnergyLimitForA[5] * particle->GetAtomicMass();
242 highEnergyLimit[carbon] = 1e6* particle->GetAtomicMass() * MeV;
243
244 // Cross section
245
247 eV,
248 scaleFactor );
249 tableCarbon->LoadData(fileCarbon);
250 tableData[carbon] = tableCarbon;
251
252 // **********************************************************************************************
253
254 oxygen = oxygenDef->GetParticleName();
255 tableFile[oxygen] = fileOxygen;
256
257 lowEnergyLimit[oxygen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
258 highEnergyLimit[oxygen] = 1e6* particle->GetAtomicMass()* MeV;
259
260 // Cross section
261
263 eV,
264 scaleFactor );
265 tableOxygen->LoadData(fileOxygen);
266 tableData[oxygen] = tableOxygen;
267
268 // **********************************************************************************************
269
270 nitrogen = nitrogenDef->GetParticleName();
271 tableFile[nitrogen] = fileNitrogen;
272
273 lowEnergyLimit[nitrogen] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
274 highEnergyLimit[nitrogen] = 1e6* particle->GetAtomicMass()* MeV;
275
276 // Cross section
277
279 eV,
280 scaleFactor );
281 tableNitrogen->LoadData(fileNitrogen);
282 tableData[nitrogen] = tableNitrogen;
283
284 // **********************************************************************************************
285
286 iron = ironDef->GetParticleName();
287 tableFile[iron] = fileIron;
288
289 lowEnergyLimit[iron] = lowEnergyLimitForA[5]* particle->GetAtomicMass();
290 highEnergyLimit[iron] = 1e6* particle->GetAtomicMass()* MeV;
291
292 // Cross section
293
295 eV,
296 scaleFactor );
297 tableIron->LoadData(fileIron);
298 tableData[iron] = tableIron;
299
300 // ZF Following lines can be replaced by:
301 SetLowEnergyLimit(lowEnergyLimit[particle->GetParticleName()]);
302 SetHighEnergyLimit(highEnergyLimit[particle->GetParticleName()]);
303 // at least for HZE
304
305 /*
306 if (particle==protonDef)
307 {
308 SetLowEnergyLimit(lowEnergyLimit[proton]);
309 SetHighEnergyLimit(highEnergyLimit[proton]);
310 }
311
312 if (particle==hydrogenDef)
313 {
314 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
315 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
316 }
317
318 if (particle==heliumDef)
319 {
320 SetLowEnergyLimit(lowEnergyLimit[helium]);
321 SetHighEnergyLimit(highEnergyLimit[helium]);
322 }
323
324 if (particle==alphaPlusDef)
325 {
326 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
327 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
328 }
329
330 if (particle==alphaPlusPlusDef)
331 {
332 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
333 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
334 }
335
336 if (particle==carbonDef)
337 {
338 SetLowEnergyLimit(lowEnergyLimit[carbon]);
339 SetHighEnergyLimit(highEnergyLimit[carbon]);
340 }
341
342 if (particle==oxygenDef)
343 {
344 SetLowEnergyLimit(lowEnergyLimit[oxygen]);
345 SetHighEnergyLimit(highEnergyLimit[oxygen]);
346 }*/
347
348 //----------------------------------------------------------------------
349
350 if( verboseLevel>0 )
351 {
352 G4cout << "Rudd ionisation model is initialized " << G4endl
353 << "Energy range: "
354 << LowEnergyLimit() / eV << " eV - "
355 << HighEnergyLimit() / keV << " keV for "
356 << particle->GetParticleName()
357 << G4endl;
358 }
359
360 // Initialize water density pointer
362
363 //
364
365 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
366
367 if (isInitialised) { return; }
369 isInitialised = true;
370}
371
372//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
373
375 const G4ParticleDefinition* particleDefinition,
376 G4double k,
377 G4double,
378 G4double)
379{
380 if (verboseLevel > 3)
381 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationExtendedModel" << G4endl;
382
383 // Calculate total cross section for model
384
385 G4DNAGenericIonsManager *instance;
387
388 if (
389 particleDefinition != G4Proton::ProtonDefinition()
390 &&
391 particleDefinition != instance->GetIon("hydrogen")
392 &&
393 particleDefinition != instance->GetIon("alpha++")
394 &&
395 particleDefinition != instance->GetIon("alpha+")
396 &&
397 particleDefinition != instance->GetIon("helium")
398 &&
399 particleDefinition != instance->GetIon("carbon")
400 &&
401 particleDefinition != instance->GetIon("nitrogen")
402 &&
403 particleDefinition != instance->GetIon("oxygen")
404 &&
405 particleDefinition != instance->GetIon("iron")
406 )
407
408 return 0;
409
410 G4double lowLim = 0;
411
412 if ( particleDefinition == G4Proton::ProtonDefinition()
413 || particleDefinition == instance->GetIon("hydrogen")
414 )
415
416 lowLim = lowEnergyLimitOfModelForA[1];
417
418 else if ( particleDefinition == instance->GetIon("alpha++")
419 || particleDefinition == instance->GetIon("alpha+")
420 || particleDefinition == instance->GetIon("helium")
421 )
422
423 lowLim = lowEnergyLimitOfModelForA[4];
424
425 else lowLim = lowEnergyLimitOfModelForA[5];
426
427 G4double highLim = 0;
428 G4double sigma=0;
429
430
431 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
432
433 if(waterDensity!= 0.0)
434// if (material == nistwater || material->GetBaseMaterial() == nistwater)
435 {
436 const G4String& particleName = particleDefinition->GetParticleName();
437
438 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
439 pos2 = highEnergyLimit.find(particleName);
440
441 if (pos2 != highEnergyLimit.end())
442 {
443 highLim = pos2->second;
444 }
445
446 if (k <= highLim)
447 {
448
449 //SI : XS must not be zero otherwise sampling of secondaries method ignored
450
451 if (k < lowLim) k = lowLim;
452
453 //
454
455 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
456 pos = tableData.find(particleName);
457
458 if (pos != tableData.end())
459 {
460 G4DNACrossSectionDataSet* table = pos->second;
461 if (table != 0)
462 {
463 sigma = table->FindValue(k);
464 }
465 }
466 else
467 {
468 G4Exception("G4DNARuddIonisationExtendedModel::CrossSectionPerVolume","em0002",
469 FatalException,"Model not applicable to particle type.");
470 }
471
472 } // if (k >= lowLim && k < highLim)
473
474 if (verboseLevel > 2)
475 {
476 G4cout << "__________________________________" << G4endl;
477 G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO START" << G4endl;
478 G4cout << "°°° Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
479 G4cout << "°°° Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
480 G4cout << "°°° Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
481 // G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
482 G4cout << "°°° G4DNARuddIonisationExtendedModel - XS INFO END" << G4endl;
483
484 }
485
486 } // if (waterMaterial)
487
488 return sigma*waterDensity;
489// return sigma*material->GetAtomicNumDensityVector()[1];
490
491}
492
493//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
494
495void G4DNARuddIonisationExtendedModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
496 const G4MaterialCutsCouple* /*couple*/,
497 const G4DynamicParticle* particle,
498 G4double,
499 G4double)
500{
501 if (verboseLevel > 3)
502 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationExtendedModel" << G4endl;
503
504 G4double lowLim = 0;
505 G4double highLim = 0;
506
507 // ZF: the following line summarizes the commented part
508 if(particle->GetDefinition()->GetAtomicMass() <= 4) lowLim = killBelowEnergyForA[particle->GetDefinition()->GetAtomicMass()];
509 else lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
510
511 /*
512 if(particle->GetDefinition()->GetAtomicMass() >= 5) lowLim = killBelowEnergyForA[5]*particle->GetDefinition()->GetAtomicMass();
513
514
515 if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
516 || particle->GetDefinition() == instance->GetIon("hydrogen")
517 )
518
519 lowLim = killBelowEnergyForA[1];
520
521 if ( particle->GetDefinition() == instance->GetIon("alpha++")
522 || particle->GetDefinition() == instance->GetIon("alpha+")
523 || particle->GetDefinition() == instance->GetIon("helium")
524 )
525
526 lowLim = killBelowEnergyForA[4];*/
527
528 G4double k = particle->GetKineticEnergy();
529
530 const G4String& particleName = particle->GetDefinition()->GetParticleName();
531
532 // SI - the following is useless since lowLim is already defined
533 /*
534 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
535 pos1 = lowEnergyLimit.find(particleName);
536
537 if (pos1 != lowEnergyLimit.end())
538 {
539 lowLim = pos1->second;
540 }
541 */
542
543 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
544 pos2 = highEnergyLimit.find(particleName);
545
546 if (pos2 != highEnergyLimit.end())highLim = pos2->second;
547
548 if (k >= lowLim && k < highLim)
549 {
550 G4ParticleDefinition* definition = particle->GetDefinition();
551 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
552 /*
553 G4double particleMass = definition->GetPDGMass();
554 G4double totalEnergy = k + particleMass;
555 G4double pSquare = k*(totalEnergy+particleMass);
556 G4double totalMomentum = std::sqrt(pSquare);
557 */
558
559 G4int ionizationShell = RandomSelect(k,particleName);
560
561 // sample deexcitation
562 // here we assume that H_{2}O electronic levels are the same of Oxigen.
563 // this can be considered true with a rough 10% error in energy on K-shell,
564
565 G4int secNumberInit = 0; // need to know at a certain point the enrgy of secondaries
566 G4int secNumberFinal = 0; // So I'll make the diference and then sum the energies
567 G4double bindingEnergy = 0;
568 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
569
570 if(fAtomDeexcitation) {
571 G4int Z = 8;
573
574 if (ionizationShell <5 && ionizationShell >1)
575 {
576 as = G4AtomicShellEnumerator(4-ionizationShell);
577 }
578 else if (ionizationShell <2)
579 {
581 }
582
583
584 // DEBUG
585 // if (ionizationShell == 4) {
586 //
587 // G4cout << "Z: " << Z << " as: " << as
588 // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
589 // G4cout << "Press <Enter> key to continue..." << G4endl;
590 // G4cin.ignore();
591 // }
592
593
594 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
595 secNumberInit = fvect->size();
596 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
597 secNumberFinal = fvect->size();
598 }
599
600
601 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
602
603 G4double cosTheta = 0.;
604 G4double phi = 0.;
605 RandomizeEjectedElectronDirection(definition, k,secondaryKinetic, cosTheta, phi, ionizationShell);
606
607 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta);
608 G4double dirX = sinTheta*std::cos(phi);
609 G4double dirY = sinTheta*std::sin(phi);
610 G4double dirZ = cosTheta;
611 G4ThreeVector deltaDirection(dirX,dirY,dirZ);
612 deltaDirection.rotateUz(primaryDirection);
613
614 // Ignored for ions on electrons
615 /*
616 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
617
618 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
619 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
620 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
621 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
622 finalPx /= finalMomentum;
623 finalPy /= finalMomentum;
624 finalPz /= finalMomentum;
625
626 G4ThreeVector direction;
627 direction.set(finalPx,finalPy,finalPz);
628
629 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
630 */
631
633 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
634 G4double deexSecEnergy = 0;
635 for (G4int j=secNumberInit; j < secNumberFinal; j++) {
636
637 deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
638
639 }
640
642 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
643
644 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic) ;
645 fvect->push_back(dp);
646
647 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
649 ionizationShell,
650 theIncomingTrack);
651 }
652
653 // SI - not useful since low energy of model is 0 eV
654
655 if (k < lowLim)
656 {
660 }
661
662}
663
664//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
665
666G4double G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
667 G4double k,
668 G4int shell)
669{
670 //-- Fast sampling method -----
671 G4double proposed_energy;
672 G4double random1;
673 G4double value_sampling;
674 G4double max1;
675
676 do
677 {
678 proposed_energy = ProposedSampledEnergy(particleDefinition, k, shell); // Proposed energy by inverse function sampling
679
680 max1=0.;
681
682 for(G4double en=0.; en<20.; en+=1.) if(RejectionFunction(particleDefinition, k, en, shell) > max1)
683 max1=RejectionFunction(particleDefinition, k, en, shell);
684
685 random1 = G4UniformRand()*max1;
686
687 value_sampling = RejectionFunction(particleDefinition, k, proposed_energy, shell);
688
689 } while(random1 > value_sampling);
690
691 return(proposed_energy);
692}
693
694//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
695
696
697void G4DNARuddIonisationExtendedModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
698 G4double k,
699 G4double secKinetic,
700 G4double & cosTheta,
701 G4double & phi,
702 G4int shell )
703{
704 G4double maxSecKinetic = 0.;
705 G4double maximumEnergyTransfer = 0.;
706
707 /* if (particleDefinition == G4Proton::ProtonDefinition()
708 || particleDefinition == instance->GetIon("hydrogen"))
709 {
710 if(k/MeV < 100.)maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
711 else {
712 G4double beta2 = 1.-(1.+k*);
713 maxSecKinetic =
714 }
715 }
716
717 if (particleDefinition == instance->GetIon("helium")
718 || particleDefinition == instance->GetIon("alpha+")
719 || particleDefinition == instance->GetIon("alpha++"))
720 {
721 maxSecKinetic = 4.* (0.511 / 3728) * k;
722 }
723
724 if (particleDefinition == G4Carbon::Carbon())
725 {
726 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k / 12.;
727 }*/
728
729 // ZF. generalized & relativistic version
730
731 if( (k/MeV)/(particleDefinition->GetPDGMass()/MeV) <= 0.1 )
732 {
733 maximumEnergyTransfer= 4.* (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
734 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
735 }
736 else
737 {
738 G4double approx_nuc_number = particleDefinition->GetPDGMass() / proton_mass_c2;
739 G4double en_per_nucleon = k/approx_nuc_number;
740 G4double beta2 = 1. - 1./pow( (1.+(en_per_nucleon/electron_mass_c2)*(electron_mass_c2/proton_mass_c2)), 2.);
741 G4double gamma = 1./sqrt(1.-beta2);
742 maximumEnergyTransfer = 2.*electron_mass_c2*(gamma*gamma-1.)/(1.+2.*gamma*(electron_mass_c2/particleDefinition->GetPDGMass())+pow(electron_mass_c2/particleDefinition->GetPDGMass(), 2.) );
743 maximumEnergyTransfer+=waterStructure.IonisationEnergy(shell);
744 }
745
746 maxSecKinetic = maximumEnergyTransfer-waterStructure.IonisationEnergy(shell);
747
748 phi = twopi * G4UniformRand();
749
750 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
751 else cosTheta = (2.*G4UniformRand())-1.;
752
753}
754
755//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
756
757G4double G4DNARuddIonisationExtendedModel::RejectionFunction(G4ParticleDefinition* particleDefinition,
758 G4double k,
759 G4double proposed_ws,
760 G4int ionizationLevelIndex)
761{
762 const G4int j=ionizationLevelIndex;
763 G4double Bj_energy, alphaConst;
764 G4double Ry = 13.6*eV;
765 const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
766
767 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
768
769 // Following values provided by M. Dingfelder (priv. comm)
770 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
771
772 if (j == 4)
773 {
774 alphaConst = 0.66;
775 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
776 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
777 //---
778 }
779 else
780 {
781 alphaConst = 0.64;
782 Bj_energy = Bj[ionizationLevelIndex];
783 }
784
785 G4double energyTransfer = proposed_ws + Bj_energy;
786 proposed_ws/=Bj_energy;
787 G4DNAGenericIonsManager *instance;
789 G4double tau = 0.;
790 G4double A_ion = 0.;
791 tau = (electron_mass_c2 / particleDefinition->GetPDGMass()) * k;
792 A_ion = particleDefinition->GetAtomicMass();
793
794 G4double v2;
795 G4double beta2;
796
797 if((tau/MeV)<5.447761194e-2)
798 {
799 v2 = tau / Bj_energy;
800 beta2 = 2.*tau / electron_mass_c2;
801 }
802 // Relativistic
803 else
804 {
805 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
806 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
807 }
808
809 G4double v = std::sqrt(v2);
810 G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
811 G4double rejection_term = 1.+std::exp(alphaConst*(proposed_ws - wc) / v);
812 rejection_term = (1./rejection_term)*CorrectionFactor(particleDefinition,k,ionizationLevelIndex) * Gj[j];
813 //* (S/Bj_energy) ; Not needed anymore
814
815 G4bool isHelium = false;
816
817 if ( particleDefinition == G4Proton::ProtonDefinition()
818 || particleDefinition == instance->GetIon("hydrogen")
819 )
820 {
821 return(rejection_term);
822 }
823
824 else if(particleDefinition->GetAtomicMass() > 4) // anything above Helium
825 {
826 G4double Z = particleDefinition->GetAtomicNumber();
827 G4double x = 100.*std::sqrt(beta2)/std::pow(Z,(2./3.));
828 G4double Zeffion = Z*(1.-std::exp(-1.316*x+0.112*x*x-0.0650*x*x*x));
829 rejection_term*=Zeffion*Zeffion;
830 }
831
832 else if (particleDefinition == instance->GetIon("alpha++") )
833 {
834 isHelium = true;
835 slaterEffectiveCharge[0]=0.;
836 slaterEffectiveCharge[1]=0.;
837 slaterEffectiveCharge[2]=0.;
838 sCoefficient[0]=0.;
839 sCoefficient[1]=0.;
840 sCoefficient[2]=0.;
841 }
842
843 else if (particleDefinition == instance->GetIon("alpha+") )
844 {
845 isHelium = true;
846 slaterEffectiveCharge[0]=2.0;
847 // The following values are provided by M. Dingfelder (priv. comm)
848 slaterEffectiveCharge[1]=2.0;
849 slaterEffectiveCharge[2]=2.0;
850 //
851 sCoefficient[0]=0.7;
852 sCoefficient[1]=0.15;
853 sCoefficient[2]=0.15;
854 }
855
856 else if (particleDefinition == instance->GetIon("helium") )
857 {
858 isHelium = true;
859 slaterEffectiveCharge[0]=1.7;
860 slaterEffectiveCharge[1]=1.15;
861 slaterEffectiveCharge[2]=1.15;
862 sCoefficient[0]=0.5;
863 sCoefficient[1]=0.25;
864 sCoefficient[2]=0.25;
865 }
866
867// if ( particleDefinition == instance->GetIon("helium")
868// || particleDefinition == instance->GetIon("alpha+")
869// || particleDefinition == instance->GetIon("alpha++")
870// )
871 if (isHelium)
872 {
873
874 G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
875
876 zEff -= ( sCoefficient[0] * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.) +
877 sCoefficient[1] * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.) +
878 sCoefficient[2] * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.) );
879
880 rejection_term*= zEff * zEff;
881 }
882
883 return (rejection_term);
884}
885
886
887//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
888
889
890G4double G4DNARuddIonisationExtendedModel::ProposedSampledEnergy(G4ParticleDefinition* particle,
891 G4double k,
892 G4int ionizationLevelIndex)
893{
894
895 const G4int j=ionizationLevelIndex;
896
897 G4double A1, B1, C1, D1, E1, A2, B2, C2, D2;
898 //G4double alphaConst ;
899 G4double Bj_energy;
900
901 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV}; //Ding Paper
902 // Following values provided by M. Dingfelder (priv. comm)
903
904 const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
905
906 if (j == 4)
907 {
908 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
909 A1 = 1.25;
910 B1 = 0.5;
911 C1 = 1.00;
912 D1 = 1.00;
913 E1 = 3.00;
914 A2 = 1.10;
915 B2 = 1.30;
916 C2 = 1.00;
917 D2 = 0.00;
918 //alphaConst = 0.66;
919 //---Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
920 Bj_energy = waterStructure.IonisationEnergy(ionizationLevelIndex);
921 //---
922 }
923 else
924 {
925 //Data For Liquid Water from Dingfelder (Protons in Water)
926 A1 = 1.02;
927 B1 = 82.0;
928 C1 = 0.45;
929 D1 = -0.80;
930 E1 = 0.38;
931 A2 = 1.07;
932 //B2 = 14.6; From Ding Paper
933 // Value provided by M. Dingfelder (priv. comm)
934 B2 = 11.6;
935 //
936 C2 = 0.60;
937 D2 = 0.04;
938 //alphaConst = 0.64;
939
940 Bj_energy = Bj[ionizationLevelIndex];
941 }
942
943 G4double tau = 0.;
944 G4double A_ion = 0.;
945 tau = (electron_mass_c2 / particle->GetPDGMass()) * k;
946
947 A_ion = particle->GetAtomicMass();
948
949 G4double v2;
950 G4double beta2;
951 if((tau/MeV)<5.447761194e-2)
952 {
953 v2 = tau / Bj_energy;
954 beta2 = 2.*tau / electron_mass_c2;
955 }
956 // Relativistic
957 else
958 {
959 v2 = (electron_mass_c2 / 2. / Bj_energy) * (1. - (1./ pow( (1.+ (tau/electron_mass_c2)),2) ));
960 beta2 =1. - 1./(1.+ (tau/electron_mass_c2/A_ion))/(1.+ (tau/electron_mass_c2/A_ion));
961 }
962
963 G4double v = std::sqrt(v2);
964 //G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj_energy));
965 G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
966 G4double L2 = C2*std::pow(v,(D2));
967 G4double H1 = (A1*std::log(1.+v2)) / (v2+(B1/v2));
968 G4double H2 = (A2/v2) + (B2/(v2*v2));
969 G4double F1 = L1+H1;
970 G4double F2 = (L2*H2)/(L2+H2);
971
972 // ZF. generalized & relativistic version
973 G4double maximumEnergy;
974
975 //---- maximum kinetic energy , non relativistic ------
976 if( (k/MeV)/(particle->GetPDGMass()/MeV) <= 0.1 )
977 {
978 maximumEnergy = 4.* (electron_mass_c2 / particle->GetPDGMass()) * k;
979 }
980 //---- relativistic -----------------------------------
981 else
982 {
983 G4double gamma = 1./sqrt(1.-beta2);
984 maximumEnergy = 2.*electron_mass_c2*(gamma*gamma-1.)/
985 (1.+2.*gamma*(electron_mass_c2/particle->GetPDGMass())+pow(electron_mass_c2/particle->GetPDGMass(), 2.) );
986 }
987
988 //either it is transfered energy or secondary electron energy ...
989 //maximumEnergy-=Bj_energy;
990
991 //-----------------------------------------------------
992 G4double wmax = maximumEnergy/Bj_energy;
993 G4double c = wmax*(F2*wmax+F1*(2.+wmax))/(2.*(1.+wmax)*(1.+wmax));
994 c=1./c; //!!!!!!!!!!! manual calculus leads to c=1/c
995 G4double randVal = G4UniformRand();
996 G4double proposed_ws = F1*F1*c*c + 2.*F2*c*randVal - 2.*F1*c*randVal;
997 proposed_ws = -F1*c+2.*randVal+std::sqrt(proposed_ws);
998 // proposed_ws = -F1*c+2.*randVal-std::sqrt(proposed_ws);
999 proposed_ws/= ( F1*c + F2*c - 2.*randVal );
1000 proposed_ws*=Bj_energy;
1001
1002 return(proposed_ws);
1003
1004}
1005
1006//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1007
1008G4double G4DNARuddIonisationExtendedModel::S_1s(G4double t,
1009 G4double energyTransferred,
1010 G4double slaterEffectiveChg,
1011 G4double shellNumber)
1012{
1013 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
1014 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
1015
1016 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1017 G4double value = 1. - std::exp(-2 * r) * ( ( 2. * r + 2. ) * r + 1. );
1018
1019 return value;
1020}
1021
1022//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1023
1024G4double G4DNARuddIonisationExtendedModel::S_2s(G4double t,
1025 G4double energyTransferred,
1026 G4double slaterEffectiveChg,
1027 G4double shellNumber)
1028{
1029 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
1030 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
1031
1032 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1033 G4double value = 1. - std::exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
1034
1035 return value;
1036
1037}
1038
1039//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1040
1041G4double G4DNARuddIonisationExtendedModel::S_2p(G4double t,
1042 G4double energyTransferred,
1043 G4double slaterEffectiveChg,
1044 G4double shellNumber)
1045{
1046 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
1047 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
1048
1049 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
1050 G4double value = 1. - std::exp(-2 * r) * (((( 2./3. * r + 4./3.) * r + 2.) * r + 2.) * r + 1.);
1051
1052 return value;
1053}
1054
1055//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1056
1057G4double G4DNARuddIonisationExtendedModel::R(G4double t,
1058 G4double energyTransferred,
1059 G4double slaterEffectiveChg,
1060 G4double shellNumber)
1061{
1062 // tElectron = m_electron / m_alpha * t
1063 // Dingfelder, in Chattanooga 2005 proceedings, p 4
1064
1065 G4double tElectron = 0.511/3728. * t;
1066 // The following values are provided by M. Dingfelder (priv. comm)
1067 G4double H = 2.*13.60569172 * eV;
1068 G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) * (slaterEffectiveChg/shellNumber);
1069
1070 return value;
1071}
1072
1073//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1074
1075G4double G4DNARuddIonisationExtendedModel::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k, G4int shell)
1076{
1077 // ZF Shortened
1078 G4DNAGenericIonsManager *instance;
1080
1081 if (particleDefinition == instance->GetIon("hydrogen") && shell < 4)
1082 {
1083 G4double value = (std::log10(k/eV)-4.2)/0.5;
1084 // The following values are provided by M. Dingfelder (priv. comm)
1085 return((0.6/(1+std::exp(value))) + 0.9);
1086 }
1087 else
1088 {
1089 return(1.);
1090 }
1091}
1092
1093//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1094
1095G4int G4DNARuddIonisationExtendedModel::RandomSelect(G4double k, const G4String& particle )
1096{
1097
1098 G4int level = 0;
1099
1100 // Retrieve data table corresponding to the current particle type
1101
1102 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1103 pos = tableData.find(particle);
1104
1105 if (pos != tableData.end())
1106 {
1107 G4DNACrossSectionDataSet* table = pos->second;
1108
1109 if (table != 0)
1110 {
1111 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1112
1113 const size_t n(table->NumberOfComponents());
1114 size_t i(n);
1115 G4double value = 0.;
1116
1117 while (i>0)
1118 {
1119 i--;
1120 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1121
1122 value += valuesBuffer[i];
1123 }
1124
1125 value *= G4UniformRand();
1126
1127 i = n;
1128
1129 while (i > 0)
1130 {
1131 i--;
1132
1133 if (valuesBuffer[i] > value)
1134 {
1135 delete[] valuesBuffer;
1136 return i;
1137 }
1138 value -= valuesBuffer[i];
1139 }
1140
1141 if (valuesBuffer) delete[] valuesBuffer;
1142
1143 }
1144 }
1145 else
1146 {
1147 G4Exception("G4DNARuddIonisationExtendedModel::RandomSelect","em0002",
1148 FatalException,"Model not applicable to particle type.");
1149 }
1150
1151 return level;
1152}
1153
1154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1155
1156G4double G4DNARuddIonisationExtendedModel::PartialCrossSection(const G4Track& track )
1157{
1158 G4double sigma = 0.;
1159
1160 const G4DynamicParticle* particle = track.GetDynamicParticle();
1161 G4double k = particle->GetKineticEnergy();
1162
1163 G4double lowLim = 0;
1164 G4double highLim = 0;
1165
1166 const G4String& particleName = particle->GetDefinition()->GetParticleName();
1167
1168 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
1169 pos1 = lowEnergyLimit.find(particleName);
1170
1171 if (pos1 != lowEnergyLimit.end())
1172 {
1173 lowLim = pos1->second;
1174 }
1175
1176 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
1177 pos2 = highEnergyLimit.find(particleName);
1178
1179 if (pos2 != highEnergyLimit.end())
1180 {
1181 highLim = pos2->second;
1182 }
1183
1184 if (k >= lowLim && k <= highLim)
1185 {
1186 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
1187 pos = tableData.find(particleName);
1188
1189 if (pos != tableData.end())
1190 {
1191 G4DNACrossSectionDataSet* table = pos->second;
1192 if (table != 0)
1193 {
1194 sigma = table->FindValue(k);
1195 }
1196 }
1197 else
1198 {
1199 G4Exception("G4DNARuddIonisationExtendedModel::PartialCrossSection","em0002",
1200 FatalException,"Model not applicable to particle type.");
1201 }
1202 }
1203
1204 return sigma;
1205}
1206
1207//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1208
1209G4double G4DNARuddIonisationExtendedModel::Sum(G4double /* energy */, const G4String& /* particle */)
1210{
1211 return 0;
1212}
1213
G4AtomicShellEnumerator
@ eIonizedMolecule
@ FatalException
@ fStopAndKill
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define C1
#define G4UniformRand()
Definition: Randomize.hh:53
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual size_t NumberOfComponents(void) const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
virtual G4bool LoadData(const G4String &argFileName)
static G4DNAGenericIonsManager * Instance(void)
G4ParticleDefinition * GetIon(const G4String &name)
static G4DNAMolecularMaterial * Instance()
const std::vector< double > * GetNumMolPerVolTableFor(const G4Material *) const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4DNARuddIonisationExtendedModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationExtendedModel")
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:94
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
size_t GetIndex() const
Definition: G4Material.hh:261
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:576
const G4Track * GetCurrentTrack() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4int GetAtomicNumber() const
G4int GetAtomicMass() const
G4double GetPDGCharge() const
const G4String & GetParticleName() const
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
const G4DynamicParticle * GetDynamicParticle() const
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:585
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:109
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:529
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:522
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:592
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:641
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41