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

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

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