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

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

update geant4.9.3 tag
File size: 22.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: G4VEmProcess.cc,v 1.79 2009/11/10 20:30:55 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03 $
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// 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
55//
56// Class Description:
57//
58
59// -------------------------------------------------------------------
60//
61//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63
64#include "G4VEmProcess.hh"
65#include "G4LossTableManager.hh"
66#include "G4Step.hh"
67#include "G4ParticleDefinition.hh"
68#include "G4VEmModel.hh"
69#include "G4DataVector.hh"
70#include "G4PhysicsTable.hh"
71#include "G4PhysicsVector.hh"
72#include "G4PhysicsLogVector.hh"
73#include "G4VParticleChange.hh"
74#include "G4ProductionCutsTable.hh"
75#include "G4Region.hh"
76#include "G4RegionStore.hh"
77#include "G4Gamma.hh"
78#include "G4Electron.hh"
79#include "G4Positron.hh"
80#include "G4PhysicsTableHelper.hh"
81#include "G4EmConfigurator.hh"
82
83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84
85G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
86  G4VDiscreteProcess(name, type),
87  secondaryParticle(0),
88  buildLambdaTable(true),
89  theLambdaTable(0),
90  theEnergyOfCrossSectionMax(0),
91  theCrossSectionMax(0),
92  integral(false),
93  applyCuts(false),
94  startFromNull(true),
95  useDeexcitation(false),
96  nDERegions(0),
97  idxDERegions(0),
98  currentModel(0),
99  particle(0),
100  currentCouple(0)
101{
102  SetVerboseLevel(1);
103
104  // Size of tables assuming spline
105  minKinEnergy = 0.1*keV;
106  maxKinEnergy = 10.0*TeV;
107  nLambdaBins  = 77;
108
109  // default lambda factor
110  lambdaFactor  = 0.8;
111
112  // default limit on polar angle
113  polarAngleLimit = 0.0;
114
115  // particle types
116  theGamma     = G4Gamma::Gamma();
117  theElectron  = G4Electron::Electron();
118  thePositron  = G4Positron::Positron();
119
120  pParticleChange = &fParticleChange;
121  secParticles.reserve(5);
122
123  modelManager = new G4EmModelManager();
124  (G4LossTableManager::Instance())->Register(this);
125}
126
127//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
128
129G4VEmProcess::~G4VEmProcess()
130{
131  if(1 < verboseLevel) 
132    G4cout << "G4VEmProcess destruct " << GetProcessName() 
133           << G4endl;
134  Clear();
135  if(theLambdaTable) {
136    theLambdaTable->clearAndDestroy();
137    delete theLambdaTable;
138  }
139  delete modelManager;
140  (G4LossTableManager::Instance())->DeRegister(this);
141}
142
143//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144
145void G4VEmProcess::Clear()
146{
147  delete [] theEnergyOfCrossSectionMax;
148  delete [] theCrossSectionMax;
149  delete [] idxDERegions;
150  theEnergyOfCrossSectionMax = 0;
151  theCrossSectionMax = 0;
152  idxDERegions = 0;
153  currentCouple = 0;
154  preStepLambda = 0.0;
155  mfpKinEnergy  = DBL_MAX;
156  deRegions.clear();
157  nDERegions = 0;
158}
159
160//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161
162void G4VEmProcess::AddEmModel(G4int order, G4VEmModel* p, 
163                              const G4Region* region)
164{
165  G4VEmFluctuationModel* fm = 0;
166  modelManager->AddEmModel(order, p, fm, region);
167  if(p) p->SetParticleChange(pParticleChange);
168}
169
170//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
171
172void G4VEmProcess::SetModel(G4VEmModel* p, G4int index)
173{
174  G4int n = emModels.size();
175  if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
176  emModels[index] = p;
177}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180
181G4VEmModel* G4VEmProcess::Model(G4int index)
182{
183  G4VEmModel* p = 0;
184  if(index >= 0 && index <  G4int(emModels.size())) { p = emModels[index]; }
185  return p;
186}
187
188//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189
190void G4VEmProcess::UpdateEmModel(const G4String& nam, 
191                                 G4double emin, G4double emax)
192{
193  modelManager->UpdateEmModel(nam, emin, emax);
194}
195
196//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
197
198G4VEmModel* G4VEmProcess::GetModelByIndex(G4int idx, G4bool ver)
199{
200  return modelManager->GetModel(idx, ver);
201}
202
203//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
204
205void G4VEmProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
206{
207  if(!particle) particle = &part;
208  if(1 < verboseLevel) {
209    G4cout << "G4VEmProcess::PreparePhysicsTable() for "
210           << GetProcessName()
211           << " and particle " << part.GetParticleName()
212           << " local particle " << particle->GetParticleName() 
213           << G4endl;
214  }
215
216  (G4LossTableManager::Instance())->EmConfigurator()->AddModels();
217
218  if(particle == &part) {
219    Clear();
220    InitialiseProcess(particle);
221
222    // initialisation of models
223    G4int nmod = modelManager->NumberOfModels();
224    for(G4int i=0; i<nmod; ++i) {
225      G4VEmModel* mod = modelManager->GetModel(i);
226      mod->SetPolarAngleLimit(polarAngleLimit);
227      if(mod->HighEnergyLimit() > maxKinEnergy) {
228        mod->SetHighEnergyLimit(maxKinEnergy);
229      }
230    }
231
232    theCuts = modelManager->Initialise(particle,secondaryParticle,2.,verboseLevel);
233    const G4ProductionCutsTable* theCoupleTable=
234          G4ProductionCutsTable::GetProductionCutsTable();
235    theCutsGamma    = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
236    theCutsElectron = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
237    theCutsPositron = theCoupleTable->GetEnergyCutsVector(idxG4PositronCut);
238
239    // prepare tables
240    if(buildLambdaTable){
241      theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
242    }
243  }
244  // Sub Cutoff and Deexcitation
245  if (nDERegions>0) {
246
247    const G4ProductionCutsTable* theCoupleTable=
248          G4ProductionCutsTable::GetProductionCutsTable();
249    size_t numOfCouples = theCoupleTable->GetTableSize();
250
251    idxDERegions = new G4bool[numOfCouples];
252
253    for (size_t j=0; j<numOfCouples; ++j) {
254
255      const G4MaterialCutsCouple* couple =
256        theCoupleTable->GetMaterialCutsCouple(j);
257      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
258      G4bool reg = false;
259      for(G4int i=0; i<nDERegions; ++i) {
260        if(deRegions[i]) {
261          if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
262        }
263      }
264      idxDERegions[j] = reg;
265    }
266  }
267  if (1 < verboseLevel && nDERegions>0) {
268    G4cout << " Deexcitation is activated for regions: " << G4endl;
269    for (G4int i=0; i<nDERegions; ++i) {
270      const G4Region* r = deRegions[i];
271      G4cout << "           " << r->GetName() << G4endl;
272    }
273  }
274}
275
276//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277
278void G4VEmProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
279{
280  G4String partname = part.GetParticleName();
281  if(1 < verboseLevel) {
282    G4cout << "G4VEmProcess::BuildPhysicsTable() for "
283           << GetProcessName()
284           << " and particle " << partname
285           << " buildLambdaTable= " << buildLambdaTable
286           << G4endl;
287  }
288
289  if(buildLambdaTable) {
290    BuildLambdaTable();
291    FindLambdaMax();
292  }
293
294  // reduce printout for nuclear stopping
295  G4bool gproc = true;
296  if(GetProcessName() == "nuclearStopping" && 
297     partname != "GenericIon" && partname != "alpha") { gproc = false; } 
298
299  if(gproc && 0 < verboseLevel) { PrintInfoDefinition(); }
300
301  if(1 < verboseLevel) {
302    G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
303           << GetProcessName()
304           << " and particle " << partname
305           << G4endl;
306  }
307}
308
309//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
310
311void G4VEmProcess::BuildLambdaTable()
312{
313  if(1 < verboseLevel) {
314    G4cout << "G4EmProcess::BuildLambdaTable() for process "
315           << GetProcessName() << " and particle "
316           << particle->GetParticleName()
317           << G4endl;
318  }
319
320  // Access to materials
321  const G4ProductionCutsTable* theCoupleTable=
322        G4ProductionCutsTable::GetProductionCutsTable();
323  size_t numOfCouples = theCoupleTable->GetTableSize();
324
325  G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
326
327  G4PhysicsLogVector* aVector = 0;
328  G4PhysicsLogVector* bVector = 0;
329
330  for(size_t i=0; i<numOfCouples; ++i) {
331
332    if (theLambdaTable->GetFlag(i)) {
333
334      // create physics vector and fill it
335      const G4MaterialCutsCouple* couple = 
336        theCoupleTable->GetMaterialCutsCouple(i);
337      if(!bVector) {
338        aVector = 
339          static_cast<G4PhysicsLogVector*>(LambdaPhysicsVector(couple));
340        bVector = aVector;
341      } else {
342        aVector = new G4PhysicsLogVector(*bVector);
343      }
344        //      G4PhysicsVector* aVector = LambdaPhysicsVector(couple);
345      aVector->SetSpline(splineFlag);
346      modelManager->FillLambdaVector(aVector, couple, startFromNull);
347      if(splineFlag) aVector->FillSecondDerivatives();
348      G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
349    }
350  }
351
352  if(1 < verboseLevel) {
353    G4cout << "Lambda table is built for "
354           << particle->GetParticleName()
355           << G4endl;
356    if(2 < verboseLevel) {
357      G4cout << *theLambdaTable << G4endl;
358    }
359  }
360}
361
362//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
363
364void G4VEmProcess::PrintInfoDefinition()
365{
366  if(verboseLevel > 0) {
367    G4cout << G4endl << GetProcessName() << ":   for  "
368           << particle->GetParticleName();
369    if(integral) G4cout << ", integral: 1 ";
370    if(applyCuts) G4cout << ", applyCuts: 1 ";
371    G4cout << "    SubType= " << GetProcessSubType() << G4endl;
372    if(buildLambdaTable) {
373      G4cout << "      Lambda tables from "
374             << G4BestUnit(minKinEnergy,"Energy") 
375             << " to "
376             << G4BestUnit(maxKinEnergy,"Energy")
377             << " in " << nLambdaBins << " bins, spline: " 
378             << (G4LossTableManager::Instance())->SplineFlag()
379             << G4endl;
380    }
381    PrintInfo();
382    modelManager->DumpModelList(verboseLevel);
383  }
384
385  if(verboseLevel > 2 && buildLambdaTable) {
386    G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
387    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
388  }
389}
390
391//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
392
393G4double G4VEmProcess::PostStepGetPhysicalInteractionLength(
394                             const G4Track& track,
395                             G4double   previousStepSize,
396                             G4ForceCondition* condition)
397{
398  // condition is set to "Not Forced"
399  *condition = NotForced;
400  G4double x = DBL_MAX;
401  if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
402  InitialiseStep(track);
403  if(!currentModel->IsActive(preStepKinEnergy)) return x;
404
405  if(preStepKinEnergy < mfpKinEnergy) {
406    if (integral) ComputeIntegralLambda(preStepKinEnergy);
407    else  preStepLambda = GetCurrentLambda(preStepKinEnergy);
408    if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
409  }
410
411  // non-zero cross section
412  if(preStepLambda > DBL_MIN) { 
413    if (theNumberOfInteractionLengthLeft < 0.0) {
414      // beggining of tracking (or just after DoIt of this process)
415      ResetNumberOfInteractionLengthLeft();
416    } else if(currentInteractionLength < DBL_MAX) {
417      // subtract NumberOfInteractionLengthLeft
418      SubtractNumberOfInteractionLengthLeft(previousStepSize);
419      if(theNumberOfInteractionLengthLeft < 0.)
420        theNumberOfInteractionLengthLeft = perMillion;
421    }
422
423    // get mean free path and step limit
424    currentInteractionLength = 1.0/preStepLambda;
425    x = theNumberOfInteractionLengthLeft * currentInteractionLength;
426#ifdef G4VERBOSE
427    if (verboseLevel>2){
428      G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
429      G4cout << "[ " << GetProcessName() << "]" << G4endl; 
430      G4cout << " for " << particle->GetParticleName() 
431             << " in Material  " <<  currentMaterial->GetName()
432             << " Ekin(MeV)= " << preStepKinEnergy/MeV
433             <<G4endl;
434      G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" 
435             << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
436    }
437#endif
438
439    // zero cross section case
440  } else {
441    if(theNumberOfInteractionLengthLeft > DBL_MIN && 
442       currentInteractionLength < DBL_MAX) {
443
444      // subtract NumberOfInteractionLengthLeft
445      SubtractNumberOfInteractionLengthLeft(previousStepSize);
446      if(theNumberOfInteractionLengthLeft < 0.)
447        theNumberOfInteractionLengthLeft = perMillion;
448    }
449    currentInteractionLength = DBL_MAX;
450  }
451  return x;
452}
453
454//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
455
456G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
457                                              const G4Step&)
458{
459  fParticleChange.InitializeForPostStep(track);
460
461  // Do not make anything if particle is stopped, the annihilation then
462  // should be performed by the AtRestDoIt!
463  if (track.GetTrackStatus() == fStopButAlive) return &fParticleChange;
464
465  G4double finalT = track.GetKineticEnergy();
466
467  // Integral approach
468  if (integral) {
469    G4double lx = GetLambda(finalT, currentCouple);
470    if(preStepLambda<lx && 1 < verboseLevel) {
471      G4cout << "WARING: for " << particle->GetParticleName() 
472             << " and " << GetProcessName()
473             << " E(MeV)= " << finalT/MeV
474             << " preLambda= " << preStepLambda << " < " << lx << " (postLambda) "
475             << G4endl; 
476    }
477
478    if(preStepLambda*G4UniformRand() > lx) {
479      ClearNumberOfInteractionLengthLeft();
480      return &fParticleChange;
481    }
482  }
483
484  SelectModel(finalT, currentCoupleIndex);
485  if(!currentModel->IsActive(finalT)) return &fParticleChange;
486  if(useDeexcitation) {
487    currentModel->SetDeexcitationFlag(idxDERegions[currentCoupleIndex]);
488  }
489  /* 
490  if(0 < verboseLevel) {
491    G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
492           << finalT/MeV
493           << " MeV; model= (" << currentModel->LowEnergyLimit()
494           << ", " <<  currentModel->HighEnergyLimit() << ")"
495           << G4endl;
496  }
497  */
498
499 
500  // sample secondaries
501  secParticles.clear();
502  currentModel->SampleSecondaries(&secParticles, 
503                                  currentCouple, 
504                                  track.GetDynamicParticle(),
505                                  (*theCuts)[currentCoupleIndex]);
506
507  // save secondaries
508  G4int num = secParticles.size();
509  if(num > 0) {
510
511    fParticleChange.SetNumberOfSecondaries(num);
512    G4double edep = fParticleChange.GetLocalEnergyDeposit();
513     
514    for (G4int i=0; i<num; ++i) {
515      G4DynamicParticle* dp = secParticles[i];
516      const G4ParticleDefinition* p = dp->GetDefinition();
517      G4double e = dp->GetKineticEnergy();
518      G4bool good = true;
519      if(applyCuts) {
520        if (p == theGamma) {
521          if (e < (*theCutsGamma)[currentCoupleIndex]) good = false;
522
523        } else if (p == theElectron) {
524          if (e < (*theCutsElectron)[currentCoupleIndex]) good = false;
525
526        } else if (p == thePositron) {
527          if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
528              e < (*theCutsPositron)[currentCoupleIndex]) {
529            good = false;
530            e += 2.0*electron_mass_c2;
531          }
532        }
533        if(!good) {
534          delete dp;
535          edep += e;
536        }
537      }
538      if (good) fParticleChange.AddSecondary(dp);
539    } 
540    fParticleChange.ProposeLocalEnergyDeposit(edep);
541  }
542
543  ClearNumberOfInteractionLengthLeft();
544  return &fParticleChange;
545}
546
547//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
548
549G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
550                                       const G4String& directory,
551                                             G4bool ascii)
552{
553  G4bool yes = true;
554
555  if ( theLambdaTable && part == particle) {
556    const G4String name = 
557      GetPhysicsTableFileName(part,directory,"Lambda",ascii);
558    yes = theLambdaTable->StorePhysicsTable(name,ascii);
559
560    if ( yes ) {
561      G4cout << "Physics tables are stored for " << particle->GetParticleName()
562             << " and process " << GetProcessName()
563             << " in the directory <" << directory
564             << "> " << G4endl;
565    } else {
566      G4cout << "Fail to store Physics Tables for " 
567             << particle->GetParticleName()
568             << " and process " << GetProcessName()
569             << " in the directory <" << directory
570             << "> " << G4endl;
571    }
572  }
573  return yes;
574}
575
576//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
577
578G4bool G4VEmProcess::RetrievePhysicsTable(const G4ParticleDefinition* part,
579                                          const G4String& directory,
580                                          G4bool ascii)
581{
582  if(1 < verboseLevel) {
583    G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
584           << part->GetParticleName() << " and process "
585           << GetProcessName() << G4endl;
586  }
587  G4bool yes = true;
588
589  if(!buildLambdaTable || particle != part) return yes;
590
591  const G4String particleName = part->GetParticleName();
592  G4String filename;
593
594  filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
595  yes = G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,
596                                                   filename,ascii);
597  if ( yes ) {
598    if (0 < verboseLevel) {
599      G4cout << "Lambda table for " << particleName
600             << " is Retrieved from <"
601             << filename << ">"
602             << G4endl;
603    }
604    if((G4LossTableManager::Instance())->SplineFlag()) {
605      size_t n = theLambdaTable->length();
606      for(size_t i=0; i<n; ++i) {
607        if((* theLambdaTable)[i]) {
608          (* theLambdaTable)[i]->SetSpline(true);
609        }
610      }
611    }
612  } else {
613    if (1 < verboseLevel) {
614      G4cout << "Lambda table for " << particleName << " in file <"
615             << filename << "> is not exist"
616             << G4endl;
617    }
618  }
619
620  return yes;
621}
622
623//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
624
625void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
626{
627  G4RegionStore* regionStore = G4RegionStore::GetInstance();
628  const G4Region* reg = r;
629  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
630
631  // the region is in the list
632  if (nDERegions) {
633    for (G4int i=0; i<nDERegions; ++i) {
634      if (reg == deRegions[i]) {
635        if(!val) deRegions[i] = 0;
636        return;
637      }
638    }
639  }
640
641  // new region
642  if(val) {
643    useDeexcitation = true;
644    deRegions.push_back(reg);
645    nDERegions++;
646  } else {
647    useDeexcitation = false;
648  }
649}
650
651//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
652
653G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
654                                             const G4MaterialCutsCouple* couple)
655{
656  // Cross section per atom is calculated
657  DefineMaterial(couple);
658  G4double cross = 0.0;
659  if(theLambdaTable) {
660    cross = (((*theLambdaTable)[currentCoupleIndex])->Value(kineticEnergy));
661  } else {
662    SelectModel(kineticEnergy, currentCoupleIndex);
663    cross = currentModel->CrossSectionPerVolume(currentMaterial,
664                                                particle,kineticEnergy);
665  }
666
667  return cross;
668}
669
670//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
671
672G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
673                                       G4double,
674                                       G4ForceCondition* condition)
675{
676  *condition = NotForced;
677  return G4VEmProcess::MeanFreePath(track);
678}
679
680//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
681
682void G4VEmProcess::FindLambdaMax()
683{
684  if(1 < verboseLevel) {
685    G4cout << "### G4VEmProcess::FindLambdaMax: " 
686           << particle->GetParticleName() 
687           << " and process " << GetProcessName() << G4endl; 
688  }
689  size_t n = theLambdaTable->length();
690  G4PhysicsVector* pv = (*theLambdaTable)[0];
691  G4double e, s, emax, smax;
692  theEnergyOfCrossSectionMax = new G4double [n];
693  theCrossSectionMax = new G4double [n];
694
695  for (size_t i=0; i<n; ++i) {
696    pv = (*theLambdaTable)[i];
697    emax = DBL_MAX;
698    smax = 0.0;
699    if(pv) {
700      size_t nb = pv->GetVectorLength();
701      emax = DBL_MAX;
702      smax = 0.0;
703      if(nb > 0) {
704        for (size_t j=0; j<nb; ++j) {
705          e = pv->Energy(j);
706          s = (*pv)(j);
707          if(s > smax) {
708            smax = s;
709            emax = e;
710          }
711        }
712      }
713    }
714    theEnergyOfCrossSectionMax[i] = emax;
715    theCrossSectionMax[i] = smax;
716    if(2 < verboseLevel) {
717      G4cout << "For " << particle->GetParticleName() 
718             << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
719             << " lambda= " << smax << G4endl;
720    }
721  }
722}
723
724//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
725
726G4PhysicsVector* G4VEmProcess::LambdaPhysicsVector(const G4MaterialCutsCouple*)
727{
728  G4PhysicsVector* v = 
729    new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
730  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
731  return v;
732}
733
734//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Note: See TracBrowser for help on using the repository browser.