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

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

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