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

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

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