source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNARuddIonisationExtendedModel.cc

Last change on this file was 1350, checked in by garnier, 15 years ago

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