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

Last change on this file since 1015 was 1007, checked in by garnier, 15 years ago

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