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

Last change on this file since 869 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 18.0 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.48 2007/10/29 08:38:58 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-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  selectedModel(0),                   
86  theLambdaTable(0),
87  theEnergyOfCrossSectionMax(0),
88  theCrossSectionMax(0),
89  particle(0),
90  secondaryParticle(0),
91  nLambdaBins(90),
92  lambdaFactor(0.8),
93  currentCouple(0),
94  integral(false),
95  buildLambdaTable(true),
96  applyCuts(false),
97  startFromNull(true),
98  nRegions(0)
99{
100  SetVerboseLevel(1);
101  minKinEnergy = 0.1*keV;
102  maxKinEnergy = 100.0*GeV;
103  theGamma     = G4Gamma::Gamma();
104  theElectron  = G4Electron::Electron();
105  thePositron  = G4Positron::Positron();
106
107  pParticleChange = &fParticleChange;
108  secParticles.reserve(5);
109
110  modelManager = new G4EmModelManager();
111  (G4LossTableManager::Instance())->Register(this);
112}
113
114//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
115
116G4VEmProcess::~G4VEmProcess()
117{
118  if(1 < verboseLevel) 
119    G4cout << "G4VEmProcess destruct " << GetProcessName() 
120           << G4endl;
121  Clear();
122  if(theLambdaTable) {
123    theLambdaTable->clearAndDestroy();
124    delete theLambdaTable;
125  }
126  delete modelManager;
127  (G4LossTableManager::Instance())->DeRegister(this);
128}
129
130//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131
132void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
133{
134  if(!particle) particle = &part;
135  if(1 < verboseLevel) {
136    G4cout << "G4VEmProcess::PreparePhysicsTable() for "
137           << GetProcessName()
138           << " and particle " << part.GetParticleName()
139           << " local particle " << particle->GetParticleName() 
140           << G4endl;
141  }
142
143  if(particle == &part) {
144    Clear();
145    InitialiseProcess(particle);
146    theCuts = modelManager->Initialise(particle,secondaryParticle,2.,verboseLevel);
147    const G4ProductionCutsTable* theCoupleTable=
148          G4ProductionCutsTable::GetProductionCutsTable();
149    theCutsGamma    = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
150    theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
151    theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
152    if(buildLambdaTable)
153      theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
154  }
155}
156
157//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158
159void G4VEmProcess::Clear()
160{
161  if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax;
162  if(theCrossSectionMax) delete [] theCrossSectionMax;
163  theEnergyOfCrossSectionMax = 0;
164  theCrossSectionMax = 0;
165  currentCouple = 0;
166  preStepLambda = 0.0;
167  mfpKinEnergy  = DBL_MAX;
168}
169
170//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171
172void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
173{
174  if(1 < verboseLevel) {
175    G4cout << "G4VEmProcess::BuildPhysicsTable() for "
176           << GetProcessName()
177           << " and particle " << part.GetParticleName()
178           << " buildLambdaTable= " << buildLambdaTable
179           << G4endl;
180  }
181
182  if(buildLambdaTable) {
183    BuildLambdaTable();
184    FindLambdaMax();
185  }
186  if(0 < verboseLevel) PrintInfoDefinition();
187
188  if(1 < verboseLevel) {
189    G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
190           << GetProcessName()
191           << " and particle " << part.GetParticleName()
192           << G4endl;
193  }
194}
195
196//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
197
198void G4VEmProcess::BuildLambdaTable()
199{
200  if(1 < verboseLevel) {
201    G4cout << "G4EmProcess::BuildLambdaTable() for process "
202           << GetProcessName() << " and particle "
203           << particle->GetParticleName()
204           << G4endl;
205  }
206
207  // Access to materials
208  const G4ProductionCutsTable* theCoupleTable=
209        G4ProductionCutsTable::GetProductionCutsTable();
210  size_t numOfCouples = theCoupleTable->GetTableSize();
211  for(size_t i=0; i<numOfCouples; i++) {
212
213    if (theLambdaTable->GetFlag(i)) {
214
215      // create physics vector and fill it
216      const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
217      G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
218      modelManager->FillLambdaVector(aVector, couple, startFromNull);
219      G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
220    }
221  }
222
223  if(1 < verboseLevel) {
224    G4cout << "Lambda table is built for "
225           << particle->GetParticleName()
226           << G4endl;
227    if(2 < verboseLevel) {
228      G4cout << *theLambdaTable << G4endl;
229    }
230  }
231}
232
233//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
234
235G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
236                             const G4Track& track,
237                             G4double   previousStepSize,
238                             G4ForceCondition* condition)
239{
240  // condition is set to "Not Forced"
241  *condition = NotForced;
242  G4double x = DBL_MAX;
243  if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
244  InitialiseStep(track);
245
246  if(preStepKinEnergy < mfpKinEnergy) {
247    if (integral) ComputeIntegralLambda(preStepKinEnergy);
248    else  preStepLambda = GetCurrentLambda(preStepKinEnergy);
249    if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
250  }
251
252  // non-zero cross section
253  if(preStepLambda > DBL_MIN) { 
254    if (theNumberOfInteractionLengthLeft < 0.0) {
255      // beggining of tracking (or just after DoIt of this process)
256      ResetNumberOfInteractionLengthLeft();
257    } else if(currentInteractionLength < DBL_MAX) {
258      // subtract NumberOfInteractionLengthLeft
259      SubtractNumberOfInteractionLengthLeft(previousStepSize);
260      if(theNumberOfInteractionLengthLeft < 0.)
261        theNumberOfInteractionLengthLeft = perMillion;
262    }
263
264    // get mean free path and step limit
265    currentInteractionLength = 1.0/preStepLambda;
266    x = theNumberOfInteractionLengthLeft * currentInteractionLength;
267#ifdef G4VERBOSE
268    if (verboseLevel>2){
269      G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
270      G4cout << "[ " << GetProcessName() << "]" << G4endl; 
271      G4cout << " for " << particle->GetParticleName() 
272             << " in Material  " <<  currentMaterial->GetName()
273             << " Ekin(MeV)= " << preStepKinEnergy/MeV
274             <<G4endl;
275      G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" 
276             << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
277    }
278#endif
279
280    // zero cross section case
281  } else {
282    if(theNumberOfInteractionLengthLeft > DBL_MIN && 
283       currentInteractionLength < DBL_MAX) {
284
285      // subtract NumberOfInteractionLengthLeft
286      SubtractNumberOfInteractionLengthLeft(previousStepSize);
287      if(theNumberOfInteractionLengthLeft < 0.)
288        theNumberOfInteractionLengthLeft = perMillion;
289    }
290    currentInteractionLength = DBL_MAX;
291  }
292  return x;
293}
294
295//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296
297G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
298                                       G4double,
299                                       G4ForceCondition* condition)
300{
301  *condition = NotForced;
302  return G4VEmProcess::MeanFreePath(track);
303}
304
305//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
306
307G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
308                                              const G4Step&)
309{
310  fParticleChange.InitializeForPostStep(track);
311
312  // Do not make anything if particle is stopped, the annihilation then
313  // should be performed by the AtRestDoIt!
314  if (track.GetTrackStatus() == fStopButAlive) return &fParticleChange;
315
316  G4double finalT = track.GetKineticEnergy();
317
318  // Integral approach
319  if (integral) {
320    G4double lx = GetLambda(finalT, currentCouple);
321    if(preStepLambda<lx && 1 < verboseLevel) {
322      G4cout << "WARING: for " << particle->GetParticleName() 
323             << " and " << GetProcessName()
324             << " E(MeV)= " << finalT/MeV
325             << " preLambda= " << preStepLambda << " < " << lx << " (postLambda) "
326             << G4endl; 
327    }
328
329    if(preStepLambda*G4UniformRand() > lx) {
330      ClearNumberOfInteractionLengthLeft();
331      return &fParticleChange;
332    }
333  }
334
335  G4VEmModel* currentModel = SelectModel(finalT);
336
337  /* 
338  if(0 < verboseLevel) {
339    G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
340           << finalT/MeV
341           << " MeV; model= (" << currentModel->LowEnergyLimit()
342           << ", " <<  currentModel->HighEnergyLimit() << ")"
343           << G4endl;
344  }
345  */
346
347 
348  // sample secondaries
349  secParticles.clear();
350  currentModel->SampleSecondaries(&secParticles, 
351                                  currentCouple, 
352                                  track.GetDynamicParticle(),
353                                  (*theCuts)[currentMaterialIndex]);
354
355  // save secondaries
356  G4int num = secParticles.size();
357  if(num > 0) {
358
359    fParticleChange.SetNumberOfSecondaries(num);
360    G4double edep = fParticleChange.GetLocalEnergyDeposit();
361     
362    for (G4int i=0; i<num; i++) {
363      G4DynamicParticle* dp = secParticles[i];
364      const G4ParticleDefinition* p = dp->GetDefinition();
365      G4double e = dp->GetKineticEnergy();
366      G4bool good = true;
367      if(applyCuts) {
368        if (p == theGamma) {
369          if (e < (*theCutsGamma)[currentMaterialIndex]) good = false;
370
371        } else if (p == theElectron) {
372          if (e < (*theCutsElectron)[currentMaterialIndex]) good = false;
373
374        } else if (p == thePositron) {
375          if (e < (*theCutsPositron)[currentMaterialIndex]) {
376            good = false;
377            e += 2.0*electron_mass_c2;
378          }
379        }
380        if(!good) {
381          delete dp;
382          edep += e;
383        }
384      }
385      if (good) fParticleChange.AddSecondary(dp);
386    } 
387    fParticleChange.ProposeLocalEnergyDeposit(edep);
388  }
389
390  ClearNumberOfInteractionLengthLeft();
391  return &fParticleChange;
392}
393
394//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
395
396void G4VEmProcess::PrintInfoDefinition()
397{
398  if(verboseLevel > 0) {
399    G4cout << G4endl << GetProcessName() << ": " ;
400    PrintInfo();
401    if(integral) {
402      G4cout << "      Integral mode is used  "<< G4endl;
403    }
404  }
405
406  if (!buildLambdaTable)  return;
407 
408  if(verboseLevel > 0) {
409    G4cout << "      tables are built for  "
410           << particle->GetParticleName()
411           << G4endl
412           << "      Lambda tables from "
413           << G4BestUnit(minKinEnergy,"Energy") 
414           << " to "
415           << G4BestUnit(maxKinEnergy,"Energy")
416           << " in " << nLambdaBins << " bins."
417           << G4endl;
418  }
419
420  if(verboseLevel > 1) {
421    G4cout << "Tables are built for " << particle->GetParticleName()
422           << G4endl;
423
424  if(verboseLevel > 2) {
425    G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
426    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
427    }
428  }
429}
430
431//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
432
433G4double G4VEmProcess::MicroscopicCrossSection(G4double kineticEnergy,
434                                         const G4MaterialCutsCouple* couple)
435{
436  // Cross section per atom is calculated
437  DefineMaterial(couple);
438  G4double cross = 0.0;
439  G4bool b;
440  if(theLambdaTable) {
441    cross = (((*theLambdaTable)[currentMaterialIndex])->
442                           GetValue(kineticEnergy, b));
443
444    cross /= currentMaterial->GetTotNbOfAtomsPerVolume();
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  } else {
511    if (1 < verboseLevel) {
512      G4cout << "Lambda table for " << particleName << " in file <"
513             << filename << "> is not exist"
514             << G4endl;
515    }
516  }
517
518  return yes;
519}
520
521//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
522
523void G4VEmProcess::FindLambdaMax()
524{
525  if(1 < verboseLevel) {
526    G4cout << "### G4VEmProcess::FindLambdaMax: " 
527           << particle->GetParticleName() 
528           << " and process " << GetProcessName() << G4endl; 
529  }
530  size_t n = theLambdaTable->length();
531  G4PhysicsVector* pv = (*theLambdaTable)[0];
532  G4double e, s, emax, smax;
533  theEnergyOfCrossSectionMax = new G4double [n];
534  theCrossSectionMax = new G4double [n];
535  G4bool b;
536
537  for (size_t i=0; i<n; i++) {
538    pv = (*theLambdaTable)[i];
539    emax = DBL_MAX;
540    smax = 0.0;
541    if(pv) {
542      size_t nb = pv->GetVectorLength();
543      emax = pv->GetLowEdgeEnergy(nb);
544      smax = 0.0;
545      for (size_t j=0; j<nb; j++) {
546        e = pv->GetLowEdgeEnergy(j);
547        s = pv->GetValue(e,b);
548        if(s > smax) {
549          smax = s;
550          emax = e;
551        }
552      }
553    }
554    theEnergyOfCrossSectionMax[i] = emax;
555    theCrossSectionMax[i] = smax;
556    if(2 < verboseLevel) {
557      G4cout << "For " << particle->GetParticleName() 
558             << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
559             << " lambda= " << smax << G4endl;
560    }
561  }
562}
563
564//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
565
566G4PhysicsVector* G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
567{
568  G4PhysicsVector* v = 
569    new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
570  return v;
571}
572
573//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.