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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

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