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

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

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