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

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

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