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

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
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-04-beta-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.