Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LossTableManager.cc
Go to the documentation of this file.
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$
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4LossTableManager
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 03.01.2002
38 //
39 // Modifications:
40 //
41 // 20-01-03 Migrade to cut per region (V.Ivanchenko)
42 // 15-02-03 Lambda table can be scaled (V.Ivanchenko)
43 // 17-02-03 Fix problem of store/restore tables (V.Ivanchenko)
44 // 10-03-03 Add Ion registration (V.Ivanchenko)
45 // 25-03-03 Add deregistration (V.Ivanchenko)
46 // 02-04-03 Change messenger (V.Ivanchenko)
47 // 26-04-03 Fix retrieve tables (V.Ivanchenko)
48 // 13-05-03 Add calculation of precise range (V.Ivanchenko)
49 // 23-07-03 Add exchange with G4EnergyLossTables (V.Ivanchenko)
50 // 05-10-03 Add G4VEmProcesses registration and Verbose command (V.Ivanchenko)
51 // 17-10-03 Add SetParameters method (V.Ivanchenko)
52 // 23-10-03 Add control on inactive processes (V.Ivanchenko)
53 // 04-11-03 Add checks in RetrievePhysicsTable (V.Ivanchenko)
54 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossProcess (V.Ivanchenko)
55 // 14-01-04 Activate precise range calculation (V.Ivanchenko)
56 // 10-03-04 Fix a problem of Precise Range table (V.Ivanchenko)
57 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
58 // 13-01-04 Fix problem which takes place for inactivate eIoni (V.Ivanchenko)
59 // 25-01-04 Fix initialisation problem for ions (V.Ivanchenko)
60 // 11-03-05 Shift verbose level by 1 (V.Ivantchenko)
61 // 10-01-06 PreciseRange -> CSDARange (V.Ivantchenko)
62 // 20-01-06 Introduce G4EmTableType to remove repeating code (VI)
63 // 23-03-06 Set flag isIonisation (VI)
64 // 10-05-06 Add methods SetMscStepLimitation, FacRange and MscFlag (VI)
65 // 22-05-06 Add methods Set/Get bremsTh (VI)
66 // 05-06-06 Do not clear loss_table map between runs (VI)
67 // 16-01-07 Create new energy loss table for e+,e-,mu+,mu- and
68 // left ionisation table for further usage (VI)
69 // 12-02-07 Add SetSkin, SetLinearLossLimit (V.Ivanchenko)
70 // 18-06-07 Move definition of msc parameters to G4EmProcessOptions (V.Ivanchenko)
71 // 21-02-08 Added G4EmSaturation (V.Ivanchenko)
72 // 12-04-10 Added PreparePhsyicsTables and BuildPhysicsTables entries (V.Ivanchenko)
73 //
74 // Class Description:
75 //
76 // -------------------------------------------------------------------
77 //
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
81 #include "G4LossTableManager.hh"
82 #include "G4SystemOfUnits.hh"
83 #include "G4EnergyLossMessenger.hh"
84 #include "G4PhysicsTable.hh"
85 #include "G4ParticleDefinition.hh"
86 #include "G4MaterialCutsCouple.hh"
87 #include "G4ProcessManager.hh"
88 #include "G4Electron.hh"
89 #include "G4Proton.hh"
90 #include "G4VMultipleScattering.hh"
91 #include "G4VEmProcess.hh"
92 #include "G4ProductionCutsTable.hh"
93 #include "G4PhysicsTableHelper.hh"
94 #include "G4EmCorrections.hh"
95 #include "G4EmSaturation.hh"
96 #include "G4EmConfigurator.hh"
97 #include "G4ElectronIonPair.hh"
98 #include "G4EmTableType.hh"
99 #include "G4LossTableBuilder.hh"
100 #include "G4VAtomDeexcitation.hh"
101 #include "G4Region.hh"
102 
103 G4LossTableManager* G4LossTableManager::theInstance = 0;
104 
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
106 
108 {
109  if(0 == theInstance) {
110  static G4LossTableManager manager;
111  theInstance = &manager;
112  }
113  return theInstance;
114 }
115 
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
117 
119 {
120  for (G4int i=0; i<n_loss; ++i) {
121  if( loss_vector[i] ) { delete loss_vector[i]; }
122  }
123  size_t msc = msc_vector.size();
124  for (size_t j=0; j<msc; ++j) {
125  if( msc_vector[j] ) { delete msc_vector[j]; }
126  }
127  size_t emp = emp_vector.size();
128  for (size_t k=0; k<emp; ++k) {
129  if( emp_vector[k] ) { delete emp_vector[k]; }
130  }
131  size_t mod = mod_vector.size();
132  for (size_t a=0; a<mod; ++a) {
133  if( mod_vector[a] ) { delete mod_vector[a]; }
134  }
135  size_t fmod = fmod_vector.size();
136  for (size_t b=0; b<fmod; ++b) {
137  if( fmod_vector[b] ) { delete fmod_vector[b]; }
138  }
139  Clear();
140  delete theMessenger;
141  delete tableBuilder;
142  delete emCorrections;
143  delete emSaturation;
144  delete emConfigurator;
145  delete emElectronIonPair;
146  delete atomDeexcitation;
147 }
148 
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
150 
151 G4LossTableManager::G4LossTableManager()
152 {
153  n_loss = 0;
154  run = 0;
155  startInitialisation = false;
156  all_tables_are_built = false;
157  currentLoss = 0;
158  currentParticle = 0;
159  firstParticle = 0;
160  lossFluctuationFlag = true;
161  subCutoffFlag = false;
162  rndmStepFlag = false;
163  minSubRange = 0.0;
164  maxRangeVariation = 1.0;
165  maxFinalStep = 0.0;
166  minKinEnergy = 0.1*keV;
167  maxKinEnergy = 10.0*TeV;
168  nbinsLambda = 77;
169  nbinsPerDecade = 7;
170  maxKinEnergyForMuons = 10.*TeV;
171  integral = true;
172  integralActive = false;
173  buildCSDARange = false;
174  minEnergyActive = false;
175  maxEnergyActive = false;
176  maxEnergyForMuonsActive = false;
177  stepFunctionActive = false;
178  flagLPM = true;
179  splineFlag = true;
180  bremsTh = DBL_MAX;
181  factorForAngleLimit = 1.0;
182  verbose = 1;
183  theMessenger = new G4EnergyLossMessenger();
184  theElectron = G4Electron::Electron();
185  tableBuilder = new G4LossTableBuilder();
186  emCorrections= new G4EmCorrections();
187  emSaturation = new G4EmSaturation();
188  emConfigurator = new G4EmConfigurator(verbose);
189  emElectronIonPair = new G4ElectronIonPair();
190  tableBuilder->SetSplineFlag(splineFlag);
191  atomDeexcitation = 0;
192 }
193 
194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
195 
197 {
198  all_tables_are_built = false;
199  currentLoss = 0;
200  currentParticle = 0;
201  if(n_loss)
202  {
203  dedx_vector.clear();
204  range_vector.clear();
205  inv_range_vector.clear();
206  loss_map.clear();
207  loss_vector.clear();
208  part_vector.clear();
209  base_part_vector.clear();
210  tables_are_built.clear();
211  isActive.clear();
212  n_loss = 0;
213  }
214 }
215 
216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
217 
219 {
220  if(!p) { return; }
221  for (G4int i=0; i<n_loss; ++i) {
222  if(loss_vector[i] == p) { return; }
223  }
224  if(verbose > 1) {
225  G4cout << "G4LossTableManager::Register G4VEnergyLossProcess : "
226  << p->GetProcessName() << " idx= " << n_loss << G4endl;
227  }
228  ++n_loss;
229  loss_vector.push_back(p);
230  part_vector.push_back(0);
231  base_part_vector.push_back(0);
232  dedx_vector.push_back(0);
233  range_vector.push_back(0);
234  inv_range_vector.push_back(0);
235  tables_are_built.push_back(false);
236  isActive.push_back(true);
237  all_tables_are_built = false;
238  if(!lossFluctuationFlag) { p->SetLossFluctuations(false); }
239  if(subCutoffFlag) { p->ActivateSubCutoff(true); }
240  if(rndmStepFlag) { p->SetRandomStep(true); }
241  if(stepFunctionActive) { p->SetStepFunction(maxRangeVariation, maxFinalStep); }
242  if(integralActive) { p->SetIntegral(integral); }
243  if(minEnergyActive) { p->SetMinKinEnergy(minKinEnergy); }
244  if(maxEnergyActive) { p->SetMaxKinEnergy(maxKinEnergy); }
245 }
246 
247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
248 
250 {
251  if(!p) { return; }
252  for (G4int i=0; i<n_loss; ++i) {
253  if(loss_vector[i] == p) { loss_vector[i] = 0; }
254  }
255 }
256 
257 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
258 
260 {
261  if(!p) { return; }
262  G4int n = msc_vector.size();
263  for (G4int i=0; i<n; ++i) {
264  if(msc_vector[i] == p) { return; }
265  }
266  if(verbose > 1) {
267  G4cout << "G4LossTableManager::Register G4VMultipleScattering : "
268  << p->GetProcessName() << " idx= " << msc_vector.size() << G4endl;
269  }
270  msc_vector.push_back(p);
271 }
272 
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
274 
276 {
277  if(!p) { return; }
278  size_t msc = msc_vector.size();
279  for (size_t i=0; i<msc; ++i) {
280  if(msc_vector[i] == p) { msc_vector[i] = 0; }
281  }
282 }
283 
284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
285 
287 {
288  if(!p) { return; }
289  G4int n = emp_vector.size();
290  for (G4int i=0; i<n; ++i) {
291  if(emp_vector[i] == p) { return; }
292  }
293  if(verbose > 1) {
294  G4cout << "G4LossTableManager::Register G4VEmProcess : "
295  << p->GetProcessName() << " idx= " << emp_vector.size() << G4endl;
296  }
297  emp_vector.push_back(p);
298 }
299 
300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
301 
303 {
304  if(!p) { return; }
305  size_t emp = emp_vector.size();
306  for (size_t i=0; i<emp; ++i) {
307  if(emp_vector[i] == p) { emp_vector[i] = 0; }
308  }
309 }
310 
311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
312 
314 {
315  mod_vector.push_back(p);
316  if(verbose > 1) {
317  G4cout << "G4LossTableManager::Register G4VEmModel : "
318  << p->GetName() << G4endl;
319  }
320 }
321 
322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
323 
325 {
326  size_t n = mod_vector.size();
327  for (size_t i=0; i<n; ++i) {
328  if(mod_vector[i] == p) { mod_vector[i] = 0; }
329  }
330 }
331 
332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
333 
335 {
336  fmod_vector.push_back(p);
337  if(verbose > 1) {
338  G4cout << "G4LossTableManager::Register G4VEmFluctuationModel : "
339  << p->GetName() << G4endl;
340  }
341 }
342 
343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
344 
346 {
347  size_t n = fmod_vector.size();
348  for (size_t i=0; i<n; ++i) {
349  if(fmod_vector[i] == p) { fmod_vector[i] = 0; }
350  }
351 }
352 
353 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
354 
357 {
358  loss_map[ion] = p;
359 }
360 
361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
362 
364  const G4ParticleDefinition* part,
366 {
367  if(!p || !part) { return; }
368  for (G4int i=0; i<n_loss; ++i) {
369  if(loss_vector[i] == p) { return; }
370  }
371  if(verbose > 1) {
372  G4cout << "G4LossTableManager::RegisterExtraParticle "
373  << part->GetParticleName() << " G4VEnergyLossProcess : "
374  << p->GetProcessName() << " idx= " << n_loss << G4endl;
375  }
376  ++n_loss;
377  loss_vector.push_back(p);
378  part_vector.push_back(part);
379  base_part_vector.push_back(p->BaseParticle());
380  dedx_vector.push_back(0);
381  range_vector.push_back(0);
382  inv_range_vector.push_back(0);
383  tables_are_built.push_back(false);
384  all_tables_are_built = false;
385 }
386 
387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
388 
389 void
392 {
393  if (1 < verbose) {
394  G4cout << "G4LossTableManager::PreparePhysicsTable for "
395  << particle->GetParticleName()
396  << " and " << p->GetProcessName() << " run= " << run
397  << " loss_vector " << loss_vector.size() << G4endl;
398  }
399  if(!startInitialisation) { tableBuilder->SetInitialisationFlag(false); }
400 
401  // start initialisation for the first run
402  startInitialisation = true;
403 
404  if( 0 == run ) {
405  emConfigurator->PrepareModels(particle, p);
406 
407  // initialise particles for given process
408  for (G4int j=0; j<n_loss; ++j) {
409  if (p == loss_vector[j]) {
410  if (!part_vector[j]) { part_vector[j] = particle; }
411  }
412  }
413  }
414 }
415 
416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
417 
418 void
420  G4VEmProcess* p)
421 {
422  if (1 < verbose) {
423  G4cout << "G4LossTableManager::PreparePhysicsTable for "
424  << particle->GetParticleName()
425  << " and " << p->GetProcessName() << G4endl;
426  }
427  if(!startInitialisation) { tableBuilder->SetInitialisationFlag(false); }
428 
429  // start initialisation for the first run
430  if( 0 == run ) {
431  emConfigurator->PrepareModels(particle, p);
432  }
433  startInitialisation = true;
434 }
435 
436 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
437 
438 void
441 {
442  if (1 < verbose) {
443  G4cout << "G4LossTableManager::PreparePhysicsTable for "
444  << particle->GetParticleName()
445  << " and " << p->GetProcessName() << G4endl;
446  }
447  if(!startInitialisation) { tableBuilder->SetInitialisationFlag(false); }
448 
449  // start initialisation for the first run
450  if( 0 == run ) {
451  emConfigurator->PrepareModels(particle, p);
452  }
453 }
454 
455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
456 
457 void
459 {
460  if(0 == run && startInitialisation) {
461  emConfigurator->Clear();
462  }
463 }
464 
465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
466 
468  const G4ParticleDefinition* aParticle,
470 {
471  if(1 < verbose) {
472  G4cout << "### G4LossTableManager::BuildDEDXTable() is requested for "
473  << aParticle->GetParticleName()
474  << " and process " << p->GetProcessName() << G4endl;
475  }
476  // clear configurator
477  if(0 == run && startInitialisation) {
478  emConfigurator->Clear();
479  firstParticle = aParticle;
480  }
481  if(startInitialisation && atomDeexcitation) {
482  atomDeexcitation->InitialiseAtomicDeexcitation();
483  }
484  startInitialisation = false;
485 
486  // initialisation before any table is built
487  if ( aParticle == firstParticle ) {
488  all_tables_are_built = true;
489 
490  if(1 < verbose) {
491  G4cout << "### G4LossTableManager start initilisation for first particle "
492  << firstParticle->GetParticleName()
493  << G4endl;
494  }
495  for (G4int i=0; i<n_loss; ++i) {
496  G4VEnergyLossProcess* el = loss_vector[i];
497 
498  if(el) {
499  const G4ProcessManager* pm = el->GetProcessManager();
500  isActive[i] = false;
501  if(pm) { isActive[i] = pm->GetProcessActivation(el); }
502  if(0 == run) { base_part_vector[i] = el->BaseParticle(); }
503  tables_are_built[i] = false;
504  all_tables_are_built= false;
505  if(!isActive[i]) { el->SetIonisation(false); }
506 
507  if(1 < verbose) {
508  G4cout << i <<". "<< el->GetProcessName();
509  if(el->Particle()) {
510  G4cout << " for " << el->Particle()->GetParticleName();
511  }
512  G4cout << " active= " << isActive[i]
513  << " table= " << tables_are_built[i]
514  << " isIonisation= " << el->IsIonisationProcess();
515  if(base_part_vector[i]) {
516  G4cout << " base particle " << base_part_vector[i]->GetParticleName();
517  }
518  G4cout << G4endl;
519  }
520  } else {
521  tables_are_built[i] = true;
522  part_vector[i] = 0;
523  }
524  }
525  ++run;
526  currentParticle = 0;
527  }
528 
529  // Set run time parameters
530  SetParameters(aParticle, p);
531 
532  if (all_tables_are_built) { return; }
533 
534  // Build tables for given particle
535  all_tables_are_built = true;
536 
537  for(G4int i=0; i<n_loss; ++i) {
538  if(p == loss_vector[i] && !tables_are_built[i] && !base_part_vector[i]) {
539  const G4ParticleDefinition* curr_part = part_vector[i];
540  if(1 < verbose) {
541  G4cout << "### BuildPhysicsTable for " << p->GetProcessName()
542  << " and " << curr_part->GetParticleName()
543  << " start BuildTable " << G4endl;
544  }
545  G4VEnergyLossProcess* curr_proc = BuildTables(curr_part);
546  if(curr_proc) { CopyTables(curr_part, curr_proc); }
547  }
548  if ( !tables_are_built[i] ) { all_tables_are_built = false; }
549  }
550 
551  if(1 < verbose) {
552  G4cout << "### G4LossTableManager::BuildDEDXTable end: "
553  << "all_tables_are_built= " << all_tables_are_built
554  << G4endl;
555 
556  if(all_tables_are_built) {
557  G4cout << "### All dEdx and Range tables are built #####" << G4endl;
558  }
559  }
560 }
561 
562 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
563 
564 void G4LossTableManager::CopyTables(const G4ParticleDefinition* part,
565  G4VEnergyLossProcess* base_proc)
566 {
567  for (G4int j=0; j<n_loss; ++j) {
568 
569  G4VEnergyLossProcess* proc = loss_vector[j];
570  //if(proc == base_proc || proc->Particle() == part)
571  // tables_are_built[j] = true;
572 
573  if (!tables_are_built[j] && part == base_part_vector[j]) {
574  tables_are_built[j] = true;
575  proc->SetDEDXTable(base_proc->DEDXTable(),fRestricted);
576  proc->SetDEDXTable(base_proc->DEDXTableForSubsec(),fSubRestricted);
577  proc->SetDEDXTable(base_proc->DEDXunRestrictedTable(),fTotal);
578  proc->SetCSDARangeTable(base_proc->CSDARangeTable());
579  proc->SetRangeTableForLoss(base_proc->RangeTableForLoss());
580  proc->SetInverseRangeTable(base_proc->InverseRangeTable());
581  proc->SetLambdaTable(base_proc->LambdaTable());
582  proc->SetSubLambdaTable(base_proc->SubLambdaTable());
583  proc->SetIonisation(base_proc->IsIonisationProcess());
584  loss_map[part_vector[j]] = proc;
585  if (1 < verbose) {
586  G4cout << "For " << proc->GetProcessName()
587  << " for " << part_vector[j]->GetParticleName()
588  << " base_part= " << part->GetParticleName()
589  << " tables are assigned "
590  << G4endl;
591  }
592  }
593 
594  if (theElectron == part && theElectron == proc->SecondaryParticle() ) {
595  proc->SetSecondaryRangeTable(base_proc->RangeTableForLoss());
596  }
597  }
598 }
599 
600 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
601 
602 G4VEnergyLossProcess* G4LossTableManager::BuildTables(
603  const G4ParticleDefinition* aParticle)
604 {
605  if(1 < verbose) {
606  G4cout << "G4LossTableManager::BuildTables() for "
607  << aParticle->GetParticleName() << G4endl;
608  }
609 
610  std::vector<G4PhysicsTable*> t_list;
611  std::vector<G4VEnergyLossProcess*> loss_list;
612  loss_list.clear();
613  G4VEnergyLossProcess* em = 0;
614  G4VEnergyLossProcess* p = 0;
615  G4int iem = 0;
616  G4PhysicsTable* dedx = 0;
617  G4int i;
618 
619  for (i=0; i<n_loss; ++i) {
620  p = loss_vector[i];
621  if (p && aParticle == part_vector[i] && !tables_are_built[i]) {
622  if ((p->IsIonisationProcess() && isActive[i]) ||
623  !em || (em && !isActive[iem]) ) {
624  em = p;
625  iem= i;
626  }
627  dedx = p->BuildDEDXTable(fRestricted);
628  // G4cout << "Build DEDX table for " << aParticle->GetParticleName()
629  // << " " << dedx << " " << dedx->length() << G4endl;
630  p->SetDEDXTable(dedx,fRestricted);
631  t_list.push_back(dedx);
632  loss_list.push_back(p);
633  tables_are_built[i] = true;
634  }
635  }
636 
637  G4int n_dedx = t_list.size();
638  if (0 == n_dedx || !em) {
639  G4cout << "G4LossTableManager WARNING: no DEDX processes for "
640  << aParticle->GetParticleName() << G4endl;
641  return 0;
642  }
643  G4int nSubRegions = em->NumberOfSubCutoffRegions();
644 
645  if (1 < verbose) {
646  G4cout << "G4LossTableManager::BuildTables() start to build range tables"
647  << " and the sum of " << n_dedx << " processes"
648  << " iem= " << iem << " em= " << em->GetProcessName()
649  << " buildCSDARange= " << buildCSDARange
650  << " nSubRegions= " << nSubRegions
651  << G4endl;
652  }
653 
654  dedx = em->IonisationTable();
655  if (1 < n_dedx) {
656  em->SetDEDXTable(dedx, fIsIonisation);
657  dedx = 0;
659  tableBuilder->BuildDEDXTable(dedx, t_list);
660  em->SetDEDXTable(dedx, fRestricted);
661  }
662  dedx_vector[iem] = dedx;
663 
664  G4PhysicsTable* range = em->RangeTableForLoss();
665  if(!range) range = G4PhysicsTableHelper::PreparePhysicsTable(range);
666  range_vector[iem] = range;
667 
668  G4PhysicsTable* invrange = em->InverseRangeTable();
669  if(!invrange) invrange = G4PhysicsTableHelper::PreparePhysicsTable(invrange);
670  inv_range_vector[iem] = invrange;
671 
672  G4bool flag = em->IsIonisationProcess();
673  tableBuilder->BuildRangeTable(dedx, range, flag);
674  tableBuilder->BuildInverseRangeTable(range, invrange, flag);
675 
676  // if(1<verbose) G4cout << *dedx << G4endl;
677 
678  em->SetRangeTableForLoss(range);
679  em->SetInverseRangeTable(invrange);
680 
681  // if(1<verbose) G4cout << *range << G4endl;
682 
683  std::vector<G4PhysicsTable*> listSub;
684  std::vector<G4PhysicsTable*> listCSDA;
685 
686  for (i=0; i<n_dedx; ++i) {
687  p = loss_list[i];
688  p->SetIonisation(false);
690  if (0 < nSubRegions) {
691  dedx = p->BuildDEDXTable(fSubRestricted);
692  p->SetDEDXTable(dedx,fSubRestricted);
693  listSub.push_back(dedx);
695  if(p != em) em->AddCollaborativeProcess(p);
696  }
697  if(buildCSDARange) {
698  dedx = p->BuildDEDXTable(fTotal);
699  p->SetDEDXTable(dedx,fTotal);
700  listCSDA.push_back(dedx);
701  }
702  }
703 
704  if (0 < nSubRegions) {
705  G4PhysicsTable* dedxSub = em->IonisationTableForSubsec();
706  if (1 < listSub.size()) {
707  em->SetDEDXTable(dedxSub, fIsSubIonisation);
708  dedxSub = 0;
710  tableBuilder->BuildDEDXTable(dedxSub, listSub);
711  em->SetDEDXTable(dedxSub, fSubRestricted);
712  }
713  }
714  if(buildCSDARange) {
715  G4PhysicsTable* dedxCSDA = em->DEDXunRestrictedTable();
716  if (1 < n_dedx) {
717  dedxCSDA = 0;
718  dedxCSDA = G4PhysicsTableHelper::PreparePhysicsTable(dedxCSDA);
719  tableBuilder->BuildDEDXTable(dedxCSDA, listCSDA);
720  em->SetDEDXTable(dedxCSDA,fTotal);
721  }
722  G4PhysicsTable* rCSDA = em->CSDARangeTable();
723  if(!rCSDA) rCSDA = G4PhysicsTableHelper::PreparePhysicsTable(rCSDA);
724  tableBuilder->BuildRangeTable(dedxCSDA, rCSDA, flag);
725  em->SetCSDARangeTable(rCSDA);
726  }
727 
728  em->SetIonisation(true);
729  loss_map[aParticle] = em;
730 
731  if (1 < verbose) {
732  G4cout << "G4LossTableManager::BuildTables: Tables are built for "
733  << aParticle->GetParticleName()
734  << "; ionisation process: " << em->GetProcessName()
735  << G4endl;
736  }
737  return em;
738 }
739 
740 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
741 
743 {
744  return theMessenger;
745 }
746 
747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
748 
749 void G4LossTableManager::ParticleHaveNoLoss(
750  const G4ParticleDefinition* aParticle)
751 {
753  ed << "Energy loss process not found for " << aParticle->GetParticleName()
754  << " !" << G4endl;
755  G4Exception("G4LossTableManager::ParticleHaveNoLoss", "em0001",
756  FatalException, ed);
757 
758 }
759 
760 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
761 
763 {
764  return buildCSDARange;
765 }
766 
767 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
768 
770 {
771  lossFluctuationFlag = val;
772  for(G4int i=0; i<n_loss; ++i) {
773  if(loss_vector[i]) { loss_vector[i]->SetLossFluctuations(val); }
774  }
775 }
776 
777 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
778 
780 {
781  subCutoffFlag = val;
782  for(G4int i=0; i<n_loss; ++i) {
783  if(loss_vector[i]) { loss_vector[i]->ActivateSubCutoff(val, r); }
784  }
785 }
786 
787 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
788 
790 {
791  integral = val;
792  integralActive = true;
793  for(G4int i=0; i<n_loss; ++i) {
794  if(loss_vector[i]) { loss_vector[i]->SetIntegral(val); }
795  }
796  size_t emp = emp_vector.size();
797  for (size_t k=0; k<emp; ++k) {
798  if(emp_vector[k]) { emp_vector[k]->SetIntegral(val); }
799  }
800 }
801 
802 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
803 
805 {
806  minSubRange = val;
807  for(G4int i=0; i<n_loss; ++i) {
808  if(loss_vector[i]) { loss_vector[i]->SetMinSubRange(val); }
809  }
810 }
811 
812 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
813 
815 {
816  rndmStepFlag = val;
817  for(G4int i=0; i<n_loss; ++i) {
818  if(loss_vector[i]) { loss_vector[i]->SetRandomStep(val); }
819  }
820 }
821 
822 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
823 
825 {
826  minEnergyActive = true;
827  minKinEnergy = val;
828  for(G4int i=0; i<n_loss; ++i) {
829  if(loss_vector[i]) { loss_vector[i]->SetMinKinEnergy(val); }
830  }
831  size_t emp = emp_vector.size();
832  for (size_t k=0; k<emp; ++k) {
833  if(emp_vector[k]) { emp_vector[k]->SetMinKinEnergy(val); }
834  }
835 }
836 
837 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
838 
840 {
841  maxEnergyActive = true;
842  maxKinEnergy = val;
843  for(G4int i=0; i<n_loss; ++i) {
844  if(loss_vector[i]) { loss_vector[i]->SetMaxKinEnergy(val); }
845  }
846  size_t emp = emp_vector.size();
847  for (size_t k=0; k<emp; ++k) {
848  if(emp_vector[k]) { emp_vector[k]->SetMaxKinEnergy(val); }
849  }
850 }
851 
852 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
853 
855 {
856  for(G4int i=0; i<n_loss; ++i) {
857  if(loss_vector[i]) { loss_vector[i]->SetMaxKinEnergyForCSDARange(val); }
858  }
859 }
860 
861 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
862 
864 {
865  maxEnergyForMuonsActive = true;
866  maxKinEnergyForMuons = val;
867 }
868 
869 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
870 
872 {
873  for(G4int i=0; i<n_loss; ++i) {
874  if(loss_vector[i]) { loss_vector[i]->SetDEDXBinning(val); }
875  }
876 }
877 
878 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
879 
881 {
882  for(G4int i=0; i<n_loss; ++i) {
883  if(loss_vector[i]) { loss_vector[i]->SetDEDXBinningForCSDARange(val); }
884  }
885 }
886 
887 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
888 
890 {
891  G4int n = val/G4int(std::log10(maxKinEnergy/minKinEnergy) + 0.5);
892  if(n < 5) {
893  G4cout << "G4LossTableManager::SetLambdaBinning WARNING "
894  << "too small number of bins " << val << " ignored"
895  << G4endl;
896  return;
897  }
898  nbinsLambda = val;
899  nbinsPerDecade = n;
900  size_t emp = emp_vector.size();
901  for (size_t k=0; k<emp; ++k) {
902  if(emp_vector[k]) { emp_vector[k]->SetLambdaBinning(val); }
903  }
904 }
905 
906 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
907 
909 {
910  return nbinsPerDecade;
911 }
912 
913 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
914 
916 {
917  verbose = val;
918  for(G4int i=0; i<n_loss; ++i) {
919  if(loss_vector[i]) { loss_vector[i]->SetVerboseLevel(val); }
920  }
921  size_t msc = msc_vector.size();
922  for (size_t j=0; j<msc; ++j) {
923  if(msc_vector[j]) { msc_vector[j]->SetVerboseLevel(val); }
924  }
925  size_t emp = emp_vector.size();
926  for (size_t k=0; k<emp; ++k) {
927  if(emp_vector[k]) { emp_vector[k]->SetVerboseLevel(val); }
928  }
929  emConfigurator->SetVerbose(val);
930  //tableBuilder->SetVerbose(val);
931  //emCorrections->SetVerbose(val);
932  emSaturation->SetVerbose(val);
933  emElectronIonPair->SetVerbose(val);
934  if(atomDeexcitation) { atomDeexcitation->SetVerboseLevel(val); }
935 }
936 
937 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
938 
940 {
941  stepFunctionActive = true;
942  maxRangeVariation = v1;
943  maxFinalStep = v2;
944  for(G4int i=0; i<n_loss; ++i) {
945  if(loss_vector[i]) { loss_vector[i]->SetStepFunction(v1, v2); }
946  }
947 }
948 
949 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
950 
952 {
953  for(G4int i=0; i<n_loss; ++i) {
954  if(loss_vector[i]) { loss_vector[i]->SetLinearLossLimit(val); }
955  }
956 }
957 
958 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
959 
961 {
962  buildCSDARange = val;
963 }
964 
965 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
966 
967 void
968 G4LossTableManager::SetParameters(const G4ParticleDefinition* aParticle,
970 {
971  if(stepFunctionActive) { p->SetStepFunction(maxRangeVariation, maxFinalStep); }
972  if(integralActive) { p->SetIntegral(integral); }
973  if(minEnergyActive) { p->SetMinKinEnergy(minKinEnergy); }
974  if(maxEnergyActive) { p->SetMaxKinEnergy(maxKinEnergy); }
975  p->SetVerboseLevel(verbose);
976  if(maxEnergyForMuonsActive) {
977  G4double dm = std::abs(aParticle->GetPDGMass() - 105.7*MeV);
978  if(dm < 5.*MeV) { p->SetMaxKinEnergy(maxKinEnergyForMuons); }
979  }
980 }
981 
982 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
983 
984 const std::vector<G4VEnergyLossProcess*>&
986 {
987  return loss_vector;
988 }
989 
990 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
991 
992 const std::vector<G4VEmProcess*>& G4LossTableManager::GetEmProcessVector()
993 {
994  return emp_vector;
995 }
996 
997 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
998 
999 const std::vector<G4VMultipleScattering*>&
1001 {
1002  return msc_vector;
1003 }
1004 
1005 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1006 
1008 {
1009  flagLPM = val;
1010 }
1011 
1012 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1013 
1015 {
1016  return flagLPM;
1017 }
1018 
1019 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1020 
1022 {
1023  splineFlag = val;
1024  tableBuilder->SetSplineFlag(splineFlag);
1025 }
1026 
1027 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1028 
1030 {
1031  return splineFlag;
1032 }
1033 
1034 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1035 
1037 {
1038  bremsTh = val;
1039 }
1040 
1041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1042 
1044 {
1045  return bremsTh;
1046 }
1047 
1048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1049 
1051 {
1052  if(val > 0.0) { factorForAngleLimit = val; }
1053 }
1054 
1055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1056 
1058 {
1059  return factorForAngleLimit;
1060 }
1061 
1062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1063 
1065 {
1066  return minKinEnergy;
1067 }
1068 
1069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1070 
1072 {
1073  return maxKinEnergy;
1074 }
1075 
1076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1077 
1079 {
1080  return emCorrections;
1081 }
1082 
1083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1084 
1086 {
1087  return emSaturation;
1088 }
1089 
1090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1091 
1093 {
1094  return emConfigurator;
1095 }
1096 
1097 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1098 
1100 {
1101  return emElectronIonPair;
1102 }
1103 
1104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1105 
1107 {
1108  return atomDeexcitation;
1109 }
1110 
1111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1112 
1114 {
1115  return tableBuilder;
1116 }
1117 
1118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1119 
1121 {
1122  atomDeexcitation = p;
1123 }
1124 
1125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1126 
1129 {
1130  if(aParticle != currentParticle) {
1131  currentParticle = aParticle;
1132  std::map<PD,G4VEnergyLossProcess*,std::less<PD> >::const_iterator pos;
1133  if ((pos = loss_map.find(aParticle)) != loss_map.end()) {
1134  currentLoss = (*pos).second;
1135  } else {
1136  currentLoss = 0;
1137  //ParticleHaveNoLoss(aParticle);
1138  }
1139  }
1140  return currentLoss;
1141 }
1142 
1143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1144 
1146  G4double kineticEnergy,
1147  const G4MaterialCutsCouple *couple)
1148 {
1149  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
1150  G4double x = 0.0;
1151  if(currentLoss) { x = currentLoss->GetDEDX(kineticEnergy, couple); }
1152  else { x = G4EnergyLossTables::GetDEDX(currentParticle,
1153  kineticEnergy,couple,false); }
1154  return x;
1155 }
1156 
1157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1158 
1160  G4double kineticEnergy,
1161  const G4MaterialCutsCouple *couple)
1162 {
1163  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
1164  G4double x = 0.0;
1165  if(currentLoss) { x = currentLoss->GetDEDXForSubsec(kineticEnergy, couple); }
1166  return x;
1167 }
1168 
1169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1170 
1172  G4double kineticEnergy,
1173  const G4MaterialCutsCouple *couple)
1174 {
1175  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
1176  G4double x = DBL_MAX;
1177  if(currentLoss) { x = currentLoss->GetCSDARange(kineticEnergy, couple); }
1178  return x;
1179 }
1180 
1181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1182 
1183 G4double
1185  G4double kineticEnergy,
1186  const G4MaterialCutsCouple *couple)
1187 {
1188  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
1189  G4double x;
1190  if(currentLoss) { x = currentLoss->GetRangeForLoss(kineticEnergy, couple); }
1191  else { x = G4EnergyLossTables::GetRange(currentParticle,kineticEnergy,
1192  couple,false); }
1193  return x;
1194 }
1195 
1196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
1197 
1199  G4double kineticEnergy,
1200  const G4MaterialCutsCouple *couple)
1201 {
1202  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
1203  G4double x;
1204  if(currentLoss) { x = currentLoss->GetRange(kineticEnergy, couple); }
1205  else { x = G4EnergyLossTables::GetRange(currentParticle,kineticEnergy,
1206  couple,false); }
1207  return x;
1208 }
1209 
1210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1211 
1213  G4double range,
1214  const G4MaterialCutsCouple *couple)
1215 {
1216  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
1217  G4double x;
1218  if(currentLoss) { x = currentLoss->GetKineticEnergy(range, couple); }
1219  else { x = G4EnergyLossTables::GetPreciseEnergyFromRange(currentParticle,range,
1220  couple,false); }
1221  return x;
1222 }
1223 
1224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1225 
1227  const G4DynamicParticle* dp,
1228  G4double& length)
1229 {
1230  const G4ParticleDefinition* aParticle = dp->GetParticleDefinition();
1231  if(aParticle != currentParticle) { GetEnergyLossProcess(aParticle); }
1232  G4double x = 0.0;
1233  if(currentLoss) { currentLoss->GetDEDXDispersion(couple, dp, length); }
1234  return x;
1235 }
1236 
1237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....