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

Last change on this file since 1275 was 1228, checked in by garnier, 16 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.