source: trunk/source/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc@ 1153

Last change on this file since 1153 was 1055, checked in by garnier, 17 years ago

maj sur la beta de geant 4.9.3
File size: 56.5 KB
RevLine 
[819]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//
[1055]26// $Id: G4VEnergyLossProcess.cc,v 1.149 2009/04/17 10:35:32 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
[819]28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4VEnergyLossProcess
35//
36// Author: Vladimir Ivanchenko
37//
38// Creation date: 03.01.2002
39//
40// Modifications:
41//
42// 13-11-02 Minor fix - use normalised direction (V.Ivanchenko)
43// 04-12-02 Minor change in PostStepDoIt (V.Ivanchenko)
44// 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
45// 26-12-02 Secondary production moved to derived classes (V.Ivanchenko)
46// 04-01-03 Fix problem of very small steps for ions (V.Ivanchenko)
47// 20-01-03 Migrade to cut per region (V.Ivanchenko)
48// 24-01-03 Temporarily close a control on usage of couples (V.Ivanchenko)
49// 24-01-03 Make models region aware (V.Ivanchenko)
50// 05-02-03 Fix compilation warnings (V.Ivanchenko)
51// 06-02-03 Add control on tmax in PostStepDoIt (V.Ivanchenko)
52// 13-02-03 SubCutoffProcessors defined for regions (V.Ivanchenko)
53// 15-02-03 Lambda table can be scaled (V.Ivanchenko)
54// 17-02-03 Fix problem of store/restore tables (V.Ivanchenko)
55// 18-02-03 Add control on CutCouple usage (V.Ivanchenko)
56// 26-02-03 Simplify control on GenericIons (V.Ivanchenko)
57// 06-03-03 Control on GenericIons using SubType + update verbose (V.Ivanchenko)
58// 10-03-03 Add Ion registration (V.Ivanchenko)
59// 22-03-03 Add Initialisation of cash (V.Ivanchenko)
60// 26-03-03 Remove finalRange modification (V.Ivanchenko)
61// 09-04-03 Fix problem of negative range limit for non integral (V.Ivanchenko)
62// 26-04-03 Fix retrieve tables (V.Ivanchenko)
63// 06-05-03 Set defalt finalRange = 1 mm (V.Ivanchenko)
64// 12-05-03 Update range calculations + lowKinEnergy (V.Ivanchenko)
65// 13-05-03 Add calculation of precise range (V.Ivanchenko)
66// 23-05-03 Remove tracking cuts (V.Ivanchenko)
67// 03-06-03 Fix initialisation problem for STD ionisation (V.Ivanchenko)
68// 21-07-03 Add UpdateEmModel method (V.Ivanchenko)
69// 03-11-03 Fix initialisation problem in RetrievePhysicsTable (V.Ivanchenko)
70// 04-11-03 Add checks in RetrievePhysicsTable (V.Ivanchenko)
71// 12-11-03 G4EnergyLossSTD -> G4EnergyLossProcess (V.Ivanchenko)
72// 21-01-04 Migrade to G4ParticleChangeForLoss (V.Ivanchenko)
73// 27-02-04 Fix problem of loss in low presure gases, cleanup precise range
74// calculation, use functions ForLoss in AlongStepDoIt (V.Ivanchenko)
75// 10-03-04 Fix a problem of Precise Range table (V.Ivanchenko)
76// 19-03-04 Fix a problem energy below lowestKinEnergy (V.Ivanchenko)
77// 31-03-04 Fix a problem of retrieve tables (V.Ivanchenko)
78// 21-07-04 Check weather AtRest are active or not (V.Ivanchenko)
79// 03-08-04 Add pointer of DEDX table to all processes (V.Ivanchenko)
80// 06-08-04 Clear up names of member functions (V.Ivanchenko)
81// 06-08-04 Clear up names of member functions (V.Ivanchenko)
82// 27-08-04 Add NeedBuildTables method (V.Ivanchneko)
83// 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko)
84// 11-03-05 Shift verbose level by 1 (V.Ivantchenko)
85// 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
86// 11-04-05 Use MaxSecondaryEnergy from a model (V.Ivanchenko)
87// 25-07-05 Add extra protection PostStep for non-integral mode (V.Ivanchenko)
88// 12-08-05 Integral=false; SetStepFunction(0.2, 0.1*mm) (mma)
89// 18-08-05 Return back both AlongStep and PostStep from 7.0 (V.Ivanchenko)
90// 02-09-05 Default StepFunction 0.2 1 mm + integral (V.Ivanchenko)
91// 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
92// 05-10-05 protection against 0 energy loss added (L.Urban)
93// 17-10-05 protection above has been removed (L.Urban)
94// 06-01-06 reset currentCouple when StepFunction is changed (V.Ivanchenko)
95// 10-01-06 PreciseRange -> CSDARange (V.Ivantchenko)
96// 18-01-06 Clean up subcutoff including recalculation of presafety (VI)
97// 20-01-06 Introduce G4EmTableType and reducing number of methods (VI)
98// 22-03-06 Add control on warning printout AlongStep (VI)
99// 23-03-06 Use isIonisation flag (V.Ivanchenko)
100// 07-06-06 Do not reflect AlongStep in subcutoff regime (V.Ivanchenko)
101// 14-01-07 add SetEmModel(index) and SetFluctModel() (mma)
102// 16-01-07 add IonisationTable and IonisationSubTable (V.Ivanchenko)
103// 16-02-07 set linLossLimit=1.e-6 (V.Ivanchenko)
104// 13-03-07 use SafetyHelper instead of navigator (V.Ivanchenko)
105// 10-04-07 use unique SafetyHelper (V.Ivanchenko)
106// 12-04-07 Add verbosity at destruction (V.Ivanchenko)
107// 25-04-07 move initialisation of safety helper to BuildPhysicsTable (VI)
108// 27-10-07 Virtual functions moved to source (V.Ivanchenko)
109//
110// Class Description:
111//
112// It is the unified energy loss process it calculates the continuous
113// energy loss for charged particles using a set of Energy Loss
114// models valid for different energy regions. There are a possibility
115// to create and access to dE/dx and range tables, or to calculate
116// that information on fly.
117// -------------------------------------------------------------------
118//
119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
121
122#include "G4VEnergyLossProcess.hh"
123#include "G4LossTableManager.hh"
124#include "G4Step.hh"
125#include "G4ParticleDefinition.hh"
126#include "G4VEmModel.hh"
127#include "G4VEmFluctuationModel.hh"
128#include "G4DataVector.hh"
129#include "G4PhysicsLogVector.hh"
130#include "G4VParticleChange.hh"
131#include "G4Gamma.hh"
132#include "G4Electron.hh"
133#include "G4Positron.hh"
134#include "G4Proton.hh"
135#include "G4ProcessManager.hh"
136#include "G4UnitsTable.hh"
137#include "G4GenericIon.hh"
138#include "G4ProductionCutsTable.hh"
139#include "G4Region.hh"
140#include "G4RegionStore.hh"
141#include "G4PhysicsTableHelper.hh"
142#include "G4SafetyHelper.hh"
143#include "G4TransportationManager.hh"
[1055]144#include "G4EmConfigurator.hh"
[819]145
146//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
148G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name,
[961]149 G4ProcessType type):
150 G4VContinuousDiscreteProcess(name, type),
151 secondaryParticle(0),
[819]152 nSCoffRegions(0),
[1055]153 nDERegions(0),
[819]154 idxSCoffRegions(0),
[1055]155 idxDERegions(0),
[819]156 nProcesses(0),
157 theDEDXTable(0),
158 theDEDXSubTable(0),
159 theDEDXunRestrictedTable(0),
160 theIonisationTable(0),
161 theIonisationSubTable(0),
162 theRangeTableForLoss(0),
163 theCSDARangeTable(0),
164 theSecondaryRangeTable(0),
165 theInverseRangeTable(0),
166 theLambdaTable(0),
167 theSubLambdaTable(0),
168 theDEDXAtMaxEnergy(0),
169 theRangeAtMaxEnergy(0),
170 theEnergyOfCrossSectionMax(0),
171 theCrossSectionMax(0),
172 baseParticle(0),
173 minSubRange(0.1),
174 lossFluctuationFlag(true),
175 rndmStepFlag(false),
176 tablesAreBuilt(false),
177 integral(true),
[961]178 isIon(false),
[819]179 isIonisation(true),
[961]180 useSubCutoff(false),
[1055]181 useDeexcitation(false),
[961]182 particle(0),
183 currentCouple(0),
184 nWarnings(0),
185 mfpKinEnergy(0.0)
[819]186{
187 SetVerboseLevel(1);
188
[961]189 // low energy limit
[819]190 lowestKinEnergy = 1.*eV;
[961]191
192 // Size of tables assuming spline
[819]193 minKinEnergy = 0.1*keV;
[1055]194 maxKinEnergy = 10.0*TeV;
195 nBins = 77;
[819]196 maxKinEnergyCSDA = 1.0*GeV;
[961]197 nBinsCSDA = 35;
[819]198
[961]199 // default linear loss limit for spline
200 linLossLimit = 0.01;
201
[819]202 // default dRoverRange and finalRange
203 SetStepFunction(0.2, 1.0*mm);
204
[961]205 // default lambda factor
206 lambdaFactor = 0.8;
[819]207
[961]208 // particle types
209 theElectron = G4Electron::Electron();
210 thePositron = G4Positron::Positron();
211 theGenericIon = 0;
212
[819]213 // run time objects
214 pParticleChange = &fParticleChange;
215 modelManager = new G4EmModelManager();
216 safetyHelper = G4TransportationManager::GetTransportationManager()
217 ->GetSafetyHelper();
218 aGPILSelection = CandidateForSelection;
219
220 // initialise model
221 (G4LossTableManager::Instance())->Register(this);
222 fluctModel = 0;
223
224 scTracks.reserve(5);
225 secParticles.reserve(5);
226
227 // Data for stragling of ranges from ICRU'37 report
228 const G4int nrbins = 7;
[961]229 vstrag = new G4PhysicsLogVector(keV, GeV, nrbins-1);
230 vstrag->SetSpline(true);
[819]231 G4double s[nrbins] = {-0.2, -0.85, -1.3, -1.578, -1.76, -1.85, -1.9};
232 for(G4int i=0; i<nrbins; i++) {vstrag->PutValue(i, s[i]);}
233}
234
235//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236
237G4VEnergyLossProcess::~G4VEnergyLossProcess()
238{
239 if(1 < verboseLevel)
240 G4cout << "G4VEnergyLossProcess destruct " << GetProcessName()
241 << G4endl;
242 delete vstrag;
[1055]243 Clean();
[819]244
245 if ( !baseParticle ) {
246 if(theDEDXTable && theRangeTableForLoss) {
247 if(theIonisationTable == theDEDXTable) theIonisationTable = 0;
248 theDEDXTable->clearAndDestroy();
249 delete theDEDXTable;
250 if(theDEDXSubTable) {
251 if(theIonisationSubTable == theDEDXSubTable)
252 theIonisationSubTable = 0;
253 theDEDXSubTable->clearAndDestroy();
254 delete theDEDXSubTable;
255 }
256 }
257 if(theIonisationTable) {
258 theIonisationTable->clearAndDestroy();
259 delete theIonisationTable;
260 }
261 if(theIonisationSubTable) {
262 theIonisationSubTable->clearAndDestroy();
263 delete theIonisationSubTable;
264 }
265 if(theDEDXunRestrictedTable && theCSDARangeTable) {
266 theDEDXunRestrictedTable->clearAndDestroy();
267 delete theDEDXunRestrictedTable;
268 }
269 if(theCSDARangeTable) {
270 theCSDARangeTable->clearAndDestroy();
271 delete theCSDARangeTable;
272 }
273 if(theRangeTableForLoss) {
274 theRangeTableForLoss->clearAndDestroy();
275 delete theRangeTableForLoss;
276 }
277 if(theInverseRangeTable) {
278 theInverseRangeTable->clearAndDestroy();
279 delete theInverseRangeTable;
280 }
281 if(theLambdaTable) {
282 theLambdaTable->clearAndDestroy();
283 delete theLambdaTable;
284 }
285 if(theSubLambdaTable) {
286 theSubLambdaTable->clearAndDestroy();
287 delete theSubLambdaTable;
288 }
289 }
290
291 delete modelManager;
292 (G4LossTableManager::Instance())->DeRegister(this);
293}
294
295//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296
[1055]297void G4VEnergyLossProcess::Clean()
[819]298{
[961]299 if(1 < verboseLevel) {
[819]300 G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName()
[961]301 << G4endl;
302 }
[819]303 delete [] theDEDXAtMaxEnergy;
304 delete [] theRangeAtMaxEnergy;
305 delete [] theEnergyOfCrossSectionMax;
306 delete [] theCrossSectionMax;
307 delete [] idxSCoffRegions;
[1055]308 delete [] idxDERegions;
[819]309
310 theDEDXAtMaxEnergy = 0;
311 theRangeAtMaxEnergy = 0;
312 theEnergyOfCrossSectionMax = 0;
313 theCrossSectionMax = 0;
314 tablesAreBuilt = false;
315
[961]316 //scTracks.clear();
[819]317 scProcesses.clear();
[961]318 nProcesses = 0;
[819]319}
320
321//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
322
[1055]323G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*,
324 const G4Material*,
325 G4double cut)
326{
327 return cut;
328}
329
330//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
331
332void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* p,
333 G4VEmFluctuationModel* fluc,
334 const G4Region* region)
335{
336 modelManager->AddEmModel(order, p, fluc, region);
337 if(p) p->SetParticleChange(pParticleChange, fluc);
338}
339
340//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
341
342void G4VEnergyLossProcess::UpdateEmModel(const G4String& nam,
343 G4double emin, G4double emax)
344{
345 modelManager->UpdateEmModel(nam, emin, emax);
346}
347
348//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
349
350void G4VEnergyLossProcess::SetEmModel(G4VEmModel* p, G4int index)
351{
352 G4int n = emModels.size();
353 if(index >= n) for(G4int i=n; i<index+1; i++) {emModels.push_back(0);}
354 emModels[index] = p;
355}
356
357//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
358
359G4VEmModel* G4VEnergyLossProcess::EmModel(G4int index)
360{
361 G4VEmModel* p = 0;
362 if(index >= 0 && index < G4int(emModels.size())) p = emModels[index];
363 return p;
364}
365
366//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
367
368G4VEmModel* G4VEnergyLossProcess::GetModelByIndex(G4int idx, G4bool ver)
369{
370 return modelManager->GetModel(idx, ver);
371}
372
373//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
374
375G4int G4VEnergyLossProcess::NumberOfModels()
376{
377 return modelManager->NumberOfModels();
378}
379
380//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381
[961]382void
383G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
384{
[819]385 if(1 < verboseLevel) {
386 G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for "
387 << GetProcessName()
388 << " for " << part.GetParticleName()
389 << G4endl;
390 }
391
[961]392 currentCouple = 0;
393 preStepLambda = 0.0;
394 mfpKinEnergy = DBL_MAX;
395 fRange = DBL_MAX;
396 preStepKinEnergy = 0.0;
397 chargeSqRatio = 1.0;
398 massRatio = 1.0;
399 reduceFactor = 1.0;
400
[819]401 G4LossTableManager* lManager = G4LossTableManager::Instance();
402
[961]403 // Are particle defined?
404 if( !particle ) {
405 particle = &part;
406 if(part.GetParticleType() == "nucleus") {
407 if(!theGenericIon) theGenericIon = G4GenericIon::GenericIon();
408 if(particle == theGenericIon) { isIon = true; }
409 else if(part.GetPDGCharge() > eplus) {
410 isIon = true;
411
412 // generic ions created on-fly
413 if(part.GetPDGCharge() > 2.5*eplus) {
414 particle = theGenericIon;
415 }
416 }
417 }
418 }
419
420 if( particle != &part) {
421 if(part.GetParticleType() == "nucleus") {
422 isIon = true;
423 lManager->RegisterIon(&part, this);
424 } else {
425 lManager->RegisterExtraParticle(&part, this);
426 }
[819]427 return;
428 }
429
[1055]430 Clean();
431 lManager->EmConfigurator()->AddModels();
[819]432
433 // Base particle and set of models can be defined here
434 InitialiseEnergyLossProcess(particle, baseParticle);
435
436 // Tables preparation
437 if (!baseParticle) {
438
439 theDEDXTable = G4PhysicsTableHelper::PreparePhysicsTable(theDEDXTable);
440 if (lManager->BuildCSDARange()) {
441 theDEDXunRestrictedTable =
442 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXunRestrictedTable);
443 theCSDARangeTable =
444 G4PhysicsTableHelper::PreparePhysicsTable(theCSDARangeTable);
445 }
446
447 theRangeTableForLoss =
448 G4PhysicsTableHelper::PreparePhysicsTable(theRangeTableForLoss);
449 theInverseRangeTable =
450 G4PhysicsTableHelper::PreparePhysicsTable(theInverseRangeTable);
451
452 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
453 if (nSCoffRegions) {
454 theDEDXSubTable =
455 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXSubTable);
456 theSubLambdaTable =
457 G4PhysicsTableHelper::PreparePhysicsTable(theSubLambdaTable);
458 }
459 }
460
461 G4double initialCharge = particle->GetPDGCharge();
462 G4double initialMass = particle->GetPDGMass();
463
464 if (baseParticle) {
465 massRatio = (baseParticle->GetPDGMass())/initialMass;
466 G4double q = initialCharge/baseParticle->GetPDGCharge();
467 chargeSqRatio = q*q;
468 if(chargeSqRatio > 0.0) reduceFactor = 1.0/(chargeSqRatio*massRatio);
469 }
470
471 theCuts = modelManager->Initialise(particle, secondaryParticle,
472 minSubRange, verboseLevel);
473
[1055]474 // Sub Cutoff and Deexcitation
475 if (nSCoffRegions>0 || nDERegions>0) {
[819]476 theSubCuts = modelManager->SubCutoff();
477
478 const G4ProductionCutsTable* theCoupleTable=
479 G4ProductionCutsTable::GetProductionCutsTable();
480 size_t numOfCouples = theCoupleTable->GetTableSize();
[1055]481
482 if(nSCoffRegions>0) idxSCoffRegions = new G4bool[numOfCouples];
483 if(nDERegions>0) idxDERegions = new G4bool[numOfCouples];
[819]484
485 for (size_t j=0; j<numOfCouples; j++) {
486
487 const G4MaterialCutsCouple* couple =
488 theCoupleTable->GetMaterialCutsCouple(j);
489 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
[1055]490
491 if(nSCoffRegions>0) {
492 G4bool reg = false;
493 for(G4int i=0; i<nSCoffRegions; i++) {
494 if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true;
495 }
496 idxSCoffRegions[j] = reg;
[819]497 }
[1055]498 if(nDERegions>0) {
499 G4bool reg = false;
500 for(G4int i=0; i<nDERegions; i++) {
501 if( pcuts == deRegions[i]->GetProductionCuts()) reg = true;
502 }
503 idxDERegions[j] = reg;
504 }
[819]505 }
506 }
507
508 lManager->EnergyLossProcessIsInitialised(particle, this);
509
510 if (1 < verboseLevel) {
511 G4cout << "G4VEnergyLossProcess::Initialise() is done "
[961]512 << " for local " << particle->GetParticleName()
513 << " isIon= " << isIon
[819]514 << " chargeSqRatio= " << chargeSqRatio
515 << " massRatio= " << massRatio
516 << " reduceFactor= " << reduceFactor << G4endl;
517 if (nSCoffRegions) {
518 G4cout << " SubCutoff Regime is ON for regions: " << G4endl;
519 for (G4int i=0; i<nSCoffRegions; i++) {
520 const G4Region* r = scoffRegions[i];
521 G4cout << " " << r->GetName() << G4endl;
522 }
523 }
[1055]524 if (nDERegions) {
525 G4cout << " Deexcitation is ON for regions: " << G4endl;
526 for (G4int i=0; i<nDERegions; i++) {
527 const G4Region* r = deRegions[i];
528 G4cout << " " << r->GetName() << G4endl;
529 }
530 }
[819]531 }
532}
533
534//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
535
536void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part)
537{
538 if(1 < verboseLevel) {
539 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for "
540 << GetProcessName()
541 << " and particle " << part.GetParticleName()
542 << "; local: " << particle->GetParticleName();
543 if(baseParticle) G4cout << "; base: " << baseParticle->GetParticleName();
544 G4cout << G4endl;
545 }
546
[961]547 if(&part == particle) {
548 if(!tablesAreBuilt) {
549 G4LossTableManager::Instance()->BuildPhysicsTable(particle, this);
550 }
551 if(!baseParticle) {
552 if(0 < verboseLevel) PrintInfoDefinition();
553
554 // needs to be done only once
555 safetyHelper->InitialiseHelper();
556 }
[819]557 }
558
559 if(1 < verboseLevel) {
560 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for "
561 << GetProcessName()
562 << " and particle " << part.GetParticleName();
563 if(isIonisation) G4cout << " isIonisation flag = 1";
564 G4cout << G4endl;
565 }
566}
567
568//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
569
570G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable(G4EmTableType tType)
571{
572 if(1 < verboseLevel) {
573 G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType
574 << " for " << GetProcessName()
575 << " and particle " << particle->GetParticleName()
576 << G4endl;
577 }
578 G4PhysicsTable* table = 0;
579 G4double emin = minKinEnergy;
580 G4double emax = maxKinEnergy;
581 G4int bin = nBins;
582
583 if(fTotal == tType) {
584 emax = maxKinEnergyCSDA;
585 bin = nBinsCSDA;
586 table = theDEDXunRestrictedTable;
587 } else if(fRestricted == tType) {
588 table = theDEDXTable;
589 if(theIonisationTable)
590 table = G4PhysicsTableHelper::PreparePhysicsTable(theIonisationTable);
591 } else if(fSubRestricted == tType) {
592 table = theDEDXSubTable;
593 if(theIonisationSubTable)
594 table = G4PhysicsTableHelper::PreparePhysicsTable(theIonisationSubTable);
595 } else {
596 G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type "
597 << tType << G4endl;
598 }
599
600 // Access to materials
601 const G4ProductionCutsTable* theCoupleTable=
602 G4ProductionCutsTable::GetProductionCutsTable();
603 size_t numOfCouples = theCoupleTable->GetTableSize();
604
605 if(1 < verboseLevel) {
606 G4cout << numOfCouples << " materials"
607 << " minKinEnergy= " << minKinEnergy
608 << " maxKinEnergy= " << maxKinEnergy
609 << " EmTableType= " << tType
610 << " table= " << table
611 << G4endl;
612 }
613 if(!table) return table;
614
615 for(size_t i=0; i<numOfCouples; i++) {
616
[961]617 if(1 < verboseLevel) {
[819]618 G4cout << "G4VEnergyLossProcess::BuildDEDXVector flag= "
619 << table->GetFlag(i) << G4endl;
[961]620 }
[819]621 if (table->GetFlag(i)) {
622
623 // create physics vector and fill it
624 const G4MaterialCutsCouple* couple =
625 theCoupleTable->GetMaterialCutsCouple(i);
626 G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin);
[961]627 aVector->SetSpline((G4LossTableManager::Instance())->SplineFlag());
628
[819]629 modelManager->FillDEDXVector(aVector, couple, tType);
630
631 // Insert vector for this material into the table
632 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
633 }
634 }
635
636 if(1 < verboseLevel) {
637 G4cout << "G4VEnergyLossProcess::BuildDEDXTable(): table is built for "
638 << particle->GetParticleName()
639 << " and process " << GetProcessName()
640 << G4endl;
641 // if(2 < verboseLevel) G4cout << (*table) << G4endl;
642 }
643
644 return table;
645}
646
647//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
648
649G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable(G4EmTableType tType)
650{
651 G4PhysicsTable* table = 0;
652
653 if(fRestricted == tType) {
654 table = theLambdaTable;
655 } else if(fSubRestricted == tType) {
656 table = theSubLambdaTable;
657 } else {
658 G4cout << "G4VEnergyLossProcess::BuildLambdaTable WARNING: wrong type "
659 << tType << G4endl;
660 }
661
662 if(1 < verboseLevel) {
663 G4cout << "G4VEnergyLossProcess::BuildLambdaTable() of type "
664 << tType << " for process "
665 << GetProcessName() << " and particle "
666 << particle->GetParticleName()
667 << " EmTableType= " << tType
668 << " table= " << table
669 << G4endl;
670 }
[961]671 if(!table) {return table;}
[819]672
673 // Access to materials
674 const G4ProductionCutsTable* theCoupleTable=
675 G4ProductionCutsTable::GetProductionCutsTable();
676 size_t numOfCouples = theCoupleTable->GetTableSize();
677
678 for(size_t i=0; i<numOfCouples; i++) {
679
680 if (table->GetFlag(i)) {
681
682 // create physics vector and fill it
683 const G4MaterialCutsCouple* couple =
684 theCoupleTable->GetMaterialCutsCouple(i);
685 G4double cut = (*theCuts)[i];
686 if(fSubRestricted == tType) cut = (*theSubCuts)[i];
687 G4PhysicsVector* aVector = LambdaPhysicsVector(couple, cut);
688 modelManager->FillLambdaVector(aVector, couple, true, tType);
689
690 // Insert vector for this material into the table
691 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
692 }
693 }
694
695 if(1 < verboseLevel) {
696 G4cout << "Lambda table is built for "
697 << particle->GetParticleName()
698 << G4endl;
699 }
700
701 return table;
702}
703
704//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
705
[1055]706void G4VEnergyLossProcess::PrintInfoDefinition()
[819]707{
[1055]708 if(0 < verboseLevel) {
709 G4cout << G4endl << GetProcessName() << ": for "
710 << particle->GetParticleName()
711 << " SubType= " << GetProcessSubType()
712 << G4endl
713 << " dE/dx and range tables from "
714 << G4BestUnit(minKinEnergy,"Energy")
715 << " to " << G4BestUnit(maxKinEnergy,"Energy")
716 << " in " << nBins << " bins" << G4endl
717 << " Lambda tables from threshold to "
718 << G4BestUnit(maxKinEnergy,"Energy")
719 << " in " << nBins << " bins, spline: "
720 << (G4LossTableManager::Instance())->SplineFlag()
721 << G4endl;
722 if(theRangeTableForLoss && isIonisation) {
723 G4cout << " finalRange(mm)= " << finalRange/mm
724 << ", dRoverRange= " << dRoverRange
725 << ", integral: " << integral
726 << ", fluct: " << lossFluctuationFlag
727 << ", linLossLimit= " << linLossLimit
728 << G4endl;
729 }
730 PrintInfo();
731 modelManager->DumpModelList(verboseLevel);
732 if(theCSDARangeTable && isIonisation) {
733 G4cout << " CSDA range table up"
734 << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
735 << " in " << nBinsCSDA << " bins" << G4endl;
736 }
737 if(nSCoffRegions>0 && isIonisation) {
738 G4cout << " Subcutoff sampling in " << nSCoffRegions
739 << " regions" << G4endl;
740 }
741 if(2 < verboseLevel) {
742 G4cout << " DEDXTable address= " << theDEDXTable << G4endl;
743 if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
744 G4cout << "non restricted DEDXTable address= "
745 << theDEDXunRestrictedTable << G4endl;
746 if(theDEDXunRestrictedTable && isIonisation) {
747 G4cout << (*theDEDXunRestrictedTable) << G4endl;
748 }
749 if(theDEDXSubTable && isIonisation) {
750 G4cout << (*theDEDXSubTable) << G4endl;
751 }
752 G4cout << " CSDARangeTable address= " << theCSDARangeTable
753 << G4endl;
754 if(theCSDARangeTable && isIonisation) {
755 G4cout << (*theCSDARangeTable) << G4endl;
756 }
757 G4cout << " RangeTableForLoss address= " << theRangeTableForLoss
758 << G4endl;
759 if(theRangeTableForLoss && isIonisation) {
760 G4cout << (*theRangeTableForLoss) << G4endl;
761 }
762 G4cout << " InverseRangeTable address= " << theInverseRangeTable
763 << G4endl;
764 if(theInverseRangeTable && isIonisation) {
765 G4cout << (*theInverseRangeTable) << G4endl;
766 }
767 G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
768 if(theLambdaTable && isIonisation) {
769 G4cout << (*theLambdaTable) << G4endl;
770 }
771 G4cout << " SubLambdaTable address= " << theSubLambdaTable << G4endl;
772 if(theSubLambdaTable && isIonisation) {
773 G4cout << (*theSubLambdaTable) << G4endl;
774 }
775 }
776 }
[819]777}
778
779//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
780
[1055]781void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
782{
783 G4RegionStore* regionStore = G4RegionStore::GetInstance();
784 const G4Region* reg = r;
785 if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
786
787 // the region is in the list
788 if (nSCoffRegions) {
789 for (G4int i=0; i<nSCoffRegions; i++) {
790 if (reg == scoffRegions[i]) {
791 if(!val) deRegions[i] = 0;
792 return;
793 }
794 }
795 }
796
797 // new region
798 if(val) {
799 useSubCutoff = true;
800 scoffRegions.push_back(reg);
801 nSCoffRegions++;
802 } else {
803 useSubCutoff = false;
804 }
805}
806
807//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
808
809void G4VEnergyLossProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
810{
811 G4RegionStore* regionStore = G4RegionStore::GetInstance();
812 const G4Region* reg = r;
813 if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
814
815 // the region is in the list
816 if (nDERegions) {
817 for (G4int i=0; i<nDERegions; i++) {
818 if (reg == deRegions[i]) {
819 if(!val) deRegions[i] = 0;
820 return;
821 }
822 }
823 }
824
825 // new region
826 if(val) {
827 useDeexcitation = true;
828 deRegions.push_back(reg);
829 nDERegions++;
830 } else {
831 useDeexcitation = false;
832 }
833}
834
835//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
836
[819]837G4double G4VEnergyLossProcess::AlongStepGetPhysicalInteractionLength(
838 const G4Track&,
839 G4double,
840 G4double currentMinStep,
841 G4double&,
842 G4GPILSelection* selection)
843{
844 G4double x = DBL_MAX;
845 *selection = aGPILSelection;
846 if(isIonisation) {
847 fRange = GetScaledRangeForScaledEnergy(preStepScaledEnergy)*reduceFactor;
848
849 x = fRange;
850 G4double y = x*dRoverRange;
851
852 if(x > finalRange && y < currentMinStep) {
853 x = y + finalRange*(1.0 - dRoverRange)*(2.0 - finalRange/fRange);
[961]854 } else if (rndmStepFlag) {x = SampleRange();}
855 //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
856 // <<" range= "<<fRange <<" cMinSt="<<currentMinStep
857 // << " limit= " << x <<G4endl;
[819]858 }
[961]859 //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
[819]860 // <<" stepLimit= "<<x<<G4endl;
861 return x;
862}
863
864//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
865
866G4double G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength(
867 const G4Track& track,
868 G4double previousStepSize,
869 G4ForceCondition* condition)
870{
871 // condition is set to "Not Forced"
872 *condition = NotForced;
873 G4double x = DBL_MAX;
[961]874
875 // initialisation of material, mass, charge, model at the beginning of the step
876 DefineMaterial(track.GetMaterialCutsCouple());
877
878 const G4ParticleDefinition* currPart = track.GetDefinition();
879 if(theGenericIon == particle) {
880 massRatio = proton_mass_c2/currPart->GetPDGMass();
881 }
882 preStepKinEnergy = track.GetKineticEnergy();
883 preStepScaledEnergy = preStepKinEnergy*massRatio;
884 SelectModel(preStepScaledEnergy);
885
886 if(isIon) {
887 chargeSqRatio =
888 currentModel->GetChargeSquareRatio(currPart,currentMaterial,preStepKinEnergy);
889 reduceFactor = 1.0/(chargeSqRatio*massRatio);
890 }
891 //G4cout << "q2= " << chargeSqRatio << " massRatio= " << massRatio << G4endl;
892 // initialisation for sampling of the interaction length
[819]893 if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
[961]894 if(theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
[819]895
[961]896 // compute mean free path
[819]897 if(preStepScaledEnergy < mfpKinEnergy) {
898 if (integral) ComputeLambdaForScaledEnergy(preStepScaledEnergy);
899 else preStepLambda = GetLambdaForScaledEnergy(preStepScaledEnergy);
900 if(preStepLambda <= DBL_MIN) mfpKinEnergy = 0.0;
901 }
902
903 // non-zero cross section
904 if(preStepLambda > DBL_MIN) {
905 if (theNumberOfInteractionLengthLeft < 0.0) {
906 // beggining of tracking (or just after DoIt of this process)
[961]907 //G4cout<<"G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength Reset"<<G4endl;
[819]908 ResetNumberOfInteractionLengthLeft();
909 } else if(currentInteractionLength < DBL_MAX) {
910 // subtract NumberOfInteractionLengthLeft
911 SubtractNumberOfInteractionLengthLeft(previousStepSize);
912 if(theNumberOfInteractionLengthLeft < 0.)
913 theNumberOfInteractionLengthLeft = perMillion;
914 }
915
916 // get mean free path and step limit
917 currentInteractionLength = 1.0/preStepLambda;
918 x = theNumberOfInteractionLengthLeft * currentInteractionLength;
919
920#ifdef G4VERBOSE
921 if (verboseLevel>2){
922 G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
923 G4cout << "[ " << GetProcessName() << "]" << G4endl;
[961]924 G4cout << " for " << currPart->GetParticleName()
[819]925 << " in Material " << currentMaterial->GetName()
926 << " Ekin(MeV)= " << preStepKinEnergy/MeV
927 <<G4endl;
928 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]"
929 << "InteractionLength= " << x/cm <<"[cm] " <<G4endl;
930 }
931#endif
932 // zero cross section case
933 } else {
934 if(theNumberOfInteractionLengthLeft > DBL_MIN &&
935 currentInteractionLength < DBL_MAX) {
936 // subtract NumberOfInteractionLengthLeft
937 SubtractNumberOfInteractionLengthLeft(previousStepSize);
938 if(theNumberOfInteractionLengthLeft < 0.)
939 theNumberOfInteractionLengthLeft = perMillion;
940 }
941 currentInteractionLength = DBL_MAX;
942 }
943 return x;
944}
945
946//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
947
948G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track,
949 const G4Step& step)
950{
951 fParticleChange.InitializeForAlongStep(track);
952 // The process has range table - calculate energy loss
953 if(!isIonisation) return &fParticleChange;
954
955 // Get the actual (true) Step length
956 G4double length = step.GetStepLength();
957 if(length <= DBL_MIN) return &fParticleChange;
958 G4double eloss = 0.0;
[961]959 G4double esecdep = 0.0;
960
961 /*
[819]962 if(-1 < verboseLevel) {
963 const G4ParticleDefinition* d = track.GetDefinition();
964 G4cout << "AlongStepDoIt for "
965 << GetProcessName() << " and particle "
966 << d->GetParticleName()
967 << " eScaled(MeV)= " << preStepScaledEnergy/MeV
968 << " range(mm)= " << fRange/mm
969 << " s(mm)= " << length/mm
970 << " q^2= " << chargeSqRatio
971 << " md= " << d->GetPDGMass()
972 << " status= " << track.GetTrackStatus()
973 << G4endl;
974 }
975 */
976
[961]977 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
978
[819]979 // stopping
980 if (length >= fRange) {
[961]981 eloss = preStepKinEnergy;
[1055]982 if (useDeexcitation) {
983 if(idxDERegions[currentMaterialIndex]) {
984 currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
985 if(eloss < 0.0) eloss = 0.0;
986 }
987 }
[819]988 fParticleChange.SetProposedKineticEnergy(0.0);
[1055]989 fParticleChange.ProposeLocalEnergyDeposit(eloss);
[819]990 return &fParticleChange;
991 }
992
993 // Short step
994 eloss = GetDEDXForScaledEnergy(preStepScaledEnergy)*length;
995
996 // Long step
997 //} else {
998 if(eloss > preStepKinEnergy*linLossLimit) {
999
1000 G4double x =
1001 GetScaledRangeForScaledEnergy(preStepScaledEnergy) - length/reduceFactor;
1002 eloss = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio;
[961]1003
[819]1004 /*
1005 if(-1 < verboseLevel)
1006 G4cout << "Long STEP: rPre(mm)= "
1007 << GetScaledRangeForScaledEnergy(preStepScaledEnergy)/mm
1008 << " rPost(mm)= " << x/mm
1009 << " ePre(MeV)= " << preStepScaledEnergy/MeV
1010 << " eloss(MeV)= " << eloss/MeV
1011 << " eloss0(MeV)= "
1012 << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV
[961]1013 << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV
[819]1014 << G4endl;
1015 */
1016 }
1017
[961]1018 /*
[819]1019 G4double eloss0 = eloss;
1020 if(-1 < verboseLevel ) {
1021 G4cout << "Before fluct: eloss(MeV)= " << eloss/MeV
1022 << " e-eloss= " << preStepKinEnergy-eloss
1023 << " step(mm)= " << length/mm
1024 << " range(mm)= " << fRange/mm
1025 << " fluct= " << lossFluctuationFlag
1026 << G4endl;
1027 }
1028 */
1029
1030 G4double cut = (*theCuts)[currentMaterialIndex];
1031 G4double esec = 0.0;
1032
1033 // SubCutOff
1034 if(useSubCutoff) {
1035 if(idxSCoffRegions[currentMaterialIndex]) {
1036
[961]1037 G4bool yes = false;
[819]1038 G4StepPoint* prePoint = step.GetPreStepPoint();
1039
[961]1040 // Check boundary
1041 if(prePoint->GetStepStatus() == fGeomBoundary) yes = true;
[819]1042
[961]1043 // Check PrePoint
1044 else {
1045 G4double preSafety = prePoint->GetSafety();
1046 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
[819]1047
[961]1048 // recompute presafety
1049 if(preSafety < rcut) {
[819]1050 preSafety = safetyHelper->ComputeSafety(prePoint->GetPosition());
[961]1051 }
[819]1052
[961]1053 if(preSafety < rcut) yes = true;
1054
1055 // Check PostPoint
1056 else {
1057 G4double postSafety = preSafety - length;
1058 if(postSafety < rcut) {
1059 postSafety =
1060 safetyHelper->ComputeSafety(step.GetPostStepPoint()->GetPosition());
1061 if(postSafety < rcut) yes = true;
1062 }
1063 }
[819]1064 }
[961]1065
[819]1066 // Decide to start subcut sampling
[961]1067 if(yes) {
[819]1068
1069 cut = (*theSubCuts)[currentMaterialIndex];
1070 eloss -= GetSubDEDXForScaledEnergy(preStepScaledEnergy)*length;
1071 scTracks.clear();
1072 SampleSubCutSecondaries(scTracks, step,
1073 currentModel,currentMaterialIndex,
1074 esecdep);
[961]1075 // add bremsstrahlung sampling
[819]1076 /*
1077 if(nProcesses > 0) {
1078 for(G4int i=0; i<nProcesses; i++) {
1079 (scProcesses[i])->SampleSubCutSecondaries(
1080 scTracks, step, (scProcesses[i])->
1081 SelectModelForMaterial(preStepKinEnergy, currentMaterialIndex),
1082 currentMaterialIndex,esecdep);
1083 }
1084 }
1085 */
1086 G4int n = scTracks.size();
1087 if(n>0) {
1088 G4ThreeVector mom = dynParticle->GetMomentum();
1089 fParticleChange.SetNumberOfSecondaries(n);
1090 for(G4int i=0; i<n; i++) {
1091 G4Track* t = scTracks[i];
1092 G4double e = t->GetKineticEnergy();
1093 if (t->GetDefinition() == thePositron) e += 2.0*electron_mass_c2;
1094 esec += e;
1095 pParticleChange->AddSecondary(t);
1096 }
1097 }
1098 }
1099 }
1100 }
1101
1102 // Corrections, which cannot be tabulated
[961]1103 currentModel->CorrectionsAlongStep(currentCouple, dynParticle,
1104 eloss, esecdep, length);
[819]1105
1106 // Sample fluctuations
1107 if (lossFluctuationFlag) {
1108 G4VEmFluctuationModel* fluc = currentModel->GetModelOfFluctuations();
1109 if(fluc &&
1110 (eloss + esec + esecdep + lowestKinEnergy) < preStepKinEnergy) {
1111
1112 G4double tmax =
1113 std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
[1055]1114 G4double emean = eloss;
[819]1115 eloss = fluc->SampleFluctuations(currentMaterial,dynParticle,
[1055]1116 tmax,length,emean);
[961]1117 /*
[819]1118 if(-1 < verboseLevel)
1119 G4cout << "After fluct: eloss(MeV)= " << eloss/MeV
1120 << " fluc= " << (eloss-eloss0)/MeV
[961]1121 << " ChargeSqRatio= " << chargeSqRatio
[819]1122 << " massRatio= " << massRatio
1123 << " tmax= " << tmax
1124 << G4endl;
1125 */
1126 }
1127 }
1128 // add low-energy subcutoff particles
1129 eloss += esecdep;
1130 if(eloss < 0.0) eloss = 0.0;
1131
[1055]1132 // deexcitation
1133 else if (useDeexcitation) {
1134 if(idxDERegions[currentMaterialIndex]) {
1135 currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
1136 if(eloss < 0.0) eloss = 0.0;
1137 }
1138 }
1139
[819]1140 // Energy balanse
1141 G4double finalT = preStepKinEnergy - eloss - esec;
1142 if (finalT <= lowestKinEnergy) {
1143 eloss = preStepKinEnergy - esec;
1144 finalT = 0.0;
[961]1145 } else if(isIon) {
1146 fParticleChange.SetProposedCharge(
1147 currentModel->GetParticleCharge(track.GetDefinition(),currentMaterial,finalT));
[819]1148 }
1149
1150 fParticleChange.SetProposedKineticEnergy(finalT);
1151 fParticleChange.ProposeLocalEnergyDeposit(eloss);
1152
[961]1153 /*
[819]1154 if(-1 < verboseLevel) {
1155 G4cout << "Final value eloss(MeV)= " << eloss/MeV
1156 << " preStepKinEnergy= " << preStepKinEnergy
1157 << " postStepKinEnergy= " << finalT
1158 << " lossFlag= " << lossFluctuationFlag
1159 << " status= " << track.GetTrackStatus()
1160 << G4endl;
1161 }
[961]1162 */
[819]1163
1164 return &fParticleChange;
1165}
1166
1167//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1168
1169void G4VEnergyLossProcess::SampleSubCutSecondaries(
1170 std::vector<G4Track*>& tracks,
1171 const G4Step& step,
1172 G4VEmModel* model,
1173 G4int idx,
[961]1174 G4double& /*extraEdep*/)
[819]1175{
1176 // Fast check weather subcutoff can work
1177 G4double subcut = (*theSubCuts)[idx];
1178 G4double cut = (*theCuts)[idx];
1179 if(cut <= subcut) return;
1180
1181 const G4Track* track = step.GetTrack();
1182 const G4DynamicParticle* dp = track->GetDynamicParticle();
[961]1183 G4double e = dp->GetKineticEnergy()*massRatio;
[819]1184 G4bool b;
[961]1185 G4double cross = chargeSqRatio*(((*theSubLambdaTable)[idx])->GetValue(e,b));
[819]1186 G4double length = step.GetStepLength();
1187
1188 // negligible probability to get any interaction
1189 if(length*cross < perMillion) return;
[961]1190 /*
[819]1191 if(-1 < verboseLevel)
1192 G4cout << "<<< Subcutoff for " << GetProcessName()
1193 << " cross(1/mm)= " << cross*mm << ">>>"
1194 << " e(MeV)= " << preStepScaledEnergy
1195 << " matIdx= " << currentMaterialIndex
1196 << G4endl;
1197 */
1198
1199 // Sample subcutoff secondaries
1200 G4StepPoint* preStepPoint = step.GetPreStepPoint();
[961]1201 G4StepPoint* postStepPoint = step.GetPostStepPoint();
[819]1202 G4ThreeVector prepoint = preStepPoint->GetPosition();
[961]1203 G4ThreeVector dr = postStepPoint->GetPosition() - prepoint;
[819]1204 G4double pretime = preStepPoint->GetGlobalTime();
[961]1205 G4double dt = postStepPoint->GetGlobalTime() - pretime;
1206 //G4double dt = length/preStepPoint->GetVelocity();
[819]1207 G4double fragment = 0.0;
1208
1209 do {
1210 G4double del = -std::log(G4UniformRand())/cross;
1211 fragment += del/length;
1212 if (fragment > 1.0) break;
1213
1214 // sample secondaries
1215 secParticles.clear();
1216 model->SampleSecondaries(&secParticles,track->GetMaterialCutsCouple(),
1217 dp,subcut,cut);
1218
1219 // position of subcutoff particles
1220 G4ThreeVector r = prepoint + fragment*dr;
1221 std::vector<G4DynamicParticle*>::iterator it;
1222 for(it=secParticles.begin(); it!=secParticles.end(); it++) {
1223
1224 G4bool addSec = true;
[961]1225 /*
[819]1226 // do not track very low-energy delta-electrons
1227 if(theSecondaryRangeTable && (*it)->GetDefinition() == theElectron) {
1228 G4bool b;
1229 G4double ekin = (*it)->GetKineticEnergy();
1230 G4double rg = ((*theSecondaryRangeTable)[idx]->GetValue(ekin, b));
1231 // if(rg < currentMinSafety) {
1232 if(rg < safetyHelper->ComputeSafety(r)) {
1233 extraEdep += ekin;
1234 delete (*it);
1235 addSec = false;
1236 }
1237 }
[961]1238 */
[819]1239 if(addSec) {
[961]1240 G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
1241 //G4Track* t = new G4Track((*it), pretime, r);
[819]1242 t->SetTouchableHandle(track->GetTouchableHandle());
1243 tracks.push_back(t);
1244
[961]1245 /*
1246 if(-1 < verboseLevel)
1247 G4cout << "New track " << t->GetDefinition()->GetParticleName()
1248 << " e(keV)= " << t->GetKineticEnergy()/keV
1249 << " fragment= " << fragment
1250 << G4endl;
[819]1251 */
1252 }
1253 }
1254 } while (fragment <= 1.0);
1255}
1256
1257//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1258
1259G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track,
1260 const G4Step&)
1261{
1262 fParticleChange.InitializeForPostStep(track);
1263 G4double finalT = track.GetKineticEnergy();
1264 if(finalT <= lowestKinEnergy) return &fParticleChange;
1265
1266 G4double postStepScaledEnergy = finalT*massRatio;
1267 /*
1268 if(-1 < verboseLevel) {
1269 G4cout << GetProcessName()
1270 << "::PostStepDoIt: E(MeV)= " << finalT/MeV
1271 << G4endl;
1272 }
1273 */
1274 // Integral approach
1275 if (integral) {
1276 G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy);
1277 /*
1278 if(preStepLambda<lx && 1 < verboseLevel && nWarnings<200) {
1279 G4cout << "WARNING: for " << particle->GetParticleName()
1280 << " and " << GetProcessName()
1281 << " E(MeV)= " << finalT/MeV
1282 << " preLambda= " << preStepLambda
1283 << " < " << lx << " (postLambda) "
1284 << G4endl;
1285 nWarnings++;
1286 }
1287 */
1288 if(preStepLambda*G4UniformRand() > lx) {
1289 ClearNumberOfInteractionLengthLeft();
1290 return &fParticleChange;
1291 }
1292 }
1293
[961]1294 SelectModel(postStepScaledEnergy);
[1055]1295 if(useDeexcitation) {
1296 currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
1297 }
1298
[819]1299 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
1300 G4double tcut = (*theCuts)[currentMaterialIndex];
1301
1302 // sample secondaries
1303 secParticles.clear();
1304 currentModel->SampleSecondaries(&secParticles, currentCouple, dynParticle, tcut);
1305
1306 // save secondaries
1307 G4int num = secParticles.size();
1308 if(num > 0) {
1309 fParticleChange.SetNumberOfSecondaries(num);
1310 for (G4int i=0; i<num; i++) {
1311 fParticleChange.AddSecondary(secParticles[i]);
1312 }
1313 }
1314
1315 /*
1316 if(-1 < verboseLevel) {
1317 G4cout << "::PostStepDoIt: Sample secondary; Efin= "
1318 << fParticleChange.GetProposedKineticEnergy()/MeV
1319 << " MeV; model= (" << currentModel->LowEnergyLimit()
1320 << ", " << currentModel->HighEnergyLimit() << ")"
1321 << " preStepLambda= " << preStepLambda
1322 << " dir= " << track.GetMomentumDirection()
1323 << " status= " << track.GetTrackStatus()
1324 << G4endl;
1325 }
1326 */
1327 ClearNumberOfInteractionLengthLeft();
1328 return &fParticleChange;
1329}
1330
1331//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1332
[961]1333G4bool G4VEnergyLossProcess::StorePhysicsTable(
1334 const G4ParticleDefinition* part, const G4String& directory,
1335 G4bool ascii)
[819]1336{
[961]1337 G4bool res = true;
1338 if ( baseParticle || part != particle ) return res;
1339
1340 if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX"))
1341 {res = false;}
1342
1343 if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr"))
1344 {res = false;}
1345
1346 if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX"))
1347 {res = false;}
1348
1349 if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation"))
1350 {res = false;}
1351
1352 if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation"))
1353 {res = false;}
1354
1355 if(isIonisation &&
1356 !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange"))
1357 {res = false;}
1358
1359 if(isIonisation &&
1360 !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range"))
1361 {res = false;}
1362
1363 if(isIonisation &&
1364 !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange"))
1365 {res = false;}
1366
1367 if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda"))
1368 {res = false;}
1369
1370 if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda"))
1371 {res = false;}
1372
1373 if ( res ) {
1374 if(0 < verboseLevel) {
1375 G4cout << "Physics tables are stored for " << particle->GetParticleName()
1376 << " and process " << GetProcessName()
1377 << " in the directory <" << directory
1378 << "> " << G4endl;
1379 }
1380 } else {
1381 G4cout << "Fail to store Physics Tables for "
1382 << particle->GetParticleName()
1383 << " and process " << GetProcessName()
1384 << " in the directory <" << directory
1385 << "> " << G4endl;
1386 }
1387 return res;
1388}
1389
1390//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1391
1392G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
1393 const G4ParticleDefinition* part, const G4String& directory,
1394 G4bool ascii)
1395{
1396 G4bool res = true;
1397 const G4String particleName = part->GetParticleName();
1398
1399 if(1 < verboseLevel) {
1400 G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
1401 << particleName << " and process " << GetProcessName()
1402 << "; tables_are_built= " << tablesAreBuilt
[819]1403 << G4endl;
[961]1404 }
1405 if(particle == part) {
[819]1406
[961]1407 // G4bool yes = true;
1408 if ( !baseParticle ) {
1409
1410 G4bool fpi = true;
1411 if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi))
1412 {fpi = false;}
1413
1414 if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false))
1415 {fpi = false;}
1416
1417 if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi))
1418 {res = false;}
1419
1420 if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false))
1421 {res = false;}
1422
1423 if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false))
1424 {res = false;}
1425
1426 if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi))
1427 {res = false;}
1428
1429 if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true))
1430 {res = false;}
1431
1432 G4bool yes = false;
1433 if(nSCoffRegions > 0) {yes = true;}
1434
1435 if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes))
1436 {res = false;}
1437
1438 if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes))
1439 {res = false;}
1440
1441 if(!fpi) yes = false;
1442 if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes))
1443 {res = false;}
1444 }
1445 }
1446
1447 return res;
1448}
1449
1450//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1451
[1055]1452G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part,
1453 G4PhysicsTable* aTable, G4bool ascii,
1454 const G4String& directory,
1455 const G4String& tname)
1456{
1457 G4bool res = true;
1458 if ( aTable ) {
1459 const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
1460 if( !aTable->StorePhysicsTable(name,ascii)) res = false;
1461 }
1462 return res;
1463}
1464
1465//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1466
[961]1467G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part,
1468 G4PhysicsTable* aTable, G4bool ascii,
1469 const G4String& directory,
1470 const G4String& tname,
1471 G4bool mandatory)
1472{
1473 G4bool res = true;
1474 G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
1475 G4bool yes = aTable->ExistPhysicsTable(filename);
1476 if(yes) {
1477 yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
1478 if((G4LossTableManager::Instance())->SplineFlag()) {
1479 size_t n = aTable->length();
1480 for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
1481 }
1482 }
1483 if(yes) {
1484 if (0 < verboseLevel) {
1485 G4cout << tname << " table for " << part->GetParticleName()
1486 << " is Retrieved from <" << filename << ">"
[819]1487 << G4endl;
[961]1488 }
1489 } else {
1490 if(mandatory) res = false;
1491 if(mandatory || 1 < verboseLevel) {
1492 G4cout << tname << " table for " << part->GetParticleName()
1493 << " from file <"
1494 << filename << "> is not Retrieved"
[819]1495 << G4endl;
1496 }
1497 }
[961]1498 return res;
[819]1499}
1500
1501//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
[1055]1502
1503G4double G4VEnergyLossProcess::GetDEDXDispersion(
1504 const G4MaterialCutsCouple *couple,
1505 const G4DynamicParticle* dp,
1506 G4double length)
1507{
1508 DefineMaterial(couple);
1509 G4double ekin = dp->GetKineticEnergy();
1510 SelectModel(ekin*massRatio);
1511 G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
1512 tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
1513 G4double d = 0.0;
1514 G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
1515 if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
1516 return d;
1517}
1518
1519//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1520
1521G4double G4VEnergyLossProcess::CrossSectionPerVolume(
1522 G4double kineticEnergy, const G4MaterialCutsCouple* couple)
1523{
1524 // Cross section per volume is calculated
1525 DefineMaterial(couple);
1526 G4double cross = 0.0;
1527 G4bool b;
1528 if(theLambdaTable) {
1529 cross =
1530 ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
1531 } else {
1532 SelectModel(kineticEnergy);
1533 cross =
1534 currentModel->CrossSectionPerVolume(currentMaterial,
1535 particle, kineticEnergy,
1536 (*theCuts)[currentMaterialIndex]);
1537 }
1538
1539 return cross;
1540}
1541
1542//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1543
1544G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
1545{
1546 DefineMaterial(track.GetMaterialCutsCouple());
1547 preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
1548 G4double x = DBL_MAX;
1549 if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
1550 return x;
1551}
1552
1553//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1554
1555G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track,
1556 G4double x, G4double y,
1557 G4double& z)
1558{
1559 G4GPILSelection sel;
1560 return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
1561}
1562
1563//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1564
1565G4double G4VEnergyLossProcess::GetMeanFreePath(
1566 const G4Track& track,
1567 G4double,
1568 G4ForceCondition* condition)
1569
1570{
1571 *condition = NotForced;
1572 return MeanFreePath(track);
1573}
1574
1575//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1576
1577G4double G4VEnergyLossProcess::GetContinuousStepLimit(
1578 const G4Track&,
1579 G4double, G4double, G4double&)
1580{
1581 return DBL_MAX;
1582}
1583
1584//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1585
1586G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
1587 const G4MaterialCutsCouple* couple, G4double cut)
1588{
1589 // G4double cut = (*theCuts)[couple->GetIndex()];
1590 // G4int nbins = nLambdaBins;
1591 G4double tmin =
1592 std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
1593 minKinEnergy);
1594 if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
1595 G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
1596 v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
1597
1598 return v;
1599}
1600
1601//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
[961]1602
1603void G4VEnergyLossProcess::AddCollaborativeProcess(
1604 G4VEnergyLossProcess* p)
1605{
1606 G4bool add = true;
1607 if(p->GetProcessName() != "eBrem") add = false;
1608 if(add && nProcesses > 0) {
1609 for(G4int i=0; i<nProcesses; i++) {
1610 if(p == scProcesses[i]) {
1611 add = false;
1612 break;
1613 }
1614 }
1615 }
1616 if(add) {
1617 scProcesses.push_back(p);
1618 nProcesses++;
1619 if (1 < verboseLevel) {
1620 G4cout << "### The process " << p->GetProcessName()
1621 << " is added to the list of collaborative processes of "
1622 << GetProcessName() << G4endl;
1623 }
1624 }
1625}
1626
1627//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1628
[1055]1629void G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType)
[819]1630{
[1055]1631 if(fTotal == tType && theDEDXunRestrictedTable != p) {
1632 if(theDEDXunRestrictedTable) theDEDXunRestrictedTable->clearAndDestroy();
1633 theDEDXunRestrictedTable = p;
1634 if(p) {
1635 size_t n = p->length();
1636 G4PhysicsVector* pv = (*p)[0];
1637 G4double emax = maxKinEnergyCSDA;
1638 G4bool b;
1639 theDEDXAtMaxEnergy = new G4double [n];
1640
1641 for (size_t i=0; i<n; i++) {
1642 pv = (*p)[i];
1643 G4double dedx = pv->GetValue(emax, b);
1644 theDEDXAtMaxEnergy[i] = dedx;
1645 //G4cout << "i= " << i << " emax(MeV)= " << emax/MeV<< " dedx= "
1646 //<< dedx << G4endl;
1647 }
1648 }
1649
1650 } else if(fRestricted == tType) {
1651 theDEDXTable = p;
1652 } else if(fSubRestricted == tType) {
1653 theDEDXSubTable = p;
1654 } else if(fIsIonisation == tType && theIonisationTable != p) {
1655 if(theIonisationTable) theIonisationTable->clearAndDestroy();
1656 theIonisationTable = p;
1657 } else if(fIsSubIonisation == tType && theIonisationSubTable != p) {
1658 if(theIonisationSubTable) theIonisationSubTable->clearAndDestroy();
1659 theIonisationSubTable = p;
1660 }
[819]1661}
1662
1663//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1664
[1055]1665void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p)
1666{
1667 if(theCSDARangeTable != p) theCSDARangeTable = p;
[819]1668
[1055]1669 if(p) {
1670 size_t n = p->length();
1671 G4PhysicsVector* pv = (*p)[0];
1672 G4double emax = maxKinEnergyCSDA;
1673 G4bool b;
1674 theRangeAtMaxEnergy = new G4double [n];
1675
1676 for (size_t i=0; i<n; i++) {
1677 pv = (*p)[i];
1678 G4double r2 = pv->GetValue(emax, b);
1679 theRangeAtMaxEnergy[i] = r2;
1680 //G4cout << "i= " << i << " e2(MeV)= " << emax/MeV << " r2= "
1681 //<< r2<< G4endl;
1682 }
1683 }
1684}
1685
[819]1686//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1687
[1055]1688void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p)
1689{
1690 if(theRangeTableForLoss != p) {
1691 theRangeTableForLoss = p;
1692 if(1 < verboseLevel) {
1693 G4cout << "### Set Range table " << p
1694 << " for " << particle->GetParticleName()
1695 << " and process " << GetProcessName() << G4endl;
1696 }
1697 }
1698}
1699
1700//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1701
1702void G4VEnergyLossProcess::SetSecondaryRangeTable(G4PhysicsTable* p)
1703{
1704 if(theSecondaryRangeTable != p) {
1705 theSecondaryRangeTable = p;
1706 if(1 < verboseLevel) {
1707 G4cout << "### Set SecondaryRange table " << p
1708 << " for " << particle->GetParticleName()
1709 << " and process " << GetProcessName() << G4endl;
1710 }
1711 }
1712}
1713
1714//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1715
1716void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p)
1717{
1718 if(theInverseRangeTable != p) {
1719 theInverseRangeTable = p;
1720 if(1 < verboseLevel) {
1721 G4cout << "### Set InverseRange table " << p
1722 << " for " << particle->GetParticleName()
1723 << " and process " << GetProcessName() << G4endl;
1724 }
1725 }
1726}
1727
1728//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1729
1730void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p)
1731{
1732 if(1 < verboseLevel) {
1733 G4cout << "### Set Lambda table " << p
1734 << " for " << particle->GetParticleName()
1735 << " and process " << GetProcessName() << G4endl;
1736 }
1737 if(theLambdaTable != p) theLambdaTable = p;
1738 tablesAreBuilt = true;
1739
1740 if(p) {
1741 size_t n = p->length();
1742 G4PhysicsVector* pv = (*p)[0];
1743 G4double e, s, smax, emax;
1744 theEnergyOfCrossSectionMax = new G4double [n];
1745 theCrossSectionMax = new G4double [n];
1746 G4bool b;
1747
1748 for (size_t i=0; i<n; i++) {
1749 pv = (*p)[i];
1750 emax = DBL_MAX;
1751 smax = 0.0;
1752 if(pv) {
1753 size_t nb = pv->GetVectorLength();
1754 emax = pv->GetLowEdgeEnergy(nb);
1755 for (size_t j=0; j<nb; j++) {
1756 e = pv->GetLowEdgeEnergy(j);
1757 s = pv->GetValue(e,b);
1758 if(s > smax) {
1759 smax = s;
1760 emax = e;
1761 }
1762 }
1763 }
1764 theEnergyOfCrossSectionMax[i] = emax;
1765 theCrossSectionMax[i] = smax;
1766 if(1 < verboseLevel) {
1767 G4cout << "For " << particle->GetParticleName()
1768 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
1769 << " lambda= " << smax << G4endl;
1770 }
1771 }
1772 }
1773}
1774
1775//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1776
1777void G4VEnergyLossProcess::SetSubLambdaTable(G4PhysicsTable* p)
1778{
1779 if(theSubLambdaTable != p) {
1780 theSubLambdaTable = p;
1781 if(1 < verboseLevel) {
1782 G4cout << "### Set SebLambda table " << p
1783 << " for " << particle->GetParticleName()
1784 << " and process " << GetProcessName() << G4endl;
1785 }
1786 }
1787}
1788
1789//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1790
Note: See TracBrowser for help on using the repository browser.