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

Last change on this file since 989 was 961, checked in by garnier, 15 years ago

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