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

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

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