source: trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc @ 1015

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

update to geant4.9.2
File size: 18.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: G4VEmProcess.cc,v 1.60 2008/10/17 14:46:16 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:     G4VEmProcess
35//
36// Author:        Vladimir Ivanchenko on base of Laszlo Urban code
37//
38// Creation date: 01.10.2003
39//
40// Modifications:
41// 30-06-04 make it to be pure discrete process (V.Ivanchenko)
42// 30-09-08 optimise integral option (V.Ivanchenko)
43// 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
44// 11-03-05 Shift verbose level by 1, add applyCuts and killPrimary flags (VI)
45// 14-03-05 Update logic PostStepDoIt (V.Ivanchenko)
46// 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
47// 18-04-05 Use G4ParticleChangeForGamma (V.Ivanchenko)
48// 25-07-05 Add protection: integral mode only for charged particles (VI)
49// 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
50// 11-01-06 add A to parameters of ComputeCrossSectionPerAtom (VI)
51// 12-09-06 add SetModel() (mma)
52// 12-04-07 remove double call to Clear model manager (V.Ivanchenko)
53// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
54//
55// Class Description:
56//
57
58// -------------------------------------------------------------------
59//
60//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62
63#include "G4VEmProcess.hh"
64#include "G4LossTableManager.hh"
65#include "G4Step.hh"
66#include "G4ParticleDefinition.hh"
67#include "G4VEmModel.hh"
68#include "G4DataVector.hh"
69#include "G4PhysicsTable.hh"
70#include "G4PhysicsVector.hh"
71#include "G4PhysicsLogVector.hh"
72#include "G4VParticleChange.hh"
73#include "G4ProductionCutsTable.hh"
74#include "G4Region.hh"
75#include "G4RegionStore.hh"
76#include "G4Gamma.hh"
77#include "G4Electron.hh"
78#include "G4Positron.hh"
79#include "G4PhysicsTableHelper.hh"
80
81//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82
83G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
84  G4VDiscreteProcess(name, type),
85  secondaryParticle(0),
86  buildLambdaTable(true),
87  theLambdaTable(0),
88  theEnergyOfCrossSectionMax(0),
89  theCrossSectionMax(0),
90  integral(false),
91  applyCuts(false),
92  startFromNull(true),
93  nRegions(0),
94  selectedModel(0),
95  particle(0),
96  currentCouple(0)
97{
98  SetVerboseLevel(1);
99
100  // Size of tables assuming spline
101  minKinEnergy = 0.1*keV;
102  maxKinEnergy = 100.0*TeV;
103  nLambdaBins  = 84;
104
105  // default lambda factor
106  lambdaFactor  = 0.8;
107
108  // default limit on polar angle
109  polarAngleLimit = 0.0;
110
111  // particle types
112  theGamma     = G4Gamma::Gamma();
113  theElectron  = G4Electron::Electron();
114  thePositron  = G4Positron::Positron();
115
116  pParticleChange = &fParticleChange;
117  secParticles.reserve(5);
118
119  modelManager = new G4EmModelManager();
120  (G4LossTableManager::Instance())->Register(this);
121}
122
123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124
125G4VEmProcess::~G4VEmProcess()
126{
127  if(1 < verboseLevel) 
128    G4cout << "G4VEmProcess destruct " << GetProcessName() 
129           << G4endl;
130  Clear();
131  if(theLambdaTable) {
132    theLambdaTable->clearAndDestroy();
133    delete theLambdaTable;
134  }
135  delete modelManager;
136  (G4LossTableManager::Instance())->DeRegister(this);
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140
141void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
142{
143  if(!particle) particle = &part;
144  if(1 < verboseLevel) {
145    G4cout << "G4VEmProcess::PreparePhysicsTable() for "
146           << GetProcessName()
147           << " and particle " << part.GetParticleName()
148           << " local particle " << particle->GetParticleName() 
149           << G4endl;
150  }
151
152  if(particle == &part) {
153    Clear();
154    InitialiseProcess(particle);
155    theCuts = modelManager->Initialise(particle,secondaryParticle,2.,verboseLevel);
156    const G4ProductionCutsTable* theCoupleTable=
157          G4ProductionCutsTable::GetProductionCutsTable();
158    theCutsGamma    = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
159    theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
160    theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
161    if(buildLambdaTable)
162      theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
163  }
164}
165
166//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
167
168void G4VEmProcess::Clear()
169{
170  if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax;
171  if(theCrossSectionMax) delete [] theCrossSectionMax;
172  theEnergyOfCrossSectionMax = 0;
173  theCrossSectionMax = 0;
174  currentCouple = 0;
175  preStepLambda = 0.0;
176  mfpKinEnergy  = DBL_MAX;
177}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180
181void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
182{
183  if(1 < verboseLevel) {
184    G4cout << "G4VEmProcess::BuildPhysicsTable() for "
185           << GetProcessName()
186           << " and particle " << part.GetParticleName()
187           << " buildLambdaTable= " << buildLambdaTable
188           << G4endl;
189  }
190
191  if(buildLambdaTable) {
192    BuildLambdaTable();
193    FindLambdaMax();
194  }
195  if(0 < verboseLevel) PrintInfoDefinition();
196
197  if(1 < verboseLevel) {
198    G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
199           << GetProcessName()
200           << " and particle " << part.GetParticleName()
201           << G4endl;
202  }
203}
204
205//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206
207void G4VEmProcess::BuildLambdaTable()
208{
209  if(1 < verboseLevel) {
210    G4cout << "G4EmProcess::BuildLambdaTable() for process "
211           << GetProcessName() << " and particle "
212           << particle->GetParticleName()
213           << G4endl;
214  }
215
216  // Access to materials
217  const G4ProductionCutsTable* theCoupleTable=
218        G4ProductionCutsTable::GetProductionCutsTable();
219  size_t numOfCouples = theCoupleTable->GetTableSize();
220  for(size_t i=0; i<numOfCouples; i++) {
221
222    if (theLambdaTable->GetFlag(i)) {
223
224      // create physics vector and fill it
225      const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
226      G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
227      modelManager->FillLambdaVector(aVector, couple, startFromNull);
228      G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
229    }
230  }
231
232  if(1 < verboseLevel) {
233    G4cout << "Lambda table is built for "
234           << particle->GetParticleName()
235           << G4endl;
236    if(2 < verboseLevel) {
237      G4cout << *theLambdaTable << G4endl;
238    }
239  }
240}
241
242//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
243
244G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
245                             const G4Track& track,
246                             G4double   previousStepSize,
247                             G4ForceCondition* condition)
248{
249  // condition is set to "Not Forced"
250  *condition = NotForced;
251  G4double x = DBL_MAX;
252  if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
253  InitialiseStep(track);
254
255  if(preStepKinEnergy < mfpKinEnergy) {
256    if (integral) ComputeIntegralLambda(preStepKinEnergy);
257    else  preStepLambda = GetCurrentLambda(preStepKinEnergy);
258    if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
259  }
260
261  // non-zero cross section
262  if(preStepLambda > DBL_MIN) { 
263    if (theNumberOfInteractionLengthLeft < 0.0) {
264      // beggining of tracking (or just after DoIt of this process)
265      ResetNumberOfInteractionLengthLeft();
266    } else if(currentInteractionLength < DBL_MAX) {
267      // subtract NumberOfInteractionLengthLeft
268      SubtractNumberOfInteractionLengthLeft(previousStepSize);
269      if(theNumberOfInteractionLengthLeft < 0.)
270        theNumberOfInteractionLengthLeft = perMillion;
271    }
272
273    // get mean free path and step limit
274    currentInteractionLength = 1.0/preStepLambda;
275    x = theNumberOfInteractionLengthLeft * currentInteractionLength;
276#ifdef G4VERBOSE
277    if (verboseLevel>2){
278      G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
279      G4cout << "[ " << GetProcessName() << "]" << G4endl; 
280      G4cout << " for " << particle->GetParticleName() 
281             << " in Material  " <<  currentMaterial->GetName()
282             << " Ekin(MeV)= " << preStepKinEnergy/MeV
283             <<G4endl;
284      G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" 
285             << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
286    }
287#endif
288
289    // zero cross section case
290  } else {
291    if(theNumberOfInteractionLengthLeft > DBL_MIN && 
292       currentInteractionLength < DBL_MAX) {
293
294      // subtract NumberOfInteractionLengthLeft
295      SubtractNumberOfInteractionLengthLeft(previousStepSize);
296      if(theNumberOfInteractionLengthLeft < 0.)
297        theNumberOfInteractionLengthLeft = perMillion;
298    }
299    currentInteractionLength = DBL_MAX;
300  }
301  return x;
302}
303
304//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
305
306G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
307                                       G4double,
308                                       G4ForceCondition* condition)
309{
310  *condition = NotForced;
311  return G4VEmProcess::MeanFreePath(track);
312}
313
314//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315
316G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
317                                              const G4Step&)
318{
319  fParticleChange.InitializeForPostStep(track);
320
321  // Do not make anything if particle is stopped, the annihilation then
322  // should be performed by the AtRestDoIt!
323  if (track.GetTrackStatus() == fStopButAlive) return &fParticleChange;
324
325  G4double finalT = track.GetKineticEnergy();
326
327  // Integral approach
328  if (integral) {
329    G4double lx = GetLambda(finalT, currentCouple);
330    if(preStepLambda<lx && 1 < verboseLevel) {
331      G4cout << "WARING: for " << particle->GetParticleName() 
332             << " and " << GetProcessName()
333             << " E(MeV)= " << finalT/MeV
334             << " preLambda= " << preStepLambda << " < " << lx << " (postLambda) "
335             << G4endl; 
336    }
337
338    if(preStepLambda*G4UniformRand() > lx) {
339      ClearNumberOfInteractionLengthLeft();
340      return &fParticleChange;
341    }
342  }
343
344  G4VEmModel* currentModel = SelectModel(finalT);
345
346  /* 
347  if(0 < verboseLevel) {
348    G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
349           << finalT/MeV
350           << " MeV; model= (" << currentModel->LowEnergyLimit()
351           << ", " <<  currentModel->HighEnergyLimit() << ")"
352           << G4endl;
353  }
354  */
355
356 
357  // sample secondaries
358  secParticles.clear();
359  currentModel->SampleSecondaries(&secParticles, 
360                                  currentCouple, 
361                                  track.GetDynamicParticle(),
362                                  (*theCuts)[currentMaterialIndex]);
363
364  // save secondaries
365  G4int num = secParticles.size();
366  if(num > 0) {
367
368    fParticleChange.SetNumberOfSecondaries(num);
369    G4double edep = fParticleChange.GetLocalEnergyDeposit();
370     
371    for (G4int i=0; i<num; i++) {
372      G4DynamicParticle* dp = secParticles[i];
373      const G4ParticleDefinition* p = dp->GetDefinition();
374      G4double e = dp->GetKineticEnergy();
375      G4bool good = true;
376      if(applyCuts) {
377        if (p == theGamma) {
378          if (e < (*theCutsGamma)[currentMaterialIndex]) good = false;
379
380        } else if (p == theElectron) {
381          if (e < (*theCutsElectron)[currentMaterialIndex]) good = false;
382
383        } else if (p == thePositron) {
384          if (electron_mass_c2 < (*theCutsGamma)[currentMaterialIndex] &&
385              e < (*theCutsPositron)[currentMaterialIndex]) {
386            good = false;
387            e += 2.0*electron_mass_c2;
388          }
389        }
390        if(!good) {
391          delete dp;
392          edep += e;
393        }
394      }
395      if (good) fParticleChange.AddSecondary(dp);
396    } 
397    fParticleChange.ProposeLocalEnergyDeposit(edep);
398  }
399
400  ClearNumberOfInteractionLengthLeft();
401  return &fParticleChange;
402}
403
404//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
405
406void G4VEmProcess::PrintInfoDefinition()
407{
408  if(verboseLevel > 0) {
409    G4cout << G4endl << GetProcessName() << ":   for  "
410           << particle->GetParticleName();
411    if(integral) G4cout << ", integral: 1 ";
412    if(applyCuts) G4cout << ", applyCuts: 1 ";
413    G4cout << "    SubType= " << GetProcessSubType() << G4endl;
414    if(buildLambdaTable) {
415      G4cout << "      Lambda tables from "
416             << G4BestUnit(minKinEnergy,"Energy") 
417             << " to "
418             << G4BestUnit(maxKinEnergy,"Energy")
419             << " in " << nLambdaBins << " bins, spline: " 
420             << (G4LossTableManager::Instance())->SplineFlag()
421             << G4endl;
422    }
423    PrintInfo();
424    modelManager->DumpModelList(verboseLevel);
425  }
426
427  if(verboseLevel > 2 && buildLambdaTable) {
428    G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
429    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
430  }
431}
432
433//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
434
435G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
436                                             const G4MaterialCutsCouple* couple)
437{
438  // Cross section per atom is calculated
439  DefineMaterial(couple);
440  G4double cross = 0.0;
441  G4bool b;
442  if(theLambdaTable) {
443    cross = (((*theLambdaTable)[currentMaterialIndex])->
444                           GetValue(kineticEnergy, b));
445  } else {
446    G4VEmModel* model = SelectModel(kineticEnergy);
447    cross = 
448      model->CrossSectionPerVolume(currentMaterial,particle,kineticEnergy);
449  }
450
451  return cross;
452}
453
454//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
455
456G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
457                                       const G4String& directory,
458                                             G4bool ascii)
459{
460  G4bool yes = true;
461
462  if ( theLambdaTable && part == particle) {
463    const G4String name = 
464      GetPhysicsTableFileName(part,directory,"Lambda",ascii);
465    yes = theLambdaTable->StorePhysicsTable(name,ascii);
466
467    if ( yes ) {
468      G4cout << "Physics tables are stored for " << particle->GetParticleName()
469             << " and process " << GetProcessName()
470             << " in the directory <" << directory
471             << "> " << G4endl;
472    } else {
473      G4cout << "Fail to store Physics Tables for " 
474             << particle->GetParticleName()
475             << " and process " << GetProcessName()
476             << " in the directory <" << directory
477             << "> " << G4endl;
478    }
479  }
480  return yes;
481}
482
483//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
484
485G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
486                                          const G4String& directory,
487                                                G4bool ascii)
488{
489  if(1 < verboseLevel) {
490    G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
491           << part->GetParticleName() << " and process "
492           << GetProcessName() << G4endl;
493  }
494  G4bool yes = true;
495
496  if(!buildLambdaTable || particle != part) return yes;
497
498  const G4String particleName = part->GetParticleName();
499  G4String filename;
500
501  filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
502  yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
503                                                   filename,ascii);
504  if ( yes ) {
505    if (0 < verboseLevel) {
506      G4cout << "Lambda table for " << particleName << " is Retrieved from <"
507             << filename << ">"
508             << G4endl;
509    }
510    if((G4LossTableManager::Instance())->SplineFlag()) {
511      size_t n = theLambdaTable->length();
512      for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);}
513    }
514  } else {
515    if (1 < verboseLevel) {
516      G4cout << "Lambda table for " << particleName << " in file <"
517             << filename << "> is not exist"
518             << G4endl;
519    }
520  }
521
522  return yes;
523}
524
525//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
526
527void G4VEmProcess::FindLambdaMax()
528{
529  if(1 < verboseLevel) {
530    G4cout << "### G4VEmProcess::FindLambdaMax: " 
531           << particle->GetParticleName() 
532           << " and process " << GetProcessName() << G4endl; 
533  }
534  size_t n = theLambdaTable->length();
535  G4PhysicsVector* pv = (*theLambdaTable)[0];
536  G4double e, s, emax, smax;
537  theEnergyOfCrossSectionMax = new G4double [n];
538  theCrossSectionMax = new G4double [n];
539  G4bool b;
540
541  for (size_t i=0; i<n; i++) {
542    pv = (*theLambdaTable)[i];
543    emax = DBL_MAX;
544    smax = 0.0;
545    if(pv) {
546      size_t nb = pv->GetVectorLength();
547      emax = pv->GetLowEdgeEnergy(nb);
548      smax = 0.0;
549      for (size_t j=0; j<nb; j++) {
550        e = pv->GetLowEdgeEnergy(j);
551        s = pv->GetValue(e,b);
552        if(s > smax) {
553          smax = s;
554          emax = e;
555        }
556      }
557    }
558    theEnergyOfCrossSectionMax[i] = emax;
559    theCrossSectionMax[i] = smax;
560    if(2 < verboseLevel) {
561      G4cout << "For " << particle->GetParticleName() 
562             << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
563             << " lambda= " << smax << G4endl;
564    }
565  }
566}
567
568//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
569
570G4PhysicsVector* G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
571{
572  G4PhysicsVector* v = 
573    new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
574  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
575  return v;
576}
577
578//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.