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

Last change on this file since 1129 was 1055, checked in by garnier, 15 years ago

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