source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNADingfelderChargeDecreaseModel.cc@ 1005

Last change on this file since 1005 was 968, checked in by garnier, 17 years ago

fichier ajoutes
File size: 18.2 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: G4DNADingfelderChargeDecreaseModel.cc,v 1.3 2009/02/16 11:00:11 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
28//
29
30#include "G4DNADingfelderChargeDecreaseModel.hh"
31
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
33
34using namespace std;
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37
38G4DNADingfelderChargeDecreaseModel::G4DNADingfelderChargeDecreaseModel(const G4ParticleDefinition*,
39 const G4String& nam)
40:G4VEmModel(nam),isInitialised(false)
41{
42
43 verboseLevel= 0;
44 // Verbosity scale:
45 // 0 = nothing
46 // 1 = warning for energy non-conservation
47 // 2 = details of energy budget
48 // 3 = calculation of cross sections, file openings, sampling of atoms
49 // 4 = entering in methods
50
51 G4cout << "Dingfelder charge decrease model is constructed " << G4endl;
52
53}
54
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56
57G4DNADingfelderChargeDecreaseModel::~G4DNADingfelderChargeDecreaseModel()
58{}
59
60//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61
62void G4DNADingfelderChargeDecreaseModel::Initialise(const G4ParticleDefinition* particle,
63 const G4DataVector& /*cuts*/)
64{
65
66 if (verboseLevel > 3)
67 G4cout << "Calling G4DNADingfelderChargeDecreaseModel::Initialise()" << G4endl;
68
69 // Energy limits
70
71 G4DNAGenericIonsManager *instance;
72 instance = G4DNAGenericIonsManager::Instance();
73 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
74 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
75 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
76
77 G4String proton;
78 G4String alphaPlusPlus;
79 G4String alphaPlus;
80
81 if (protonDef != 0)
82 {
83 proton = protonDef->GetParticleName();
84 lowEnergyLimit[proton] = 1. * keV;
85 highEnergyLimit[proton] = 10. * MeV;
86 }
87 else
88 {
89 G4Exception("G4DNADingfelderChargeDecreaseModel::Initialise: proton is not defined");
90 }
91
92 if (alphaPlusPlusDef != 0)
93 {
94 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
95 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
96 highEnergyLimit[alphaPlusPlus] = 10. * MeV;
97 }
98 else
99 {
100 G4Exception("G4DNADingfelderChargeDecreaseModel::Initialise: alphaPlusPlus is not defined");
101 }
102
103 if (alphaPlusDef != 0)
104 {
105 alphaPlus = alphaPlusDef->GetParticleName();
106 lowEnergyLimit[alphaPlus] = 1. * keV;
107 highEnergyLimit[alphaPlus] = 10. * MeV;
108 }
109 else
110 {
111 G4Exception("G4DNADingfelderChargeDecreaseModel::Initialise: alphaPlus is not defined");
112 }
113
114 if (particle==protonDef)
115 {
116 SetLowEnergyLimit(lowEnergyLimit[proton]);
117 SetHighEnergyLimit(highEnergyLimit[proton]);
118 }
119
120 if (particle==alphaPlusPlusDef)
121 {
122 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
123 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
124 }
125
126 if (particle==alphaPlusDef)
127 {
128 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
129 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
130 }
131
132 // Final state
133
134 //PROTON
135 f0[0][0]=1.;
136 a0[0][0]=-0.180;
137 a1[0][0]=-3.600;
138 b0[0][0]=-18.22;
139 b1[0][0]=-1.997;
140 c0[0][0]=0.215;
141 d0[0][0]=3.550;
142 x0[0][0]=3.450;
143 x1[0][0]=5.251;
144
145 numberOfPartialCrossSections[0] = 1;
146
147 //ALPHA++
148 f0[0][1]=1.; a0[0][1]=0.95;
149 a1[0][1]=-2.75;
150 b0[0][1]=-23.00;
151 c0[0][1]=0.215;
152 d0[0][1]=2.95;
153 x0[0][1]=3.50;
154
155 f0[1][1]=1.;
156 a0[1][1]=0.95;
157 a1[1][1]=-2.75;
158 b0[1][1]=-23.73;
159 c0[1][1]=0.250;
160 d0[1][1]=3.55;
161 x0[1][1]=3.72;
162
163 x1[0][1]=-1.;
164 b1[0][1]=-1.;
165
166 x1[1][1]=-1.;
167 b1[1][1]=-1.;
168
169 numberOfPartialCrossSections[1] = 2;
170
171 // ALPHA+
172 f0[0][2]=1.;
173 a0[0][2]=0.65;
174 a1[0][2]=-2.75;
175 b0[0][2]=-21.81;
176 c0[0][2]=0.232;
177 d0[0][2]=2.95;
178 x0[0][2]=3.53;
179
180 x1[0][2]=-1.;
181 b1[0][2]=-1.;
182
183 numberOfPartialCrossSections[2] = 1;
184
185 //
186
187 G4cout << "Dingfelder charge decrease model is initialized " << G4endl
188 << "Energy range: "
189 << LowEnergyLimit() / keV << " keV - "
190 << HighEnergyLimit() / MeV << " MeV for "
191 << particle->GetParticleName()
192 << G4endl;
193
194 if(!isInitialised)
195 {
196 isInitialised = true;
197
198 if(pParticleChange)
199 fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
200 else
201 fParticleChangeForGamma = new G4ParticleChangeForGamma();
202 }
203
204 // InitialiseElementSelectors(particle,cuts);
205
206 // Test if water material
207
208 flagMaterialIsWater= false;
209 densityWater = 0;
210
211 const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
212
213 if(theCoupleTable)
214 {
215 G4int numOfCouples = theCoupleTable->GetTableSize();
216
217 if(numOfCouples>0)
218 {
219 for (G4int i=0; i<numOfCouples; i++)
220 {
221 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
222 const G4Material* material = couple->GetMaterial();
223
224 size_t j = material->GetNumberOfElements();
225 while (j>0)
226 {
227 j--;
228 const G4Element* element(material->GetElement(j));
229 if (element->GetZ() == 8.)
230 {
231 G4double density = material->GetAtomicNumDensityVector()[j];
232 if (density > 0.)
233 {
234 flagMaterialIsWater = true;
235 densityWater = density;
236
237 if (verboseLevel > 3)
238 G4cout << "Water material is found with density(cm^-3)=" << density/(cm*cm*cm) << G4endl;
239 }
240 }
241 }
242
243 }
244 } // if(numOfCouples>0)
245
246 } // if (theCoupleTable)
247
248}
249
250//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251
252G4double G4DNADingfelderChargeDecreaseModel::CrossSectionPerVolume(const G4Material*,
253 const G4ParticleDefinition* particleDefinition,
254 G4double k,
255 G4double,
256 G4double)
257{
258 if (verboseLevel > 3)
259 G4cout << "Calling CrossSectionPerVolume() of G4DNADingfelderChargeDecreaseModel" << G4endl;
260
261 // Calculate total cross section for model
262
263 G4DNAGenericIonsManager *instance;
264 instance = G4DNAGenericIonsManager::Instance();
265
266 if (
267 particleDefinition != G4Proton::ProtonDefinition()
268 &&
269 particleDefinition != instance->GetIon("alpha++")
270 &&
271 particleDefinition != instance->GetIon("alpha+")
272 )
273
274 return 0;
275
276 G4double lowLim = 0;
277 G4double highLim = 0;
278 G4double crossSection = 0.;
279
280 if (flagMaterialIsWater)
281 {
282 const G4String& particleName = particleDefinition->GetParticleName();
283
284 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
285 pos1 = lowEnergyLimit.find(particleName);
286
287 if (pos1 != lowEnergyLimit.end())
288 {
289 lowLim = pos1->second;
290 }
291
292 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
293 pos2 = highEnergyLimit.find(particleName);
294
295 if (pos2 != highEnergyLimit.end())
296 {
297 highLim = pos2->second;
298 }
299
300 if (k >= lowLim && k < highLim)
301 {
302 crossSection = Sum(k,particleDefinition);
303 }
304
305 if (verboseLevel > 3)
306 {
307 G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
308 G4cout << " - Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
309 G4cout << " - Cross section per water molecule (cm^-1)=" << crossSection*densityWater/(1./cm) << G4endl;
310 }
311
312 } // if (flagMaterialIsWater)
313
314 return crossSection*densityWater;
315
316}
317
318//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
319
320void G4DNADingfelderChargeDecreaseModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
321 const G4MaterialCutsCouple* /*couple*/,
322 const G4DynamicParticle* aDynamicParticle,
323 G4double,
324 G4double)
325{
326 if (verboseLevel > 3)
327 G4cout << "Calling SampleSecondaries() of G4DNADingfelderChargeDecreaseModel" << G4endl;
328
329 G4double inK = aDynamicParticle->GetKineticEnergy();
330
331 G4ParticleDefinition* definition = aDynamicParticle->GetDefinition();
332
333 G4int finalStateIndex = RandomSelect(inK,definition);
334
335 G4int n = NumberOfFinalStates(definition, finalStateIndex);
336 G4double waterBindingEnergy = WaterBindingEnergyConstant(definition, finalStateIndex);
337 G4double outgoingParticleBindingEnergy = OutgoingParticleBindingEnergyConstant(definition, finalStateIndex);
338
339 G4double outK = 0.;
340 if (definition==G4Proton::Proton())
341 outK = inK - n*(inK*electron_mass_c2/proton_mass_c2) - waterBindingEnergy + outgoingParticleBindingEnergy;
342 else
343 outK = inK - n*(inK*electron_mass_c2/(3728*MeV)) - waterBindingEnergy + outgoingParticleBindingEnergy;
344
345 if (outK<0)
346 {
347 G4String message;
348 message="G4DNADingfelderChargeDecreaseModel::SampleSecondaries: final kinetic energy is below 0! Process ";
349 G4Exception(message);
350 }
351
352 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
353 fParticleChangeForGamma->ProposeLocalEnergyDeposit(waterBindingEnergy);
354
355 G4DynamicParticle* dp = new G4DynamicParticle (OutgoingParticleDefinition(definition, finalStateIndex),
356 aDynamicParticle->GetMomentumDirection(),
357 outK) ;
358 fvect->push_back(dp);
359
360}
361
362//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
363
364G4int G4DNADingfelderChargeDecreaseModel::NumberOfFinalStates(G4ParticleDefinition* particleDefinition,
365 G4int finalStateIndex )
366
367{
368 if (particleDefinition == G4Proton::Proton()) return 1;
369
370 G4DNAGenericIonsManager*instance;
371 instance = G4DNAGenericIonsManager::Instance();
372
373 if (particleDefinition == instance->GetIon("alpha++") )
374 {
375 if (finalStateIndex==0) return 1;
376 return 2;
377 }
378
379 if (particleDefinition == instance->GetIon("alpha+") ) return 1;
380
381 return 0;
382}
383
384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385
386G4ParticleDefinition* G4DNADingfelderChargeDecreaseModel::OutgoingParticleDefinition (G4ParticleDefinition* particleDefinition,
387 G4int finalStateIndex)
388{
389 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
390
391 if (particleDefinition == G4Proton::Proton()) return instance->GetIon("hydrogen");
392
393 if (particleDefinition == instance->GetIon("alpha++") )
394 {
395 if (finalStateIndex == 0) return instance->GetIon("alpha+");
396 return instance->GetIon("helium");
397 }
398
399 if (particleDefinition == instance->GetIon("alpha+") ) return instance->GetIon("helium");
400
401 return 0;
402}
403
404//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
405
406G4double G4DNADingfelderChargeDecreaseModel::WaterBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
407 G4int finalStateIndex )
408{
409 // Ionization energy of first water shell
410 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
411 // W + 10.79 eV -> W+ + e-
412
413 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
414
415 if(particleDefinition == G4Proton::Proton()) return 10.79*eV;
416
417 if (particleDefinition == instance->GetIon("alpha++") )
418 {
419 // Binding energy for W+ -> W++ + e- 10.79 eV
420 // Binding energy for W -> W+ + e- 10.79 eV
421 //
422 // Ionization energy of first water shell
423 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
424
425 if (finalStateIndex == 0) return 10.79*eV;
426
427 return 10.79*2*eV;
428 }
429
430 if (particleDefinition == instance->GetIon("alpha+") )
431 {
432 // Binding energy for W+ -> W++ + e- 10.79 eV
433 // Binding energy for W -> W+ + e- 10.79 eV
434 //
435 // Ionization energy of first water shell
436 // Rad. Phys. Chem. 59 p.267 by Dingf. et al.
437
438 return 10.79*eV;
439 }
440
441 return 0;
442}
443
444//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
445
446G4double G4DNADingfelderChargeDecreaseModel::OutgoingParticleBindingEnergyConstant(G4ParticleDefinition* particleDefinition,
447 G4int finalStateIndex)
448{
449 G4DNAGenericIonsManager * instance(G4DNAGenericIonsManager::Instance());
450
451 if(particleDefinition == G4Proton::Proton()) return 13.6*eV;
452
453 if (particleDefinition == instance->GetIon("alpha++") )
454 {
455 // Binding energy for He+ -> He++ + e- 54.509 eV
456 // Binding energy for He -> He+ + e- 24.587 eV
457
458 if (finalStateIndex==0) return 54.509*eV;
459
460 return (54.509 + 24.587)*eV;
461 }
462
463 if (particleDefinition == instance->GetIon("alpha+") )
464 {
465 // Binding energy for He+ -> He++ + e- 54.509 eV
466 // Binding energy for He -> He+ + e- 24.587 eV
467
468 return 24.587*eV;
469 }
470
471 return 0;
472}
473
474//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
475
476G4double G4DNADingfelderChargeDecreaseModel::PartialCrossSection(G4double k, G4int index,
477 const G4ParticleDefinition* particleDefinition)
478{
479 G4int particleTypeIndex = 0;
480 G4DNAGenericIonsManager* instance;
481 instance = G4DNAGenericIonsManager::Instance();
482
483 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
484 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
485 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
486
487 //
488 // sigma(T) = f0 10 ^ y(log10(T/eV))
489 //
490 // / a0 x + b0 x < x0
491 // |
492 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
493 // |
494 // \ a1 x + b1 x >= x1
495 //
496 //
497 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
498 //
499 // f0 has been added to the code in order to manage partial (shell-dependent) cross sections (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1)
500 //
501 // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
502 // Inelastic-collision cross sections of liquid water for interactions of energetic proton
503 //
504
505 if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
506 {
507 //
508 // if x1 < x0 means that x1 and b1 will be calculated with the following formula (this piece of code is run on all alphas and not on protons)
509 //
510 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
511 //
512 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
513 //
514
515 x1[index][particleTypeIndex]=x0[index][particleTypeIndex] + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex])
516 / (c0[index][particleTypeIndex] * d0[index][particleTypeIndex]), 1. / (d0[index][particleTypeIndex] - 1.));
517 b1[index][particleTypeIndex]=(a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex]
518 + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x1[index][particleTypeIndex]
519 - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
520 }
521
522 G4double x(std::log10(k/eV));
523 G4double y;
524
525 if (x<x0[index][particleTypeIndex])
526 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
527 else if (x<x1[index][particleTypeIndex])
528 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex]
529 * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
530 else
531 y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
532
533 return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
534
535}
536
537G4int G4DNADingfelderChargeDecreaseModel::RandomSelect(G4double k,
538 const G4ParticleDefinition* particleDefinition)
539{
540 G4int particleTypeIndex = 0;
541 G4DNAGenericIonsManager *instance;
542 instance = G4DNAGenericIonsManager::Instance();
543
544 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex = 0;
545 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex = 1;
546 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex = 2;
547
548 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
549 G4double* values(new G4double[n]);
550 G4double value(0);
551 G4int i = n;
552
553 while (i>0)
554 {
555 i--;
556 values[i]=PartialCrossSection(k, i, particleDefinition);
557 value+=values[i];
558 }
559
560 value*=G4UniformRand();
561
562 i=n;
563 while (i>0)
564 {
565 i--;
566
567 if (values[i]>value)
568 break;
569
570 value-=values[i];
571 }
572
573 delete[] values;
574
575 return i;
576}
577
578//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
579
580G4double G4DNADingfelderChargeDecreaseModel::Sum(G4double k, const G4ParticleDefinition* particleDefinition)
581{
582 G4int particleTypeIndex = 0;
583 G4DNAGenericIonsManager* instance;
584 instance = G4DNAGenericIonsManager::Instance();
585
586 if (particleDefinition == G4Proton::ProtonDefinition()) particleTypeIndex=0;
587 if (particleDefinition == instance->GetIon("alpha++")) particleTypeIndex=1;
588 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=2;
589
590 G4double totalCrossSection = 0.;
591
592 for (G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
593 {
594 totalCrossSection += PartialCrossSection(k,i,particleDefinition);
595 }
596 return totalCrossSection;
597}
598
Note: See TracBrowser for help on using the repository browser.