source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlung.cc@ 968

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

update processes
File size: 15.6 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: G4PenelopeBremsstrahlung.cc,v 1.18 2006/06/29 19:40:35 gunter Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
28//
29// --------------------------------------------------------------
30//
31// File name: G4PenelopeBremsstrahlung
32//
33// Author: Luciano Pandola
34//
35// Creation date: February 2003
36//
37// Modifications:
38// 24.04.2003 V.Ivanchenko - Cut per region mfpt
39// 20.05.2003 MGP - Removed compilation warnings
40// Restored NotForced in GetMeanFreePath
41// 23.05.2003 MGP - Removed memory leak (fix in destructor)
42// 07.11.2003 L.Pandola - Bug fixed in LoadAngularData()
43// 11.11.2003 L.Pandola - Code review: use std::map for angular data
44// 01.06.2004 L.Pandola - StopButAlive for positrons on PostStepDoIt
45//
46//----------------------------------------------------------------
47
48#include "G4PenelopeBremsstrahlung.hh"
49#include "G4PenelopeBremsstrahlungContinuous.hh"
50#include "G4eBremsstrahlungSpectrum.hh"
51#include "G4BremsstrahlungCrossSectionHandler.hh"
52#include "G4VDataSetAlgorithm.hh"
53#include "G4LogLogInterpolation.hh"
54#include "G4VEMDataSet.hh"
55#include "G4EnergyLossTables.hh"
56#include "G4UnitsTable.hh"
57#include "G4Electron.hh"
58#include "G4Gamma.hh"
59#include "G4MaterialCutsCouple.hh"
60#include "G4DataVector.hh"
61#include "G4ProductionCutsTable.hh"
62#include "G4ProcessManager.hh"
63
64G4PenelopeBremsstrahlung::G4PenelopeBremsstrahlung(const G4String& nam)
65 : G4eLowEnergyLoss(nam),
66 crossSectionHandler(0),
67 theMeanFreePath(0),
68 energySpectrum(0)
69{
70 angularData = new std::map<G4int,G4PenelopeBremsstrahlungAngular*>;
71 cutForPhotons = 0.;
72 verboseLevel = 0;
73}
74
75
76G4PenelopeBremsstrahlung::~G4PenelopeBremsstrahlung()
77{
78 delete crossSectionHandler;
79 delete energySpectrum;
80 delete theMeanFreePath;
81 for (G4int Z=1;Z<100;Z++){
82 if (angularData->count(Z)) delete (angularData->find(Z)->second);
83 }
84 delete angularData;
85}
86
87
88void G4PenelopeBremsstrahlung::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
89{
90 if(verboseLevel > 0) {
91 G4cout << "G4PenelopeBremsstrahlung::BuildPhysicsTable start"
92 << G4endl;
93 }
94
95 cutForSecondaryPhotons.clear();
96 LoadAngularData();
97 if(verboseLevel > 0) {
98 G4cout << "G4PenelopeBremsstrahlung: Angular data loaded" << G4endl;
99 }
100 // Create and fill BremsstrahlungParameters once
101 if ( energySpectrum != 0 ) delete energySpectrum;
102 //grid of reduced energy bins for photons
103
104 G4DataVector eBins;
105
106 eBins.push_back(1.0e-12);
107 eBins.push_back(0.05);
108 eBins.push_back(0.075);
109 eBins.push_back(0.1);
110 eBins.push_back(0.125);
111 eBins.push_back(0.15);
112 eBins.push_back(0.2);
113 eBins.push_back(0.25);
114 eBins.push_back(0.3);
115 eBins.push_back(0.35);
116 eBins.push_back(0.40);
117 eBins.push_back(0.45);
118 eBins.push_back(0.50);
119 eBins.push_back(0.55);
120 eBins.push_back(0.60);
121 eBins.push_back(0.65);
122 eBins.push_back(0.70);
123 eBins.push_back(0.75);
124 eBins.push_back(0.80);
125 eBins.push_back(0.85);
126 eBins.push_back(0.90);
127 eBins.push_back(0.925);
128 eBins.push_back(0.95);
129 eBins.push_back(0.97);
130 eBins.push_back(0.99);
131 eBins.push_back(0.995);
132 eBins.push_back(0.999);
133 eBins.push_back(0.9995);
134 eBins.push_back(0.9999);
135 eBins.push_back(0.99995);
136 eBins.push_back(0.99999);
137 eBins.push_back(1.0);
138
139
140 const G4String dataName("/penelope/br-sp-pen.dat");
141 energySpectrum = new G4eBremsstrahlungSpectrum(eBins,dataName);
142
143 //the shape of the energy spectrum for positron is the same used for the electrons,
144 //as the differential cross section is scaled of a factor f(E,Z) which is independent
145 //on the energy of the gamma
146
147 if(verboseLevel > 0) {
148 G4cout << "G4PenelopeBremsstrahlungSpectrum is initialized"
149 << G4endl;
150 }
151
152 // Create and fill G4CrossSectionHandler once
153
154 if( crossSectionHandler != 0 ) delete crossSectionHandler;
155 G4VDataSetAlgorithm* interpolation = new G4LogLogInterpolation();
156 G4double lowKineticEnergy = GetLowerBoundEloss();
157 G4double highKineticEnergy = GetUpperBoundEloss();
158 G4int totBin = GetNbinEloss();
159 crossSectionHandler = new G4BremsstrahlungCrossSectionHandler(energySpectrum, interpolation);
160 crossSectionHandler->Initialise(0,lowKineticEnergy, highKineticEnergy, totBin);
161 if (&aParticleType==G4Electron::Electron())
162 {
163 crossSectionHandler->LoadShellData("brem/br-cs-");
164 }
165 else
166 {
167 crossSectionHandler->LoadShellData("penelope/br-cs-pos-"); //cross section for positrons
168 }
169
170 if (verboseLevel > 0) {
171 G4cout << GetProcessName()
172 << " is created; Cross section data: "
173 << G4endl;
174 crossSectionHandler->PrintData();
175 G4cout << "Parameters: "
176 << G4endl;
177 energySpectrum->PrintData();
178 }
179
180 // Build loss table for Bremsstrahlung
181
182 BuildLossTable(aParticleType);
183
184 if(verboseLevel > 0) {
185 G4cout << "The loss table is built"
186 << G4endl;
187 }
188
189 if (&aParticleType==G4Electron::Electron()) {
190
191 RecorderOfElectronProcess[CounterOfElectronProcess] = (*this).theLossTable;
192 CounterOfElectronProcess++;
193 PrintInfoDefinition();
194
195 } else {
196
197 RecorderOfPositronProcess[CounterOfPositronProcess] = (*this).theLossTable;
198 CounterOfPositronProcess++;
199 }
200 // Build mean free path data using cut values
201
202 if( theMeanFreePath != 0 ) delete theMeanFreePath;
203 theMeanFreePath = crossSectionHandler->
204 BuildMeanFreePathForMaterials(&cutForSecondaryPhotons);
205
206 if(verboseLevel > 0) {
207 G4cout << "The MeanFreePath table is built"
208 << G4endl;
209 }
210
211 // Build common DEDX table for all ionisation processes
212
213 BuildDEDXTable(aParticleType);
214
215 if(verboseLevel > 0) {
216 G4cout << "G4PenelopeBremsstrahlung::BuildPhysicsTable end"
217 << G4endl;
218 }
219}
220
221
222void G4PenelopeBremsstrahlung::BuildLossTable(const G4ParticleDefinition& aParticleType)
223{
224 // Build table for energy loss due to soft brems
225 // the tables are built for *MATERIALS* binning is taken from LowEnergyLoss
226 G4double lowKineticEnergy = GetLowerBoundEloss();
227 G4double highKineticEnergy = GetUpperBoundEloss();
228 size_t totBin = GetNbinEloss();
229
230 // create table
231
232 if (theLossTable) {
233 theLossTable->clearAndDestroy();
234 delete theLossTable;
235 }
236
237 const G4ProductionCutsTable* theCoupleTable=
238 G4ProductionCutsTable::GetProductionCutsTable();
239 size_t numOfCouples = theCoupleTable->GetTableSize();
240 theLossTable = new G4PhysicsTable(numOfCouples);
241
242
243 // Clean up the vector of cuts
244 cutForSecondaryPhotons.clear();
245
246 // Loop for materials
247 for (size_t j=0; j<numOfCouples; j++) {
248
249 // create physics vector and fill it
250 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(lowKineticEnergy,
251 highKineticEnergy,
252 totBin);
253
254 // get material parameters needed for the energy loss calculation
255 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
256 const G4Material* material= couple->GetMaterial();
257 // the cut cannot be below lowest limit
258 G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
259 tCut = std::min(highKineticEnergy, tCut);
260 cutForSecondaryPhotons.push_back(tCut);
261
262 const G4ElementVector* theElementVector = material->GetElementVector();
263 size_t NumberOfElements = material->GetNumberOfElements() ;
264 const G4double* theAtomicNumDensityVector =
265 material->GetAtomicNumDensityVector();
266 if(verboseLevel > 1) {
267 G4cout << "Energy loss for material # " << j
268 << " tCut(keV)= " << tCut/keV
269 << G4endl;
270 }
271
272 G4DataVector* ionloss = new G4DataVector();
273 for (size_t i = 0; i<totBin; i++) {
274 ionloss->push_back(0.0);
275 }
276
277 const G4String partName = aParticleType.GetParticleName();
278 // loop for elements in the material
279 for (size_t iel=0; iel<NumberOfElements; iel++ ) {
280 G4int Z = (G4int)((*theElementVector)[iel]->GetZ());
281 G4PenelopeBremsstrahlungContinuous* ContLoss = new G4PenelopeBremsstrahlungContinuous(Z,tCut,GetLowerBoundEloss(),
282 GetUpperBoundEloss(),
283 partName);
284 // the initialization must be repeated for each Z
285 // now comes the loop for the kinetic energy values
286 for (size_t k = 0; k<totBin; k++) {
287 G4double lowEdgeEnergy = aVector->GetLowEdgeEnergy(k);
288 // method for calcutation of continuous loss
289 (*ionloss)[k] += ContLoss->CalculateStopping(lowEdgeEnergy) * theAtomicNumDensityVector[iel];
290 //va chiamato una volta per ogni energia
291 //per ogni energia k, somma su tutti gli elementi del materiale
292 }
293 delete ContLoss;
294 }
295
296 for (size_t ibin = 0; ibin<totBin; ibin++) {
297 aVector->PutValue(ibin,(*ionloss)[ibin]); //un valore per ogni energia (somma sugli elementi del mate)
298 }
299 delete ionloss;
300 theLossTable->insert(aVector);
301 }
302}
303
304
305G4VParticleChange* G4PenelopeBremsstrahlung::PostStepDoIt(const G4Track& track,
306 const G4Step& step)
307{
308 aParticleChange.Initialize(track);
309
310 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
311 const G4Material* material = couple->GetMaterial();
312 G4double kineticEnergy = track.GetKineticEnergy();
313 G4int index = couple->GetIndex();
314 G4double tCut = cutForSecondaryPhotons[index];
315
316 // Control limits
317 if(tCut >= kineticEnergy)
318 return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
319
320 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
321 G4double tGamma = energySpectrum->SampleEnergy(Z, tCut, kineticEnergy, kineticEnergy);
322 //Check if the Z has been inserted in the map
323 if (!(angularData->count(Z))) {
324 G4String excep = "Not found the angular data for material " + material->GetName();
325 G4Exception(excep);
326 }
327
328 //Check if the loaded angular data are right
329 //G4cout << "Material Z: " << angularData->find(Z)->second->GetAtomicNumber() << " !!" << G4endl;
330
331 // Sample gamma angle (Z - axis along the parent particle).
332 G4double dirZ = angularData->find(Z)->second->ExtractCosTheta(kineticEnergy,tGamma);
333 G4double totalEnergy = kineticEnergy + electron_mass_c2;
334 G4double phi = twopi * G4UniformRand();
335 G4double sinTheta = std::sqrt(1. - dirZ*dirZ);
336 G4double dirX = sinTheta*std::cos(phi);
337 G4double dirY = sinTheta*std::sin(phi);
338
339 G4ThreeVector gammaDirection (dirX, dirY, dirZ);
340 G4ThreeVector electronDirection = track.GetMomentumDirection();
341
342 gammaDirection.rotateUz(electronDirection);
343
344 //
345 // Update the incident particle
346 //
347
348 G4double finalEnergy = kineticEnergy - tGamma;
349
350 // Kinematic problem
351 if (finalEnergy < 0.) {
352 tGamma += finalEnergy;
353 finalEnergy = 0.0;
354 }
355
356 G4double momentum = std::sqrt((totalEnergy + electron_mass_c2)*kineticEnergy);
357
358 G4double finalX = momentum*electronDirection.x() - tGamma*gammaDirection.x();
359 G4double finalY = momentum*electronDirection.y() - tGamma*gammaDirection.y();
360 G4double finalZ = momentum*electronDirection.z() - tGamma*gammaDirection.z();
361
362 aParticleChange.SetNumberOfSecondaries(1);
363 G4double norm = 1./std::sqrt(finalX*finalX + finalY*finalY + finalZ*finalZ);
364 aParticleChange.ProposeMomentumDirection(finalX*norm, finalY*norm, finalZ*norm);
365
366 const G4ParticleDefinition* particle = track.GetDefinition();
367
368 if (finalEnergy > 0.)
369 {
370 aParticleChange.ProposeEnergy(finalEnergy) ;
371 }
372 else
373 {
374 aParticleChange.ProposeEnergy(0.) ;
375 if (particle->GetProcessManager()->GetAtRestProcessVector()->size())
376 //In this case there is at least one AtRest process
377 {
378 aParticleChange.ProposeTrackStatus(fStopButAlive);
379 }
380 else
381 {
382 aParticleChange.ProposeTrackStatus(fStopAndKill);
383 }
384 }
385
386
387 // create G4DynamicParticle object for the gamma
388 G4DynamicParticle* aGamma= new G4DynamicParticle (G4Gamma::Gamma(),
389 gammaDirection, tGamma);
390 aParticleChange.AddSecondary(aGamma);
391
392 return G4VContinuousDiscreteProcess::PostStepDoIt(track, step);
393}
394
395
396void G4PenelopeBremsstrahlung::PrintInfoDefinition()
397{
398 G4String comments = "Total cross sections from EEDL database for electrons.";
399 comments += "\n Total cross section for positrons calculated from the electrons";
400 comments += "\n through an empirical scaling function.";
401 comments += "\n Gamma energy sampled from a data-driven histogram.";
402 comments += "\n Implementation of the continuous dE/dx part.";
403 comments += "\n It can be used for electrons and positrons";
404 comments += " in the energy range [250eV,100GeV].";
405 comments += "\n The process must work with G4PenelopeIonisation.";
406
407 G4cout << G4endl << GetProcessName() << ": " << comments << G4endl;
408}
409
410G4bool G4PenelopeBremsstrahlung::IsApplicable(const G4ParticleDefinition& particle)
411{
412 return ( (&particle == G4Electron::Electron()) || (&particle == G4Positron::Positron()) );
413}
414
415
416G4double G4PenelopeBremsstrahlung::GetMeanFreePath(const G4Track& track,
417 G4double, // previousStepSize
418 G4ForceCondition* cond)
419{
420 *cond = NotForced;
421 G4int index = (track.GetMaterialCutsCouple())->GetIndex();
422 const G4VEMDataSet* data = theMeanFreePath->GetComponent(index);
423 G4double meanFreePath = data->FindValue(track.GetKineticEnergy());
424 return meanFreePath;
425}
426
427void G4PenelopeBremsstrahlung::SetCutForLowEnSecPhotons(G4double cut)
428{
429 cutForPhotons = cut;
430}
431
432void G4PenelopeBremsstrahlung::LoadAngularData()
433{
434 const G4ProductionCutsTable* theCoupleTable=
435 G4ProductionCutsTable::GetProductionCutsTable();
436 size_t numOfCouples = theCoupleTable->GetTableSize();
437 angularData->clear();
438 for (size_t j=0; j<numOfCouples; j++) {
439 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
440 const G4Material* material= couple->GetMaterial();
441 const G4ElementVector* theElementVector = material->GetElementVector();
442 size_t NumberOfElements = material->GetNumberOfElements();
443 //loop for elements in the material
444 for (size_t iel=0; iel<NumberOfElements; iel++ ) {
445 G4int Z = (G4int)((*theElementVector)[iel]->GetZ());
446 //if the material is not present yet --> insert it in the map
447 if (!(angularData->count(Z))) {
448 angularData->insert(std::make_pair(Z,new G4PenelopeBremsstrahlungAngular(Z)));
449 //G4cout << "Loaded......... Z= " << Z << G4endl;
450 }
451 }
452 }
453}
454
Note: See TracBrowser for help on using the repository browser.