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

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

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