source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungModel.cc@ 992

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

fichier ajoutes
File size: 20.8 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: G4PenelopeBremsstrahlungModel.cc,v 1.2 2008/12/18 11:39:23 pandola Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
28//
29// Author: Luciano Pandola
30// --------
31// 05 Dec 2008 L Pandola Migration from process to model
32//
33#include "G4PenelopeBremsstrahlungModel.hh"
34#include "G4PenelopeBremsstrahlungContinuous.hh"
35#include "G4PenelopeBremsstrahlungAngular.hh"
36#include "G4eBremsstrahlungSpectrum.hh"
37#include "G4CrossSectionHandler.hh"
38#include "G4VEMDataSet.hh"
39#include "G4DataVector.hh"
40#include "G4Positron.hh"
41#include "G4Electron.hh"
42#include "G4Gamma.hh"
43#include "G4MaterialCutsCouple.hh"
44#include "G4ProductionCutsTable.hh"
45#include "G4ProcessManager.hh"
46#include "G4LogLogInterpolation.hh"
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50G4PenelopeBremsstrahlungModel::G4PenelopeBremsstrahlungModel(const G4ParticleDefinition*,
51 const G4String& nam)
52 :G4VEmModel(nam),isInitialised(false),energySpectrum(0),
53 angularData(0),stoppingPowerData(0),crossSectionHandler(0)
54{
55 fIntrinsicLowEnergyLimit = 100.0*eV;
56 fIntrinsicHighEnergyLimit = 100.0*GeV;
57 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
58 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
59 //
60
61 verboseLevel= 0;
62
63 // Verbosity scale:
64 // 0 = nothing
65 // 1 = warning for energy non-conservation
66 // 2 = details of energy budget
67 // 3 = calculation of cross sections, file openings, sampling of atoms
68 // 4 = entering in methods
69
70 //These vectors do not change when materials or cut change.
71 //Therefore I can read it at the constructor
72 angularData = new std::map<G4int,G4PenelopeBremsstrahlungAngular*>;
73
74 //These data do not depend on materials and cuts.
75 G4DataVector eBins;
76
77 eBins.push_back(1.0e-12);
78 eBins.push_back(0.05);
79 eBins.push_back(0.075);
80 eBins.push_back(0.1);
81 eBins.push_back(0.125);
82 eBins.push_back(0.15);
83 eBins.push_back(0.2);
84 eBins.push_back(0.25);
85 eBins.push_back(0.3);
86 eBins.push_back(0.35);
87 eBins.push_back(0.40);
88 eBins.push_back(0.45);
89 eBins.push_back(0.50);
90 eBins.push_back(0.55);
91 eBins.push_back(0.60);
92 eBins.push_back(0.65);
93 eBins.push_back(0.70);
94 eBins.push_back(0.75);
95 eBins.push_back(0.80);
96 eBins.push_back(0.85);
97 eBins.push_back(0.90);
98 eBins.push_back(0.925);
99 eBins.push_back(0.95);
100 eBins.push_back(0.97);
101 eBins.push_back(0.99);
102 eBins.push_back(0.995);
103 eBins.push_back(0.999);
104 eBins.push_back(0.9995);
105 eBins.push_back(0.9999);
106 eBins.push_back(0.99995);
107 eBins.push_back(0.99999);
108 eBins.push_back(1.0);
109
110 const G4String dataName("/penelope/br-sp-pen.dat");
111 energySpectrum = new G4eBremsstrahlungSpectrum(eBins,dataName);
112}
113
114//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115
116G4PenelopeBremsstrahlungModel::~G4PenelopeBremsstrahlungModel()
117{
118 if (crossSectionHandler)
119 delete crossSectionHandler;
120
121 if (energySpectrum)
122 delete energySpectrum;
123
124 if (angularData)
125 {
126 std::map <G4int,G4PenelopeBremsstrahlungAngular*>::iterator i;
127 for (i=angularData->begin();i != angularData->end();i++)
128 if (i->second) delete i->second;
129 delete angularData;
130 }
131
132 if (stoppingPowerData)
133 {
134 std::map <std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>::iterator j;
135 for (j=stoppingPowerData->begin();j != stoppingPowerData->end();j++)
136 if (j->second) delete j->second;
137 delete stoppingPowerData;
138 }
139}
140
141//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142
143void G4PenelopeBremsstrahlungModel::Initialise(const G4ParticleDefinition* particle,
144 const G4DataVector& cuts)
145{
146 if (verboseLevel > 3)
147 G4cout << "Calling G4PenelopeBremsstrahlungModel::Initialise()" << G4endl;
148
149 // Delete everything, but angular data (do not depend on cuts)
150 if (crossSectionHandler)
151 {
152 crossSectionHandler->Clear();
153 delete crossSectionHandler;
154 }
155
156 if (stoppingPowerData)
157 {
158 std::map <std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>::iterator j;
159 for (j=stoppingPowerData->begin();j != stoppingPowerData->end();j++)
160 if (j->second) delete j->second;
161 delete stoppingPowerData;
162 }
163 //Done.
164
165 if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
166 {
167 G4cout << "G4PenelopeBremsstrahlungModel: low energy limit increased from " <<
168 LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
169 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
170 }
171 if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
172 {
173 G4cout << "G4PenelopeBremsstrahlungModel: high energy limit decreased from " <<
174 HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
175 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
176 }
177
178 crossSectionHandler = new G4CrossSectionHandler();
179 crossSectionHandler->Clear();
180 //
181 if (particle==G4Electron::Electron())
182 crossSectionHandler->LoadData("brem/br-cs-");
183 else
184 crossSectionHandler->LoadData("penelope/br-cs-pos-"); //cross section for positrons
185
186 //This is used to retrieve cross section values later on
187 crossSectionHandler->BuildMeanFreePathForMaterials();
188
189 if (verboseLevel > 2)
190 G4cout << "Loaded cross section files for PenelopeBremsstrahlungModel" << G4endl;
191
192
193 G4cout << "Penelope Bremsstrahlung model is initialized " << G4endl
194 << "Energy range: "
195 << LowEnergyLimit() / keV << " keV - "
196 << HighEnergyLimit() / GeV << " GeV"
197 << G4endl;
198
199 //This has to be invoked AFTER the crossSectionHandler has been created,
200 //because it makes use of ComputeCrossSectionPerAtom()
201 InitialiseElementSelectors(particle,cuts);
202
203 if(isInitialised) return;
204
205 if(pParticleChange)
206 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
207 else
208 fParticleChange = new G4ParticleChangeForLoss();
209 isInitialised = true;
210}
211
212//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
213
214G4double G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
215 G4double kinEnergy,
216 G4double Z,
217 G4double,
218 G4double cutEnergy,
219 G4double)
220{
221 // Penelope model to calculate cross section for hard bremsstrahlung emission
222 // (gamma energy > cutEnergy).
223 //
224 // The total bremsstrahlung cross section is read from database, following data
225 // reported in the EEDL library
226 // D.E.Cullen et al., Report UCRL-50400 (Lawrence Livermore National Laboratory) (1989)
227 // The probability to have photon emission above a given threshold is calculated
228 // analytically using the differential cross section model dSigma/dW = F(x)/x, where
229 // W is the outgoing photon energy and x = W/E is the ratio of the photon energy to the
230 // incident energy. The function F(x) is tabulated (for all elements) using 32 points in x
231 // ranging from 1e-12 to 1. Data are derived from
232 // S.M.Seltzer and M.J.Berger, At.Data Nucl.Data Tables 35,345 (1986)
233 // Differential cross sections for electrons and positrons dSigma/dW are assumed to scale
234 // with a function S(Z,E) which does not depend on W; therefore, only overall cross section
235 // changes but not the shape of the photon energy spectrum.
236 //
237
238 if (verboseLevel > 3)
239 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopeBremsstrahlungModel" << G4endl;
240
241 G4int iZ = (G4int) Z;
242
243 if (!crossSectionHandler)
244 {
245 G4cout << "G4PenelopeBremsstrahlungModel::ComputeCrossSectionPerAtom" << G4endl;
246 G4cout << "The cross section handler is not correctly initialized" << G4endl;
247 G4Exception();
248 }
249 G4double totalCs = crossSectionHandler->FindValue(iZ,kinEnergy);
250 G4double cs = totalCs * energySpectrum->Probability(iZ,cutEnergy,kinEnergy,kinEnergy);
251
252 if (verboseLevel > 2)
253 {
254 G4cout << "Bremsstrahlung cross section at " << kinEnergy/MeV << " MeV for Z=" << Z <<
255 " and energy > " << cutEnergy/keV << " keV --> " << cs/barn << " barn" << G4endl;
256 G4cout << "Total bremsstrahlung cross section at " << kinEnergy/MeV << " MeV for Z=" <<
257 Z << " --> " << totalCs/barn << " barn" << G4endl;
258 }
259 return cs;
260
261}
262
263//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
264G4double G4PenelopeBremsstrahlungModel::ComputeDEDXPerVolume(const G4Material* theMaterial,
265 const G4ParticleDefinition* theParticle,
266 G4double kineticEnergy,
267 G4double cutEnergy)
268{
269 // Penelope model to calculate the stopping power (in [Energy]/[Length]) for soft
270 // bremsstrahlung emission (gamma energy < cutEnergy).
271 //
272 // The actual calculation is performed by the helper class
273 // G4PenelopeBremsstrahlungContinuous and its method CalculateStopping(). Notice:
274 // CalculateStopping() gives the stopping cross section, namely the first momentum of
275 // dSigma/dW, restricted to W < cut (W = gamma energy) This is dimensionally:
276 // [Energy]*[Surface]
277 // The calculation is performed by interpolation (in E = incident energy and
278 // x=W/E) from the tabulated data derived from
279 // M.J.Berger and S.M.Seltzer, Report NBSIR 82-2550 (National Bureau of Standards) (1982);
280 // for electrons.
281 // For positrons, dSigma/dW are assumed to scale with a function S(Z,E) with respect to electrons.
282 // An analytical approximation for the scaling function S(Z,E) is given in
283 // L.Kim et al., Phys.Rev.A 33,3002 (1986)
284 //
285 if (!stoppingPowerData)
286 stoppingPowerData = new std::map<std::pair<G4int,G4double>,
287 G4PenelopeBremsstrahlungContinuous*>;
288
289 const G4ElementVector* theElementVector = theMaterial->GetElementVector();
290 const G4double* theAtomicNumDensityVector = theMaterial->GetAtomicNumDensityVector();
291
292 G4double sPower = 0.0;
293
294 //Loop on the elements of the material
295 for (size_t iel=0;iel<theMaterial->GetNumberOfElements();iel++)
296 {
297 G4int iZ = (G4int) ((*theElementVector)[iel]->GetZ());
298 G4PenelopeBremsstrahlungContinuous* theContinuousCalculator =
299 GetStoppingPowerData(iZ,cutEnergy,theParticle);
300 sPower += theContinuousCalculator->CalculateStopping(kineticEnergy)*
301 theAtomicNumDensityVector[iel];
302 }
303
304 if (verboseLevel > 2)
305 {
306 G4cout << "Bremsstrahlung stopping power at " << kineticEnergy/MeV << " MeV for material " <<
307 theMaterial->GetName() << " and energy < " << cutEnergy/keV << " keV --> " <<
308 sPower/(keV/mm) << " keV/mm" << G4endl;
309 }
310
311 return sPower;
312}
313
314//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315
316void G4PenelopeBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
317 const G4MaterialCutsCouple* couple,
318 const G4DynamicParticle* aDynamicParticle,
319 G4double,G4double)
320{
321 // Penelope model to sample the final state for hard bremsstrahlung emission
322 // (gamma energy < cutEnergy).
323 // The energy distributionof the emitted photons is sampled according to the F(x)
324 // function tabulated in the database from
325 // S.M.Seltzer and M.J.Berger, At.Data Nucl.Data Tables 35,345 (1986)
326 // The database contains the function F(x) (32 points) for 57 energies of the
327 // incident electron between 1 keV and 100 GeV. For other primary energies,
328 // logarithmic interpolation is used to obtain the values of the function F(x).
329 // The double differential cross section dSigma/(dW dOmega), with W=photon energy,
330 // is described as
331 // dSigma/(dW dOmega) = dSigma/dW * p(Z,E,x,cosTheta)
332 // where the shape function p depends on atomic number, incident energy and x=W/E.
333 // Numerical values of the shape function, calculated by partial-waves methods, have been
334 // reported in
335 // L.Kissel et al., At.Data Nucl.Data.Tab. 28,381 (1983);
336 // for Z=2,8,13,47,79 and 92; E=1,5,10,50,100 and 500 keV; x=0,0.6,0.8 and 0.95. The
337 // function p(Z,E,x,cosTheta) is approximated by a Lorentz-dipole function as reported in
338 // Penelope - A Code System for Monte Carlo Simulation of Electron and Photon Transport,
339 // Workshop Proceedings Issy-les-Moulineaux, France, 5-7 November 2001, AEN-NEA;
340 // The analytical function contains two adjustable parameters that are obtained by fitting
341 // the complete solution from
342 // L.Kissel et al., At.Data Nucl.Data.Tab. 28,381 (1983);
343 // This allows the evaluation of p(Z,E,x,cosTheta) for any choice of Z, E and x. The actual
344 // sampling of cos(theta) is performed in the helper class
345 // G4PenelopeBremsstrahlungAngular, method ExtractCosTheta()
346 // Energy and direction of the primary particle are updated according to energy-momentum
347 // conservation. For positrons, it is sampled the same final state as for electrons.
348 //
349 if (verboseLevel > 3)
350 G4cout << "Calling SampleSecondaries() of G4PenelopeBremsstrahlungModel" << G4endl;
351
352 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
353 const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
354
355 if (kineticEnergy <= LowEnergyLimit())
356 {
357 fParticleChange->SetProposedKineticEnergy(0.);
358 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
359 //Check if there are AtRest processes
360 if (theParticle->GetProcessManager()->GetAtRestProcessVector()->size())
361 //In this case there is at least one AtRest process
362 fParticleChange->ProposeTrackStatus(fStopButAlive);
363 else
364 fParticleChange->ProposeTrackStatus(fStopAndKill);
365 return ;
366 }
367
368 G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
369 //This is the momentum
370 G4ThreeVector initialMomentum = aDynamicParticle->GetMomentum();
371
372 //One can use Vladimir's selector!
373 if (verboseLevel > 2)
374 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
375 // atom can be selected effitiantly if element selectors are initialised
376 const G4Element* anElement = SelectRandomAtom(couple,theParticle,kineticEnergy);
377 G4int iZ = (G4int) anElement->GetZ();
378 if (verboseLevel > 2)
379 G4cout << "Selected " << anElement->GetName() << G4endl;
380 //
381 //VERIFICA SE SERVE LA CROSS SECTION TABLE COME PER IONISATION.
382 G4double cutForLowEnergySecondaryParticles = 250.0*eV;
383 const G4ProductionCutsTable* theCoupleTable=
384 G4ProductionCutsTable::GetProductionCutsTable();
385 size_t indx = couple->GetIndex();
386 G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
387 //Production cut for gamma
388 cutG = std::max(cutForLowEnergySecondaryParticles,cutG);
389
390 //Sample gamma's energy according to the spectrum
391 G4double gammaEnergy = energySpectrum->SampleEnergy(iZ,cutG,kineticEnergy,kineticEnergy);
392
393 //Now sample cosTheta for the Gamma
394 G4double cosThetaPrimary = GetAngularDataForZ(iZ)->ExtractCosTheta(kineticEnergy,gammaEnergy);
395
396 G4double residualPrimaryEnergy = kineticEnergy-gammaEnergy;
397 if (residualPrimaryEnergy < 0)
398 {
399 //Ok we have a problem, all energy goes with the gamma
400 gammaEnergy += residualPrimaryEnergy;
401 residualPrimaryEnergy = 0.0;
402 }
403
404 //Get primary kinematics
405 G4double sinTheta = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
406 G4double phi = twopi * G4UniformRand();
407 G4ThreeVector gammaDirection1(sinTheta* std::cos(phi),
408 sinTheta* std::sin(phi),
409 cosThetaPrimary);
410
411 gammaDirection1.rotateUz(particleDirection0);
412
413 //Produce final state according to momentum conservation
414 G4ThreeVector particleDirection1 = initialMomentum - gammaEnergy*gammaDirection1;
415 particleDirection1.unit(); //normalize
416
417 //Update the primary particle
418 if (residualPrimaryEnergy > 0)
419 {
420 fParticleChange->ProposeMomentumDirection(particleDirection1);
421 fParticleChange->SetProposedKineticEnergy(residualPrimaryEnergy);
422 }
423 else
424 {
425 fParticleChange->ProposeMomentumDirection(particleDirection1);
426 fParticleChange->SetProposedKineticEnergy(0.*eV);
427 if (theParticle->GetProcessManager()->GetAtRestProcessVector()->size())
428 //In this case there is at least one AtRest process
429 fParticleChange->ProposeTrackStatus(fStopButAlive);
430 else
431 fParticleChange->ProposeTrackStatus(fStopAndKill);
432 }
433 fParticleChange->ProposeLocalEnergyDeposit(0.*eV);
434
435 //Now produce the photon
436 G4DynamicParticle* theGamma = new G4DynamicParticle(G4Gamma::Gamma(),
437 gammaDirection1,
438 gammaEnergy);
439 fvect->push_back(theGamma);
440
441 if (verboseLevel > 1)
442 {
443 G4cout << "-----------------------------------------------------------" << G4endl;
444 G4cout << "Energy balance from G4PenelopeBremsstrahlung" << G4endl;
445 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
446 G4cout << "-----------------------------------------------------------" << G4endl;
447 G4cout << "Outgoing primary energy: " << residualPrimaryEnergy/keV << " keV" << G4endl;
448 G4cout << "Bremsstrahlung photon " << gammaEnergy/keV << " keV" << G4endl;
449 G4cout << "Total final state: " << (residualPrimaryEnergy+gammaEnergy)/keV
450 << " keV" << G4endl;
451 G4cout << "-----------------------------------------------------------" << G4endl;
452 }
453 if (verboseLevel > 0)
454 {
455 G4double energyDiff = std::fabs(residualPrimaryEnergy+gammaEnergy-kineticEnergy);
456 if (energyDiff > 0.05*keV)
457 G4cout << "Warning from G4PenelopeBremsstrahlung: problem with energy conservation: " <<
458 (residualPrimaryEnergy+gammaEnergy)/keV <<
459 " keV (final) vs. " <<
460 kineticEnergy/keV << " keV (initial)" << G4endl;
461 }
462}
463
464//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465
466G4PenelopeBremsstrahlungAngular* G4PenelopeBremsstrahlungModel::GetAngularDataForZ(G4int iZ)
467{
468 if (!angularData)
469 angularData = new std::map<G4int,G4PenelopeBremsstrahlungAngular*>;
470
471 if (angularData->count(iZ)) //the material already exists
472 return angularData->find(iZ)->second;
473
474 //Otherwise create a new object, store it and return it
475 G4PenelopeBremsstrahlungAngular* theAngular = new G4PenelopeBremsstrahlungAngular(iZ);
476 angularData->insert(std::make_pair(iZ,theAngular));
477
478 if (angularData->count(iZ)) //the material should exist now
479 return angularData->find(iZ)->second;
480 else
481 {
482 G4Exception("Problem in G4PenelopeBremsstrahlungModel::GetAngularDataForZ()");
483 return 0;
484 }
485}
486
487//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
488
489G4PenelopeBremsstrahlungContinuous* G4PenelopeBremsstrahlungModel::GetStoppingPowerData(G4int iZ,G4double energyCut,
490 const G4ParticleDefinition*
491 theParticle)
492{
493 if (!stoppingPowerData)
494 stoppingPowerData = new std::map<std::pair<G4int,G4double>,G4PenelopeBremsstrahlungContinuous*>;
495
496 std::pair<G4int,G4double> theKey = std::make_pair(iZ,energyCut);
497
498 if (stoppingPowerData->count(theKey)) //the material already exists
499 return stoppingPowerData->find(theKey)->second;
500
501 //Otherwise create a new object, store it and return it
502 G4String theParticleName = theParticle->GetParticleName();
503 G4PenelopeBremsstrahlungContinuous* theContinuous = new
504 G4PenelopeBremsstrahlungContinuous(iZ,energyCut,LowEnergyLimit(),HighEnergyLimit(),theParticleName);
505 stoppingPowerData->insert(std::make_pair(theKey,theContinuous));
506
507 if (stoppingPowerData->count(theKey)) //the material should exist now
508 return stoppingPowerData->find(theKey)->second;
509 else
510 {
511 G4Exception("Problem in G4PenelopeBremsstrahlungModel::GetStoppingPowerData()");
512 return 0;
513 }
514}
Note: See TracBrowser for help on using the repository browser.