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

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

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