Geant4  10.00.p01
G4VEmProcess.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: G4VEmProcess.cc 79268 2014-02-20 16:46:31Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4VEmProcess
34 //
35 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
36 //
37 // Creation date: 01.10.2003
38 //
39 // Modifications:
40 // 30-06-04 make it to be pure discrete process (V.Ivanchenko)
41 // 30-09-08 optimise integral option (V.Ivanchenko)
42 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
43 // 11-03-05 Shift verbose level by 1, add applyCuts and killPrimary flags (VI)
44 // 14-03-05 Update logic PostStepDoIt (V.Ivanchenko)
45 // 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
46 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivanchenko)
47 // 25-07-05 Add protection: integral mode only for charged particles (VI)
48 // 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko)
49 // 11-01-06 add A to parameters of ComputeCrossSectionPerAtom (VI)
50 // 12-09-06 add SetModel() (mma)
51 // 12-04-07 remove double call to Clear model manager (V.Ivanchenko)
52 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
53 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
54 // 17-02-10 Added pointer currentParticle (VI)
55 // 30-05-12 allow Russian roulette, brem splitting (D. Sawkey)
56 //
57 // Class Description:
58 //
59 
60 // -------------------------------------------------------------------
61 //
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 
65 #include "G4VEmProcess.hh"
66 #include "G4PhysicalConstants.hh"
67 #include "G4SystemOfUnits.hh"
68 #include "G4ProcessManager.hh"
69 #include "G4LossTableManager.hh"
70 #include "G4LossTableBuilder.hh"
71 #include "G4Step.hh"
72 #include "G4ParticleDefinition.hh"
73 #include "G4VEmModel.hh"
74 #include "G4DataVector.hh"
75 #include "G4PhysicsTable.hh"
76 #include "G4PhysicsLogVector.hh"
77 #include "G4VParticleChange.hh"
78 #include "G4ProductionCutsTable.hh"
79 #include "G4Region.hh"
80 #include "G4Gamma.hh"
81 #include "G4Electron.hh"
82 #include "G4Positron.hh"
83 #include "G4PhysicsTableHelper.hh"
84 #include "G4EmBiasingManager.hh"
85 #include "G4GenericIon.hh"
86 #include "G4Log.hh"
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
91  G4VDiscreteProcess(name, type),
92  secondaryParticle(0),
93  buildLambdaTable(true),
94  numberOfModels(0),
95  theLambdaTable(0),
96  theLambdaTablePrim(0),
97  theDensityFactor(0),
98  theDensityIdx(0),
99  integral(false),
100  applyCuts(false),
101  startFromNull(false),
102  splineFlag(true),
103  currentModel(0),
104  particle(0),
105  currentParticle(0),
106  currentCouple(0)
107 {
108  SetVerboseLevel(1);
109 
110  // Size of tables assuming spline
111  minKinEnergy = 0.1*keV;
112  maxKinEnergy = 10.0*TeV;
113  nLambdaBins = 77;
115 
116  // default lambda factor
117  lambdaFactor = 0.8;
118 
119  // default limit on polar angle
120  polarAngleLimit = 0.0;
121  biasFactor = 1.0;
122 
123  // particle types
127 
130  secParticles.reserve(5);
131 
132  preStepLambda = 0.0;
134 
135  idxLambda = idxLambdaPrim = 0;
136 
138  biasManager = 0;
139  biasFlag = false;
140  weightFlag = false;
142  lManager->Register(this);
143  secID = fluoID = augerID = biasID = -1;
144  mainSecondaries = 100;
145  if("phot" == GetProcessName() || "compt" == GetProcessName()) {
146  mainSecondaries = 1;
147  }
148 }
149 
150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 
153 {
154  if(1 < verboseLevel) {
155  G4cout << "G4VEmProcess destruct " << GetProcessName()
156  << " " << this << " " << theLambdaTable <<G4endl;
157  }
158  if(lManager->IsMaster()) {
159  delete theLambdaTable;
160  delete theLambdaTablePrim;
161  }
162  delete modelManager;
163  delete biasManager;
164  lManager->DeRegister(this);
165  //G4cout << "G4VEmProcess removed " << G4endl;
166 }
167 
168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
169 
171 {
172  currentCouple = 0;
173  preStepLambda = 0.0;
175  idxLambda = idxLambdaPrim = 0;
176 }
177 
178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179 
181  const G4Material*)
182 {
183  return 0.0;
184 }
185 
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187 
189  const G4Region* region)
190 {
192  modelManager->AddEmModel(order, p, fm, region);
193  if(p) { p->SetParticleChange(pParticleChange); }
194 }
195 
196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
197 
199 {
200  G4int n = emModels.size();
201  if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
202  emModels[index] = p;
203 }
204 
205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206 
208 {
209  G4VEmModel* p = 0;
210  if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; }
211  return p;
212 }
213 
214 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
215 
217  G4double emin, G4double emax)
218 {
219  modelManager->UpdateEmModel(nam, emin, emax);
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
223 
225 {
226  return modelManager->GetModel(idx, ver);
227 }
228 
229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
230 
232 {
233  G4bool isMaster = false;
234  if(GetMasterProcess() == this) { isMaster = true; }
235  if(!particle) { SetParticle(&part); }
236 
237  if(part.GetParticleType() == "nucleus" &&
238  part.GetParticleSubType() == "generic") {
239 
240  G4String pname = part.GetParticleName();
241  if(pname != "deuteron" && pname != "triton" &&
242  pname != "alpha" && pname != "He3" &&
243  pname != "alpha+" && pname != "helium" &&
244  pname != "hydrogen") {
245 
247  }
248  }
249 
250  if(1 < verboseLevel) {
251  G4cout << "G4VEmProcess::PreparePhysicsTable() for "
252  << GetProcessName()
253  << " and particle " << part.GetParticleName()
254  << " local particle " << particle->GetParticleName()
255  << G4endl;
256  }
257 
258  if(particle != &part) { return; }
259 
261 
262  lManager->PreparePhysicsTable(&part, this, isMaster);
263 
264  Clear();
266 
267  const G4ProductionCutsTable* theCoupleTable=
269  size_t n = theCoupleTable->GetTableSize();
270 
271  theEnergyOfCrossSectionMax.resize(n, 0.0);
272  theCrossSectionMax.resize(n, DBL_MAX);
273 
274  // initialisation of models
276  for(G4int i=0; i<numberOfModels; ++i) {
277  G4VEmModel* mod = modelManager->GetModel(i);
278  if(0 == i) { currentModel = mod; }
280  mod->SetMasterThread(isMaster);
281  if(mod->HighEnergyLimit() > maxKinEnergy) {
283  }
284  }
285 
288  2.,verboseLevel);
292 
293  // prepare tables
294  if(buildLambdaTable && isMaster){
295  theLambdaTable =
298  }
299  // high energy table
300  if(isMaster && minKinEnergyPrim < maxKinEnergy){
304  }
305  // forced biasing
306  if(biasManager) {
308  biasFlag = false;
309  }
310  // defined ID of secondary particles
311  if(isMaster) {
312  G4String nam1 = GetProcessName();
313  G4String nam2 = nam1 + "_fluo" ;
314  G4String nam3 = nam1 + "_auger";
315  G4String nam4 = nam1 + "_split";
320  }
321 }
322 
323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
324 
326 {
327  const G4VEmProcess* masterProc = 0;
328  if(GetMasterProcess() != this) {
329  masterProc = static_cast<const G4VEmProcess*>(GetMasterProcess());
330  }
331 
332  G4String num = part.GetParticleName();
333  if(1 < verboseLevel) {
334  G4cout << "G4VEmProcess::BuildPhysicsTable() for "
335  << GetProcessName()
336  << " and particle " << num
337  << " buildLambdaTable= " << buildLambdaTable
338  << G4endl;
339  }
340 
341  if(particle == &part) {
342 
344 
345  // worker initialisation
346  if(masterProc) {
347  theLambdaTable = masterProc->LambdaTable();
348  theLambdaTablePrim = masterProc->LambdaTablePrim();
349 
350  if(theLambdaTable) {
352  } else if(theLambdaTablePrim) {
354  }
357  if(theLambdaTable) { FindLambdaMax(); }
358 
359  // local initialisation of models
360  G4bool printing = true;
362  for(G4int i=0; i<numberOfModels; ++i) {
363  G4VEmModel* mod = GetModelByIndex(i, printing);
364  G4VEmModel* mod0= masterProc->GetModelByIndex(i, printing);
365  mod->InitialiseLocal(particle, mod0);
366  }
367  // master thread
368  } else {
373  }
374  }
375  }
376 
377  // explicitly defined printout by particle name
378  if(1 < verboseLevel ||
379  (0 < verboseLevel && (num == "gamma" || num == "e-" ||
380  num == "e+" || num == "mu+" ||
381  num == "mu-" || num == "proton"||
382  num == "pi+" || num == "pi-" ||
383  num == "kaon+" || num == "kaon-" ||
384  num == "alpha" || num == "anti_proton" ||
385  num == "GenericIon")))
386  {
387  PrintInfoProcess(part);
388  }
389 
390  if(1 < verboseLevel) {
391  G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
392  << GetProcessName()
393  << " and particle " << num
394  << G4endl;
395  }
396 }
397 
398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
399 
401 {
402  if(1 < verboseLevel) {
403  G4cout << "G4EmProcess::BuildLambdaTable() for process "
404  << GetProcessName() << " and particle "
405  << particle->GetParticleName() << " " << this
406  << G4endl;
407  }
408 
409  // Access to materials
410  const G4ProductionCutsTable* theCoupleTable=
412  size_t numOfCouples = theCoupleTable->GetTableSize();
413 
415 
416  G4PhysicsLogVector* aVector = 0;
417  G4PhysicsLogVector* aVectorPrim = 0;
418  G4PhysicsLogVector* bVectorPrim = 0;
419 
421  G4double emax1 = maxKinEnergy;
423 
424  for(size_t i=0; i<numOfCouples; ++i) {
425 
426  if (bld->GetFlag(i)) {
427 
428  // create physics vector and fill it
429  const G4MaterialCutsCouple* couple =
430  theCoupleTable->GetMaterialCutsCouple(i);
431 
432  // build main table
433  if(buildLambdaTable) {
434  delete (*theLambdaTable)[i];
435 
436  // if start from zero then change the scale
437  G4double emin = minKinEnergy;
438  G4bool startNull = false;
439  if(startFromNull) {
441  if(e >= emin) {
442  emin = e;
443  startNull = true;
444  }
445  }
446  G4double emax = emax1;
447  if(emax <= emin) { emax = 2*emin; }
448  G4int bin = G4lrint(nLambdaBins*G4Log(emax/emin)/scale);
449  if(bin < 3) { bin = 3; }
450  aVector = new G4PhysicsLogVector(emin, emax, bin);
451  aVector->SetSpline(splineFlag);
452  modelManager->FillLambdaVector(aVector, couple, startNull);
453  if(splineFlag) { aVector->FillSecondDerivatives(); }
455  }
456  // build high energy table
458  delete (*theLambdaTablePrim)[i];
459 
460  // start not from zero
461  if(!bVectorPrim) {
462  G4int bin =
464  if(bin < 3) { bin = 3; }
465  aVectorPrim =
467  bVectorPrim = aVectorPrim;
468  } else {
469  aVectorPrim = new G4PhysicsLogVector(*bVectorPrim);
470  }
471  // always use spline
472  aVectorPrim->SetSpline(true);
473  modelManager->FillLambdaVector(aVectorPrim, couple, false,
475  aVectorPrim->FillSecondDerivatives();
477  aVectorPrim);
478  }
479  }
480  }
481 
483 
484  if(1 < verboseLevel) {
485  G4cout << "Lambda table is built for "
487  << G4endl;
488  }
489 }
490 
491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
492 
494 {
495  if(verboseLevel > 0) {
496  G4cout << std::setprecision(6);
497  G4cout << G4endl << GetProcessName() << ": for "
498  << part.GetParticleName();
499  if(integral) { G4cout << ", integral: 1 "; }
500  if(applyCuts) { G4cout << ", applyCuts: 1 "; }
501  G4cout << " SubType= " << GetProcessSubType();;
502  if(biasFactor != 1.0) { G4cout << " BiasingFactor= " << biasFactor; }
503  G4cout << G4endl;
504  if(buildLambdaTable) {
505  size_t length = theLambdaTable->length();
506  for(size_t i=0; i<length; ++i) {
507  G4PhysicsVector* v = (*theLambdaTable)[i];
508  if(v) {
509  G4cout << " Lambda table from "
510  << G4BestUnit(minKinEnergy,"Energy")
511  << " to "
512  << G4BestUnit(v->GetMaxEnergy(),"Energy")
513  << " in " << v->GetVectorLength()-1
514  << " bins, spline: "
515  << splineFlag
516  << G4endl;
517  break;
518  }
519  }
520  }
522  size_t length = theLambdaTablePrim->length();
523  for(size_t i=0; i<length; ++i) {
524  G4PhysicsVector* v = (*theLambdaTablePrim)[i];
525  if(v) {
526  G4cout << " LambdaPrime table from "
527  << G4BestUnit(minKinEnergyPrim,"Energy")
528  << " to "
529  << G4BestUnit(maxKinEnergy,"Energy")
530  << " in " << v->GetVectorLength()-1
531  << " bins "
532  << G4endl;
533  break;
534  }
535  }
536  }
537  PrintInfo();
539  }
540 
541  if(verboseLevel > 2 && buildLambdaTable) {
542  G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
543  if(theLambdaTable) { G4cout << (*theLambdaTable) << G4endl; }
544  }
545 }
546 
547 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
548 
550 {
551  // reset parameters for the new track
554  //currentInteractionLength = -1.0;
555  // theInitialNumberOfInteractionLength=-1.0;
557 
558  // forced biasing only for primary particles
559  if(biasManager) {
560  if(0 == track->GetParentID()) {
561  // primary particle
562  biasFlag = true;
564  }
565  }
566 }
567 
568 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
569 
571  const G4Track& track,
572  G4double previousStepSize,
574 {
575  *condition = NotForced;
576  G4double x = DBL_MAX;
577 
581 
584  return x;
585  }
586 
587  // forced biasing only for primary particles
588  if(biasManager) {
589  if(0 == track.GetParentID()) {
590  if(biasFlag &&
592  return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
593  }
594  }
595  }
596 
597  // compute mean free path
601 
602  // zero cross section
603  if(preStepLambda <= 0.0) {
606  }
607  }
608 
609  // non-zero cross section
610  if(preStepLambda > 0.0) {
611 
613 
614  // beggining of tracking (or just after DoIt of this process)
615  // ResetNumberOfInteractionLengthLeft();
616 
619 
620  } else if(currentInteractionLength < DBL_MAX) {
621 
622  // subtract NumberOfInteractionLengthLeft using previous step
624  previousStepSize/currentInteractionLength;
625  //SubtractNumberOfInteractionLengthLeft(previousStepSize);
628  }
629  }
630 
631  // new mean free path and step limit for the next step
634 #ifdef G4VERBOSE
635  if (verboseLevel>2){
636  G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
637  G4cout << "[ " << GetProcessName() << "]" << G4endl;
638  G4cout << " for " << currentParticle->GetParticleName()
639  << " in Material " << currentMaterial->GetName()
640  << " Ekin(MeV)= " << preStepKinEnergy/MeV
641  <<G4endl;
642  G4cout << " MeanFreePath = " << currentInteractionLength/cm << "[cm]"
643  << " InteractionLength= " << x/cm <<"[cm] " <<G4endl;
644  }
645 #endif
646  }
647  return x;
648 }
649 
650 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
651 
653  const G4Step& step)
654 {
655  // In all cases clear number of interaction lengths
658 
660 
661  // Do not make anything if particle is stopped, the annihilation then
662  // should be performed by the AtRestDoIt!
663  if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
664 
665  G4double finalT = track.GetKineticEnergy();
666 
667  // forced process - should happen only once per track
668  if(biasFlag) {
670  biasFlag = false;
671  }
672  }
673 
674  // Integral approach
675  if (integral) {
676  G4double lx = GetLambda(finalT, currentCouple);
677  if(preStepLambda<lx && 1 < verboseLevel) {
678  G4cout << "WARNING: for " << currentParticle->GetParticleName()
679  << " and " << GetProcessName()
680  << " E(MeV)= " << finalT/MeV
681  << " preLambda= " << preStepLambda << " < "
682  << lx << " (postLambda) "
683  << G4endl;
684  }
685 
686  if(preStepLambda*G4UniformRand() > lx) {
688  return &fParticleChange;
689  }
690  }
691 
693  if(!currentModel->IsActive(finalT)) { return &fParticleChange; }
694 
695  // define new weight for primary and secondaries
697  if(weightFlag) {
698  weight /= biasFactor;
700  }
701 
702  /*
703  if(0 < verboseLevel) {
704  G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
705  << finalT/MeV
706  << " MeV; model= (" << currentModel->LowEnergyLimit()
707  << ", " << currentModel->HighEnergyLimit() << ")"
708  << G4endl;
709  }
710  */
711 
712  // sample secondaries
713  secParticles.clear();
715  currentCouple,
716  track.GetDynamicParticle(),
717  (*theCuts)[currentCoupleIndex]);
718 
719  G4int num0 = secParticles.size();
720 
721  // splitting or Russian roulette
722  if(biasManager) {
724  G4double eloss = 0.0;
726  secParticles, track, currentModel, &fParticleChange, eloss,
728  step.GetPostStepPoint()->GetSafety());
729  if(eloss > 0.0) {
732  }
733  }
734  }
735 
736  // save secondaries
737  G4int num = secParticles.size();
738  if(num > 0) {
739 
742  G4double time = track.GetGlobalTime();
743 
744  for (G4int i=0; i<num; ++i) {
745  if (secParticles[i]) {
748  G4double e = dp->GetKineticEnergy();
749  G4bool good = true;
750  if(applyCuts) {
751  if (p == theGamma) {
752  if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; }
753 
754  } else if (p == theElectron) {
755  if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; }
756 
757  } else if (p == thePositron) {
758  if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
760  good = false;
761  e += 2.0*electron_mass_c2;
762  }
763  }
764  // added secondary if it is good
765  }
766  if (good) {
767  G4Track* t = new G4Track(dp, time, track.GetPosition());
769  t->SetWeight(weight);
771 
772  // define type of secondary
774  else if(i < num0) {
775  if(p == theGamma) {
777  } else {
779  }
780  } else {
782  }
783 
784  //G4cout << "Secondary(post step) has weight " << t->GetWeight()
785  // << ", Ekin= " << t->GetKineticEnergy()/MeV << " MeV" <<G4endl;
786  } else {
787  delete dp;
788  edep += e;
789  }
790  }
791  }
793  }
794 
800  }
801 
802  // ClearNumberOfInteractionLengthLeft();
803  return &fParticleChange;
804 }
805 
806 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
807 
809  const G4String& directory,
810  G4bool ascii)
811 {
812  G4bool yes = true;
813 
814  if ( theLambdaTable && part == particle) {
815  const G4String name =
816  GetPhysicsTableFileName(part,directory,"Lambda",ascii);
817  yes = theLambdaTable->StorePhysicsTable(name,ascii);
818 
819  if ( yes ) {
820  G4cout << "Physics table is stored for " << particle->GetParticleName()
821  << " and process " << GetProcessName()
822  << " in the directory <" << directory
823  << "> " << G4endl;
824  } else {
825  G4cout << "Fail to store Physics Table for "
827  << " and process " << GetProcessName()
828  << " in the directory <" << directory
829  << "> " << G4endl;
830  }
831  }
832  if ( theLambdaTablePrim && part == particle) {
833  const G4String name =
834  GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
835  yes = theLambdaTablePrim->StorePhysicsTable(name,ascii);
836 
837  if ( yes ) {
838  G4cout << "Physics table prim is stored for "
840  << " and process " << GetProcessName()
841  << " in the directory <" << directory
842  << "> " << G4endl;
843  } else {
844  G4cout << "Fail to store Physics Table Prim for "
846  << " and process " << GetProcessName()
847  << " in the directory <" << directory
848  << "> " << G4endl;
849  }
850  }
851  return yes;
852 }
853 
854 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
855 
857  const G4String& directory,
858  G4bool ascii)
859 {
860  if(1 < verboseLevel) {
861  G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
862  << part->GetParticleName() << " and process "
863  << GetProcessName() << G4endl;
864  }
865  G4bool yes = true;
866 
868  || particle != part) { return yes; }
869 
870  const G4String particleName = part->GetParticleName();
871  G4String filename;
872 
873  if(buildLambdaTable) {
874  filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
876  filename,ascii);
877  if ( yes ) {
878  if (0 < verboseLevel) {
879  G4cout << "Lambda table for " << particleName
880  << " is Retrieved from <"
881  << filename << ">"
882  << G4endl;
883  }
884  if(lManager->SplineFlag()) {
885  size_t n = theLambdaTable->length();
886  for(size_t i=0; i<n; ++i) {
887  if((* theLambdaTable)[i]) {
888  (* theLambdaTable)[i]->SetSpline(true);
889  }
890  }
891  }
892  } else {
893  if (1 < verboseLevel) {
894  G4cout << "Lambda table for " << particleName << " in file <"
895  << filename << "> is not exist"
896  << G4endl;
897  }
898  }
899  }
901  filename = GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
903  filename,ascii);
904  if ( yes ) {
905  if (0 < verboseLevel) {
906  G4cout << "Lambda table prim for " << particleName
907  << " is Retrieved from <"
908  << filename << ">"
909  << G4endl;
910  }
911  if(lManager->SplineFlag()) {
912  size_t n = theLambdaTablePrim->length();
913  for(size_t i=0; i<n; ++i) {
914  if((* theLambdaTablePrim)[i]) {
915  (* theLambdaTablePrim)[i]->SetSpline(true);
916  }
917  }
918  }
919  } else {
920  if (1 < verboseLevel) {
921  G4cout << "Lambda table prim for " << particleName << " in file <"
922  << filename << "> is not exist"
923  << G4endl;
924  }
925  }
926  }
927 
928  return yes;
929 }
930 
931 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
932 
933 G4double
935  const G4MaterialCutsCouple* couple)
936 {
937  // Cross section per atom is calculated
938  DefineMaterial(couple);
939  G4double cross = 0.0;
941  cross = GetCurrentLambda(kineticEnergy);
942  } else {
943  SelectModel(kineticEnergy, currentCoupleIndex);
946  kineticEnergy);
947  }
948 
949  if(cross < 0.0) { cross = 0.0; }
950  return cross;
951 }
952 
953 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
954 
956  G4double,
958 {
959  *condition = NotForced;
960  return G4VEmProcess::MeanFreePath(track);
961 }
962 
963 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
964 
966 {
969  G4double x = DBL_MAX;
970  if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
971  return x;
972 }
973 
974 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
975 
976 G4double
978  G4double Z, G4double A, G4double cut)
979 {
980  SelectModel(kineticEnergy, currentCoupleIndex);
981  G4double x = 0.0;
982  if(currentModel) {
984  Z,A,cut);
985  }
986  return x;
987 }
988 
989 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
990 
992 {
993  if(1 < verboseLevel) {
994  G4cout << "### G4VEmProcess::FindLambdaMax: "
996  << " and process " << GetProcessName() << " " << G4endl;
997  }
998  size_t n = theLambdaTable->length();
999  G4PhysicsVector* pv;
1000  G4double e, ss, emax, smax;
1001 
1002  size_t i;
1003 
1004  // first loop on existing vectors
1005  for (i=0; i<n; ++i) {
1006  pv = (*theLambdaTable)[i];
1007  if(pv) {
1008  size_t nb = pv->GetVectorLength();
1009  emax = DBL_MAX;
1010  smax = 0.0;
1011  if(nb > 0) {
1012  for (size_t j=0; j<nb; ++j) {
1013  e = pv->Energy(j);
1014  ss = (*pv)(j);
1015  if(ss > smax) {
1016  smax = ss;
1017  emax = e;
1018  }
1019  }
1020  }
1021  theEnergyOfCrossSectionMax[i] = emax;
1022  theCrossSectionMax[i] = smax;
1023  if(1 < verboseLevel) {
1024  G4cout << "For " << particle->GetParticleName()
1025  << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
1026  << " lambda= " << smax << G4endl;
1027  }
1028  }
1029  }
1030  // second loop using base materials
1031  for (i=0; i<n; ++i) {
1032  pv = (*theLambdaTable)[i];
1033  if(!pv){
1034  G4int j = (*theDensityIdx)[i];
1036  theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
1037  }
1038  }
1039 }
1040 
1041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1042 
1045 {
1046  G4PhysicsVector* v =
1048  v->SetSpline(lManager->SplineFlag());
1049  return v;
1050 }
1051 
1052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1053 
1055 {
1056  const G4Element* elm = 0;
1058  return elm;
1059 }
1060 
1061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1062 
1064 {
1065  if(f > 0.0) {
1066  biasFactor = f;
1067  weightFlag = flag;
1068  if(1 < verboseLevel) {
1069  G4cout << "### SetCrossSectionBiasingFactor: for "
1070  << particle->GetParticleName()
1071  << " and process " << GetProcessName()
1072  << " biasFactor= " << f << " weightFlag= " << flag
1073  << G4endl;
1074  }
1075  }
1076 }
1077 
1078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1079 
1080 void
1082  G4bool flag)
1083 {
1084  if(!biasManager) { biasManager = new G4EmBiasingManager(); }
1085  if(1 < verboseLevel) {
1086  G4cout << "### ActivateForcedInteraction: for "
1087  << particle->GetParticleName()
1088  << " and process " << GetProcessName()
1089  << " length(mm)= " << length/mm
1090  << " in G4Region <" << r
1091  << "> weightFlag= " << flag
1092  << G4endl;
1093  }
1094  weightFlag = flag;
1096 }
1097 
1098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1099 
1100 void
1102  G4double factor,
1103  G4double energyLimit)
1104 {
1105  if (0.0 <= factor) {
1106 
1107  // Range cut can be applied only for e-
1108  if(0.0 == factor && secondaryParticle != G4Electron::Electron())
1109  { return; }
1110 
1111  if(!biasManager) { biasManager = new G4EmBiasingManager(); }
1112  biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
1113  if(1 < verboseLevel) {
1114  G4cout << "### ActivateSecondaryBiasing: for "
1115  << " process " << GetProcessName()
1116  << " factor= " << factor
1117  << " in G4Region <" << region
1118  << "> energyLimit(MeV)= " << energyLimit/MeV
1119  << G4endl;
1120  }
1121  }
1122 }
1123 
1124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1125 
1127 {
1129  /std::log(maxKinEnergy/minKinEnergy));
1130  minKinEnergy = e;
1131 }
1132 
1133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1134 
1136 {
1138  /std::log(maxKinEnergy/minKinEnergy));
1139  maxKinEnergy = e;
1140 }
1141 
1142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void ActivateForcedInteraction(G4double length=0.0, const G4String &r="", G4bool flag=true)
static const double cm
Definition: G4SIunits.hh:106
G4double condition(const G4ErrorSymMatrix &m)
G4PhysicsTable * LambdaTable() const
const std::vector< G4double > * theCutsGamma
const G4VProcess * GetMasterProcess() const
Definition: G4VProcess.hh:538
std::vector< G4double > theCrossSectionMax
void DefineMaterial(const G4MaterialCutsCouple *couple)
virtual void PrintInfo()=0
G4int GetParentID() const
G4bool SplineFlag() const
G4double GetLambda(G4double &kinEnergy, const G4MaterialCutsCouple *couple)
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:245
static const double MeV
Definition: G4SIunits.hh:193
G4PhysicsTable * LambdaTablePrim() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4VEmModel * SelectModel(G4double &kinEnergy, size_t index)
const std::vector< G4double > * GetDensityFactors()
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
G4bool weightFlag
static G4LossTableManager * Instance()
void ActivateForcedInteraction(G4double length=0.0, const G4String &r="")
G4double GetMaxEnergy() const
G4double ApplySecondaryBiasing(std::vector< G4DynamicParticle * > &, const G4Track &track, G4VEmModel *currentModel, G4ParticleChangeForGamma *pParticleChange, G4double &eloss, G4int coupleIdx, G4double tcut, G4double safety=0.0)
G4int verboseLevel
Definition: G4VProcess.hh:368
void UpdateEmModel(const G4String &, G4double, G4double)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
G4PhysicsTable * theLambdaTablePrim
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
G4VEmProcess(const G4String &name, G4ProcessType type=fElectromagnetic)
Definition: G4VEmProcess.cc:90
const G4DynamicParticle * GetDynamicParticle() const
G4VEmModel * currentModel
G4bool SecondaryBiasingRegion(G4int coupleIdx)
G4double ComputeCrossSectionPerAtom(G4double kineticEnergy, G4double Z, G4double A=0., G4double cut=0.0)
G4VEmModel * EmModel(G4int index=1) const
G4VEmModel * GetModel(G4int, G4bool ver=false)
void DeRegister(G4VEnergyLossProcess *p)
G4String name
Definition: TRTMaterials.hh:40
const G4String & GetName() const
Definition: G4Material.hh:176
const G4ThreeVector & GetPosition() const
G4ParticleChangeForGamma fParticleChange
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
Definition: G4VEmModel.cc:206
G4TrackStatus GetTrackStatus() const
G4bool ForcedInteractionRegion(G4int coupleIdx)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
const std::vector< G4double > * theCutsPositron
virtual void InitialiseProcess(const G4ParticleDefinition *)=0
virtual ~G4VEmProcess()
G4double theNumberOfInteractionLengthLeft
Definition: G4VProcess.hh:293
const G4MaterialCutsCouple * currentCouple
G4bool GetFlag(size_t idx) const
void SetTouchableHandle(const G4TouchableHandle &apValue)
void AddEmModel(G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
G4double GetParentWeight() const
G4double polarAngleLimit
G4double lambdaFactor
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void BuildPhysicsTable(const G4ParticleDefinition &)
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
size_t GetVectorLength() const
const G4String & GetParticleSubType() const
G4double fFactor
G4PhysicsTable * theLambdaTable
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:447
void FillLambdaVector(G4PhysicsVector *, const G4MaterialCutsCouple *, G4bool startFromNull=true, G4EmTableType t=fRestricted)
G4double minKinEnergyPrim
G4ProcessManager * GetProcessManager() const
int G4int
Definition: G4Types.hh:78
void FillSecondDerivatives()
void UpdateEmModel(const G4String &, G4double, G4double)
G4double GetCurrentLambda(G4double kinEnergy)
const G4String & GetPhysicsTableFileName(const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)
Definition: G4VProcess.cc:186
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
const G4int smax
const G4String & GetParticleName() const
G4bool applyCuts
void ActivateSecondaryBiasing(const G4String &region, G4double factor, G4double energyLimit)
G4LossTableBuilder * GetTableBuilder()
void SetWeight(G4double aValue)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetCreatorModelIndex(G4int idx)
Definition: G4Track.hh:245
void SetSpline(G4bool)
void Initialise(const G4ParticleDefinition &part, const G4String &procName, G4int verbose)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
G4double mfpKinEnergy
void SetCrossSectionBiasingFactor(G4double f, G4bool flag=true)
virtual G4double MinPrimaryEnergy(const G4ParticleDefinition *, const G4Material *)
G4double MeanFreePath(const G4Track &track)
const G4ParticleDefinition * currentParticle
void ProposeWeight(G4double finalWeight)
G4double minKinEnergy
void SetSecondaryWeightByProcess(G4bool)
void SetEmModel(G4VEmModel *, G4int index=1)
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4double GetKineticEnergy() const
G4bool splineFlag
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4bool RetrievePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii)
void SetParticleChange(G4VParticleChange *, G4VEmFluctuationModel *f=0)
Definition: G4VEmModel.cc:387
G4bool startFromNull
bool G4bool
Definition: G4Types.hh:79
const std::vector< G4int > * theDensityIdx
void SetParticle(const G4ParticleDefinition *p)
G4double GetStepLimit(G4int coupleIdx, G4double previousStep)
G4double currentInteractionLength
Definition: G4VProcess.hh:297
void PrintInfoProcess(const G4ParticleDefinition &)
std::vector< G4DynamicParticle * > secParticles
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
const G4ParticleDefinition * GetParticleDefinition() const
G4bool IsMaster() const
G4LossTableManager * lManager
const G4String & GetParticleType() const
void Register(G4VEnergyLossProcess *p)
const G4ParticleDefinition * secondaryParticle
Definition: G4Step.hh:76
std::vector< G4double > theEnergyOfCrossSectionMax
const G4ParticleDefinition * GetParticleDefinition() const
const G4int n
G4double GetGlobalTime() const
const G4ParticleDefinition * particle
void PreparePhysicsTable(const G4ParticleDefinition &)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double Energy(size_t index) const
const std::vector< G4double > * theCuts
static const G4double A[nN]
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
void AddSecondary(G4Track *aSecondary)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:300
const G4TouchableHandle & GetTouchableHandle() const
G4bool IsActive(G4double kinEnergy)
Definition: G4VEmModel.hh:711
void DumpModelList(G4int verb)
const G4DataVector * Initialise(const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
void SetMasterThread(G4bool val)
Definition: G4VEmModel.hh:669
size_t length() const
G4double G4Log(G4double x)
Definition: G4Log.hh:227
std::vector< G4VEmModel * > emModels
G4int size() const
void ActivateSecondaryBiasing(const G4String &region, G4double factor, G4double energyLimit)
static G4ProductionCutsTable * GetProductionCutsTable()
G4PhysicsVector * LambdaPhysicsVector(const G4MaterialCutsCouple *)
#define fm
static G4Positron * Positron()
Definition: G4Positron.cc:94
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:93
G4int NumberOfModels() const
const G4Material * currentMaterial
void SetMaxKinEnergy(G4double e)
size_t idxLambdaPrim
void FindLambdaMax()
const G4ParticleDefinition * theElectron
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void AddEmModel(G4int, G4VEmModel *, const G4Region *region=0)
int G4lrint(double ad)
Definition: templates.hh:163
G4double GetProposedKineticEnergy() const
void InitialiseBaseMaterials(G4PhysicsTable *table)
size_t idxLambda
void SetNumberOfSecondaries(G4int totSecondaries)
G4double GetLocalEnergyDeposit() const
G4StepPoint * GetPostStepPoint() const
G4EmModelManager * modelManager
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
const G4Element * GetCurrentElement() const
G4double preStepLambda
void ComputeIntegralLambda(G4double kinEnergy)
G4EmBiasingManager * biasManager
G4double GetSafety() const
void PreparePhysicsTable(const G4ParticleDefinition *aParticle, G4VEnergyLossProcess *p, G4bool theMaster)
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
void SetFluoFlag(G4bool val)
static const double TeV
Definition: G4SIunits.hh:197
G4bool StorePhysicsTable(const G4String &filename, G4bool ascii=false)
G4TrackStatus GetTrackStatus() const
G4VEmModel * GetModelByIndex(G4int idx=0, G4bool ver=false) const
G4int numberOfModels
const std::vector< G4double > * theDensityFactor
G4VAtomDeexcitation * AtomDeexcitation()
static const double keV
Definition: G4SIunits.hh:195
G4double biasFactor
size_t currentCoupleIndex
void SetMinKinEnergy(G4double e)
G4double theInitialNumberOfInteractionLength
Definition: G4VProcess.hh:300
static G4int Register(G4String &)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
const std::vector< G4double > * theCutsElectron
void InitializeForPostStep(const G4Track &)
G4ForceCondition
void StartTracking(G4Track *)
static G4bool RetrievePhysicsTable(G4PhysicsTable *physTable, const G4String &fileName, G4bool ascii)
#define DBL_MAX
Definition: templates.hh:83
G4double CrossSectionPerVolume(G4double kineticEnergy, const G4MaterialCutsCouple *couple)
G4bool StorePhysicsTable(const G4ParticleDefinition *, const G4String &directory, G4bool ascii=false)
static const double mm
Definition: G4SIunits.hh:102
G4int GetProcessSubType() const
Definition: G4VProcess.hh:426
G4double maxKinEnergy
G4double preStepKinEnergy
G4bool buildLambdaTable
void SetPolarAngleLimit(G4double)
Definition: G4VEmModel.hh:718
G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
const G4ParticleDefinition * theGamma
void SetVerboseLevel(G4int value)
Definition: G4VProcess.hh:437
const G4ParticleDefinition * thePositron
const G4Material * GetMaterial() const
G4int mainSecondaries
void BuildLambdaTable()
const std::vector< G4int > * GetCoupleIndexes()
G4ProcessType
const G4Element * GetCurrentElement() const
Definition: G4VEmModel.hh:440