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

Last change on this file since 869 was 819, checked in by garnier, 16 years ago

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