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

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

update
File size: 15.5 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 $
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.