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

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

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