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

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

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