Geant4  10.00.p02
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 81864 2014-06-06 11:30:54Z 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 = fFactor = 1.0;
122 
123  // particle types
127 
129 
132  secParticles.reserve(5);
133 
135 
136  preStepLambda = 0.0;
138 
140 
142  biasManager = 0;
143  biasFlag = false;
144  weightFlag = false;
146  lManager->Register(this);
147  secID = fluoID = augerID = biasID = -1;
148  mainSecondaries = 100;
149  if("phot" == GetProcessName() || "compt" == GetProcessName()) {
150  mainSecondaries = 1;
151  }
152 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155 
157 {
158  if(1 < verboseLevel) {
159  G4cout << "G4VEmProcess destruct " << GetProcessName()
160  << " " << this << " " << theLambdaTable <<G4endl;
161  }
162  if(lManager->IsMaster()) {
163  delete theLambdaTable;
164  delete theLambdaTablePrim;
165  }
166  delete modelManager;
167  delete biasManager;
168  lManager->DeRegister(this);
169  //G4cout << "G4VEmProcess removed " << G4endl;
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173 
175 {
176  currentCouple = 0;
177  preStepLambda = 0.0;
179  idxLambda = idxLambdaPrim = 0;
180 }
181 
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183 
185  const G4Material*)
186 {
187  return 0.0;
188 }
189 
190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
191 
193  const G4Region* region)
194 {
196  modelManager->AddEmModel(order, p, fm, region);
197  if(p) { p->SetParticleChange(pParticleChange); }
198 }
199 
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
201 
203 {
204  G4int n = emModels.size();
205  if(index >= n) { for(G4int i=n; i<=index; ++i) {emModels.push_back(0);} }
206  emModels[index] = p;
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210 
212 {
213  G4VEmModel* p = 0;
214  if(index >= 0 && index < G4int(emModels.size())) { p = emModels[index]; }
215  return p;
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
219 
221  G4double emin, G4double emax)
222 {
223  modelManager->UpdateEmModel(nam, emin, emax);
224 }
225 
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227 
229 {
230  return modelManager->GetModel(idx, ver);
231 }
232 
233 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
234 
236 {
237  G4bool isMaster = true;
238  const G4VEmProcess* masterProcess =
239  static_cast<const G4VEmProcess*>(GetMasterProcess());
240  if(masterProcess && masterProcess != this) { isMaster = false; }
241 
242  if(!particle) { SetParticle(&part); }
243 
244  if(part.GetParticleType() == "nucleus" &&
245  part.GetParticleSubType() == "generic") {
246 
247  G4String pname = part.GetParticleName();
248  if(pname != "deuteron" && pname != "triton" &&
249  pname != "alpha" && pname != "He3" &&
250  pname != "alpha+" && pname != "helium" &&
251  pname != "hydrogen") {
252 
254  }
255  }
256 
257  if(1 < verboseLevel) {
258  G4cout << "G4VEmProcess::PreparePhysicsTable() for "
259  << GetProcessName()
260  << " and particle " << part.GetParticleName()
261  << " local particle " << particle->GetParticleName()
262  << G4endl;
263  }
264 
265  if(particle != &part) { return; }
266 
268 
269  lManager->PreparePhysicsTable(&part, this, isMaster);
270 
271  Clear();
273 
274  const G4ProductionCutsTable* theCoupleTable=
276  size_t n = theCoupleTable->GetTableSize();
277 
278  theEnergyOfCrossSectionMax.resize(n, 0.0);
279  theCrossSectionMax.resize(n, DBL_MAX);
280 
281  // initialisation of models
283  for(G4int i=0; i<numberOfModels; ++i) {
284  G4VEmModel* mod = modelManager->GetModel(i);
285  if(0 == i) { currentModel = mod; }
287  mod->SetMasterThread(isMaster);
288  if(mod->HighEnergyLimit() > maxKinEnergy) {
290  }
291  }
292 
295  2.,verboseLevel);
299 
300  // prepare tables
301  if(buildLambdaTable && isMaster){
302  theLambdaTable =
305  }
306  // high energy table
307  if(isMaster && minKinEnergyPrim < maxKinEnergy){
311  }
312  // forced biasing
313  if(biasManager) {
315  biasFlag = false;
316  }
317  // defined ID of secondary particles
318  if(isMaster) {
319  G4String nam1 = GetProcessName();
320  G4String nam2 = nam1 + "_fluo" ;
321  G4String nam3 = nam1 + "_auger";
322  G4String nam4 = nam1 + "_split";
327  }
328 }
329 
330 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
331 
333 {
334  G4bool isMaster = true;
335  const G4VEmProcess* masterProc =
336  static_cast<const G4VEmProcess*>(GetMasterProcess());
337  if(masterProc && masterProc != this) { isMaster = false; }
338 
339  G4String num = part.GetParticleName();
340  if(1 < verboseLevel) {
341  G4cout << "G4VEmProcess::BuildPhysicsTable() for "
342  << GetProcessName()
343  << " and particle " << num
344  << " buildLambdaTable= " << buildLambdaTable
345  << G4endl;
346  }
347 
348  if(particle == &part) {
349 
351 
352  // worker initialisation
353  if(!isMaster) {
354  theLambdaTable = masterProc->LambdaTable();
355  theLambdaTablePrim = masterProc->LambdaTablePrim();
356 
357  if(theLambdaTable) {
359  } else if(theLambdaTablePrim) {
361  }
364  if(theLambdaTable) { FindLambdaMax(); }
365 
366  // local initialisation of models
367  G4bool printing = true;
369  for(G4int i=0; i<numberOfModels; ++i) {
370  G4VEmModel* mod = GetModelByIndex(i, printing);
371  G4VEmModel* mod0= masterProc->GetModelByIndex(i, printing);
372  mod->InitialiseLocal(particle, mod0);
373  }
374  // master thread
375  } else {
380  }
381  }
382  }
383 
384  // explicitly defined printout by particle name
385  if(1 < verboseLevel ||
386  (0 < verboseLevel && (num == "gamma" || num == "e-" ||
387  num == "e+" || num == "mu+" ||
388  num == "mu-" || num == "proton"||
389  num == "pi+" || num == "pi-" ||
390  num == "kaon+" || num == "kaon-" ||
391  num == "alpha" || num == "anti_proton" ||
392  num == "GenericIon")))
393  {
394  PrintInfoProcess(part);
395  }
396 
397  if(1 < verboseLevel) {
398  G4cout << "G4VEmProcess::BuildPhysicsTable() done for "
399  << GetProcessName()
400  << " and particle " << num
401  << G4endl;
402  }
403 }
404 
405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406 
408 {
409  if(1 < verboseLevel) {
410  G4cout << "G4EmProcess::BuildLambdaTable() for process "
411  << GetProcessName() << " and particle "
412  << particle->GetParticleName() << " " << this
413  << G4endl;
414  }
415 
416  // Access to materials
417  const G4ProductionCutsTable* theCoupleTable=
419  size_t numOfCouples = theCoupleTable->GetTableSize();
420 
422 
423  G4PhysicsLogVector* aVector = 0;
424  G4PhysicsLogVector* aVectorPrim = 0;
425  G4PhysicsLogVector* bVectorPrim = 0;
426 
428  G4double emax1 = maxKinEnergy;
430 
431  for(size_t i=0; i<numOfCouples; ++i) {
432 
433  if (bld->GetFlag(i)) {
434 
435  // create physics vector and fill it
436  const G4MaterialCutsCouple* couple =
437  theCoupleTable->GetMaterialCutsCouple(i);
438 
439  // build main table
440  if(buildLambdaTable) {
441  delete (*theLambdaTable)[i];
442 
443  // if start from zero then change the scale
444  G4double emin = minKinEnergy;
445  G4bool startNull = false;
446  if(startFromNull) {
448  if(e >= emin) {
449  emin = e;
450  startNull = true;
451  }
452  }
453  G4double emax = emax1;
454  if(emax <= emin) { emax = 2*emin; }
455  G4int bin = G4lrint(nLambdaBins*G4Log(emax/emin)/scale);
456  if(bin < 3) { bin = 3; }
457  aVector = new G4PhysicsLogVector(emin, emax, bin);
458  aVector->SetSpline(splineFlag);
459  modelManager->FillLambdaVector(aVector, couple, startNull);
460  if(splineFlag) { aVector->FillSecondDerivatives(); }
462  }
463  // build high energy table
465  delete (*theLambdaTablePrim)[i];
466 
467  // start not from zero
468  if(!bVectorPrim) {
469  G4int bin =
471  if(bin < 3) { bin = 3; }
472  aVectorPrim =
474  bVectorPrim = aVectorPrim;
475  } else {
476  aVectorPrim = new G4PhysicsLogVector(*bVectorPrim);
477  }
478  // always use spline
479  aVectorPrim->SetSpline(true);
480  modelManager->FillLambdaVector(aVectorPrim, couple, false,
482  aVectorPrim->FillSecondDerivatives();
484  aVectorPrim);
485  }
486  }
487  }
488 
490 
491  if(1 < verboseLevel) {
492  G4cout << "Lambda table is built for "
494  << G4endl;
495  }
496 }
497 
498 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
499 
501 {
502  if(verboseLevel > 0) {
503  G4cout << std::setprecision(6);
504  G4cout << G4endl << GetProcessName() << ": for "
505  << part.GetParticleName();
506  if(integral) { G4cout << ", integral: 1 "; }
507  if(applyCuts) { G4cout << ", applyCuts: 1 "; }
508  G4cout << " SubType= " << GetProcessSubType();;
509  if(biasFactor != 1.0) { G4cout << " BiasingFactor= " << biasFactor; }
510  G4cout << G4endl;
511  if(buildLambdaTable) {
512  size_t length = theLambdaTable->length();
513  for(size_t i=0; i<length; ++i) {
514  G4PhysicsVector* v = (*theLambdaTable)[i];
515  if(v) {
516  G4cout << " Lambda table from "
517  << G4BestUnit(minKinEnergy,"Energy")
518  << " to "
519  << G4BestUnit(v->GetMaxEnergy(),"Energy")
520  << " in " << v->GetVectorLength()-1
521  << " bins, spline: "
522  << splineFlag
523  << G4endl;
524  break;
525  }
526  }
527  }
529  size_t length = theLambdaTablePrim->length();
530  for(size_t i=0; i<length; ++i) {
531  G4PhysicsVector* v = (*theLambdaTablePrim)[i];
532  if(v) {
533  G4cout << " LambdaPrime table from "
534  << G4BestUnit(minKinEnergyPrim,"Energy")
535  << " to "
536  << G4BestUnit(maxKinEnergy,"Energy")
537  << " in " << v->GetVectorLength()-1
538  << " bins "
539  << G4endl;
540  break;
541  }
542  }
543  }
544  PrintInfo();
546  }
547 
548  if(verboseLevel > 2 && buildLambdaTable) {
549  G4cout << " LambdaTable address= " << theLambdaTable << G4endl;
550  if(theLambdaTable) { G4cout << (*theLambdaTable) << G4endl; }
551  }
552 }
553 
554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
555 
557 {
558  // reset parameters for the new track
561  //currentInteractionLength = -1.0;
562  // theInitialNumberOfInteractionLength=-1.0;
564 
565  // forced biasing only for primary particles
566  if(biasManager) {
567  if(0 == track->GetParentID()) {
568  // primary particle
569  biasFlag = true;
571  }
572  }
573 }
574 
575 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
576 
578  const G4Track& track,
579  G4double previousStepSize,
581 {
582  *condition = NotForced;
583  G4double x = DBL_MAX;
584 
588 
591  return x;
592  }
593 
594  // forced biasing only for primary particles
595  if(biasManager) {
596  if(0 == track.GetParentID()) {
597  if(biasFlag &&
599  return biasManager->GetStepLimit(currentCoupleIndex, previousStepSize);
600  }
601  }
602  }
603 
604  // compute mean free path
608 
609  // zero cross section
610  if(preStepLambda <= 0.0) {
613  }
614  }
615 
616  // non-zero cross section
617  if(preStepLambda > 0.0) {
618 
620 
621  // beggining of tracking (or just after DoIt of this process)
622  // ResetNumberOfInteractionLengthLeft();
623 
626 
627  } else if(currentInteractionLength < DBL_MAX) {
628 
629  // subtract NumberOfInteractionLengthLeft using previous step
631  previousStepSize/currentInteractionLength;
632  //SubtractNumberOfInteractionLengthLeft(previousStepSize);
635  }
636  }
637 
638  // new mean free path and step limit for the next step
641 #ifdef G4VERBOSE
642  if (verboseLevel>2){
643  G4cout << "G4VEmProcess::PostStepGetPhysicalInteractionLength ";
644  G4cout << "[ " << GetProcessName() << "]" << G4endl;
645  G4cout << " for " << currentParticle->GetParticleName()
646  << " in Material " << currentMaterial->GetName()
647  << " Ekin(MeV)= " << preStepKinEnergy/MeV
648  <<G4endl;
649  G4cout << " MeanFreePath = " << currentInteractionLength/cm << "[cm]"
650  << " InteractionLength= " << x/cm <<"[cm] " <<G4endl;
651  }
652 #endif
653  }
654  return x;
655 }
656 
657 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
658 
660  const G4Step& step)
661 {
662  // In all cases clear number of interaction lengths
665 
667 
668  // Do not make anything if particle is stopped, the annihilation then
669  // should be performed by the AtRestDoIt!
670  if (track.GetTrackStatus() == fStopButAlive) { return &fParticleChange; }
671 
672  G4double finalT = track.GetKineticEnergy();
673 
674  // forced process - should happen only once per track
675  if(biasFlag) {
677  biasFlag = false;
678  }
679  }
680 
681  // Integral approach
682  if (integral) {
683  G4double lx = GetLambda(finalT, currentCouple);
684  if(preStepLambda<lx && 1 < verboseLevel) {
685  G4cout << "WARNING: for " << currentParticle->GetParticleName()
686  << " and " << GetProcessName()
687  << " E(MeV)= " << finalT/MeV
688  << " preLambda= " << preStepLambda << " < "
689  << lx << " (postLambda) "
690  << G4endl;
691  }
692 
693  if(preStepLambda*G4UniformRand() > lx) {
695  return &fParticleChange;
696  }
697  }
698 
700  if(!currentModel->IsActive(finalT)) { return &fParticleChange; }
701 
702  // define new weight for primary and secondaries
704  if(weightFlag) {
705  weight /= biasFactor;
707  }
708 
709  /*
710  if(0 < verboseLevel) {
711  G4cout << "G4VEmProcess::PostStepDoIt: Sample secondary; E= "
712  << finalT/MeV
713  << " MeV; model= (" << currentModel->LowEnergyLimit()
714  << ", " << currentModel->HighEnergyLimit() << ")"
715  << G4endl;
716  }
717  */
718 
719  // sample secondaries
720  secParticles.clear();
722  currentCouple,
723  track.GetDynamicParticle(),
724  (*theCuts)[currentCoupleIndex]);
725 
726  G4int num0 = secParticles.size();
727 
728  // splitting or Russian roulette
729  if(biasManager) {
731  G4double eloss = 0.0;
733  secParticles, track, currentModel, &fParticleChange, eloss,
735  step.GetPostStepPoint()->GetSafety());
736  if(eloss > 0.0) {
739  }
740  }
741  }
742 
743  // save secondaries
744  G4int num = secParticles.size();
745  if(num > 0) {
746 
749  G4double time = track.GetGlobalTime();
750 
751  for (G4int i=0; i<num; ++i) {
752  if (secParticles[i]) {
755  G4double e = dp->GetKineticEnergy();
756  G4bool good = true;
757  if(applyCuts) {
758  if (p == theGamma) {
759  if (e < (*theCutsGamma)[currentCoupleIndex]) { good = false; }
760 
761  } else if (p == theElectron) {
762  if (e < (*theCutsElectron)[currentCoupleIndex]) { good = false; }
763 
764  } else if (p == thePositron) {
765  if (electron_mass_c2 < (*theCutsGamma)[currentCoupleIndex] &&
767  good = false;
768  e += 2.0*electron_mass_c2;
769  }
770  }
771  // added secondary if it is good
772  }
773  if (good) {
774  G4Track* t = new G4Track(dp, time, track.GetPosition());
776  t->SetWeight(weight);
778 
779  // define type of secondary
781  else if(i < num0) {
782  if(p == theGamma) {
784  } else {
786  }
787  } else {
789  }
790 
791  //G4cout << "Secondary(post step) has weight " << t->GetWeight()
792  // << ", Ekin= " << t->GetKineticEnergy()/MeV << " MeV" <<G4endl;
793  } else {
794  delete dp;
795  edep += e;
796  }
797  }
798  }
800  }
801 
807  }
808 
809  // ClearNumberOfInteractionLengthLeft();
810  return &fParticleChange;
811 }
812 
813 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
814 
816  const G4String& directory,
817  G4bool ascii)
818 {
819  G4bool yes = true;
820 
821  if ( theLambdaTable && part == particle) {
822  const G4String name =
823  GetPhysicsTableFileName(part,directory,"Lambda",ascii);
824  yes = theLambdaTable->StorePhysicsTable(name,ascii);
825 
826  if ( yes ) {
827  G4cout << "Physics table is stored for " << particle->GetParticleName()
828  << " and process " << GetProcessName()
829  << " in the directory <" << directory
830  << "> " << G4endl;
831  } else {
832  G4cout << "Fail to store Physics Table for "
834  << " and process " << GetProcessName()
835  << " in the directory <" << directory
836  << "> " << G4endl;
837  }
838  }
839  if ( theLambdaTablePrim && part == particle) {
840  const G4String name =
841  GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
842  yes = theLambdaTablePrim->StorePhysicsTable(name,ascii);
843 
844  if ( yes ) {
845  G4cout << "Physics table prim is stored for "
847  << " and process " << GetProcessName()
848  << " in the directory <" << directory
849  << "> " << G4endl;
850  } else {
851  G4cout << "Fail to store Physics Table Prim for "
853  << " and process " << GetProcessName()
854  << " in the directory <" << directory
855  << "> " << G4endl;
856  }
857  }
858  return yes;
859 }
860 
861 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
862 
864  const G4String& directory,
865  G4bool ascii)
866 {
867  if(1 < verboseLevel) {
868  G4cout << "G4VEmProcess::RetrievePhysicsTable() for "
869  << part->GetParticleName() << " and process "
870  << GetProcessName() << G4endl;
871  }
872  G4bool yes = true;
873 
875  || particle != part) { return yes; }
876 
877  const G4String particleName = part->GetParticleName();
878  G4String filename;
879 
880  if(buildLambdaTable) {
881  filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
883  filename,ascii);
884  if ( yes ) {
885  if (0 < verboseLevel) {
886  G4cout << "Lambda table for " << particleName
887  << " is Retrieved from <"
888  << filename << ">"
889  << G4endl;
890  }
891  if(lManager->SplineFlag()) {
892  size_t n = theLambdaTable->length();
893  for(size_t i=0; i<n; ++i) {
894  if((* theLambdaTable)[i]) {
895  (* theLambdaTable)[i]->SetSpline(true);
896  }
897  }
898  }
899  } else {
900  if (1 < verboseLevel) {
901  G4cout << "Lambda table for " << particleName << " in file <"
902  << filename << "> is not exist"
903  << G4endl;
904  }
905  }
906  }
908  filename = GetPhysicsTableFileName(part,directory,"LambdaPrim",ascii);
910  filename,ascii);
911  if ( yes ) {
912  if (0 < verboseLevel) {
913  G4cout << "Lambda table prim for " << particleName
914  << " is Retrieved from <"
915  << filename << ">"
916  << G4endl;
917  }
918  if(lManager->SplineFlag()) {
919  size_t n = theLambdaTablePrim->length();
920  for(size_t i=0; i<n; ++i) {
921  if((* theLambdaTablePrim)[i]) {
922  (* theLambdaTablePrim)[i]->SetSpline(true);
923  }
924  }
925  }
926  } else {
927  if (1 < verboseLevel) {
928  G4cout << "Lambda table prim for " << particleName << " in file <"
929  << filename << "> is not exist"
930  << G4endl;
931  }
932  }
933  }
934 
935  return yes;
936 }
937 
938 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
939 
940 G4double
942  const G4MaterialCutsCouple* couple)
943 {
944  // Cross section per atom is calculated
945  DefineMaterial(couple);
946  G4double cross = 0.0;
948  cross = GetCurrentLambda(kineticEnergy);
949  } else {
950  SelectModel(kineticEnergy, currentCoupleIndex);
953  kineticEnergy);
954  }
955 
956  if(cross < 0.0) { cross = 0.0; }
957  return cross;
958 }
959 
960 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
961 
963  G4double,
965 {
966  *condition = NotForced;
967  return G4VEmProcess::MeanFreePath(track);
968 }
969 
970 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
971 
973 {
976  G4double x = DBL_MAX;
977  if(0.0 < preStepLambda) { x = 1.0/preStepLambda; }
978  return x;
979 }
980 
981 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
982 
983 G4double
985  G4double Z, G4double A, G4double cut)
986 {
987  SelectModel(kineticEnergy, currentCoupleIndex);
988  G4double x = 0.0;
989  if(currentModel) {
991  Z,A,cut);
992  }
993  return x;
994 }
995 
996 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
997 
999 {
1000  if(1 < verboseLevel) {
1001  G4cout << "### G4VEmProcess::FindLambdaMax: "
1002  << particle->GetParticleName()
1003  << " and process " << GetProcessName() << " " << G4endl;
1004  }
1005  size_t n = theLambdaTable->length();
1006  G4PhysicsVector* pv;
1007  G4double e, ss, emax, smax;
1008 
1009  size_t i;
1010 
1011  // first loop on existing vectors
1012  for (i=0; i<n; ++i) {
1013  pv = (*theLambdaTable)[i];
1014  if(pv) {
1015  size_t nb = pv->GetVectorLength();
1016  emax = DBL_MAX;
1017  smax = 0.0;
1018  if(nb > 0) {
1019  for (size_t j=0; j<nb; ++j) {
1020  e = pv->Energy(j);
1021  ss = (*pv)(j);
1022  if(ss > smax) {
1023  smax = ss;
1024  emax = e;
1025  }
1026  }
1027  }
1028  theEnergyOfCrossSectionMax[i] = emax;
1029  theCrossSectionMax[i] = smax;
1030  if(1 < verboseLevel) {
1031  G4cout << "For " << particle->GetParticleName()
1032  << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV
1033  << " lambda= " << smax << G4endl;
1034  }
1035  }
1036  }
1037  // second loop using base materials
1038  for (i=0; i<n; ++i) {
1039  pv = (*theLambdaTable)[i];
1040  if(!pv){
1041  G4int j = (*theDensityIdx)[i];
1043  theCrossSectionMax[i] = (*theDensityFactor)[i]*theCrossSectionMax[j];
1044  }
1045  }
1046 }
1047 
1048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1049 
1052 {
1053  G4PhysicsVector* v =
1055  v->SetSpline(lManager->SplineFlag());
1056  return v;
1057 }
1058 
1059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1060 
1062 {
1063  const G4Element* elm = 0;
1065  return elm;
1066 }
1067 
1068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1069 
1071 {
1072  if(f > 0.0) {
1073  biasFactor = f;
1074  weightFlag = flag;
1075  if(1 < verboseLevel) {
1076  G4cout << "### SetCrossSectionBiasingFactor: for "
1077  << particle->GetParticleName()
1078  << " and process " << GetProcessName()
1079  << " biasFactor= " << f << " weightFlag= " << flag
1080  << G4endl;
1081  }
1082  }
1083 }
1084 
1085 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1086 
1087 void
1089  G4bool flag)
1090 {
1091  if(!biasManager) { biasManager = new G4EmBiasingManager(); }
1092  if(1 < verboseLevel) {
1093  G4cout << "### ActivateForcedInteraction: for "
1094  << particle->GetParticleName()
1095  << " and process " << GetProcessName()
1096  << " length(mm)= " << length/mm
1097  << " in G4Region <" << r
1098  << "> weightFlag= " << flag
1099  << G4endl;
1100  }
1101  weightFlag = flag;
1103 }
1104 
1105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1106 
1107 void
1109  G4double factor,
1110  G4double energyLimit)
1111 {
1112  if (0.0 <= factor) {
1113 
1114  // Range cut can be applied only for e-
1115  if(0.0 == factor && secondaryParticle != G4Electron::Electron())
1116  { return; }
1117 
1118  if(!biasManager) { biasManager = new G4EmBiasingManager(); }
1119  biasManager->ActivateSecondaryBiasing(region, factor, energyLimit);
1120  if(1 < verboseLevel) {
1121  G4cout << "### ActivateSecondaryBiasing: for "
1122  << " process " << GetProcessName()
1123  << " factor= " << factor
1124  << " in G4Region <" << region
1125  << "> energyLimit(MeV)= " << energyLimit/MeV
1126  << G4endl;
1127  }
1128  }
1129 }
1130 
1131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1132 
1134 {
1137  minKinEnergy = e;
1138 }
1139 
1140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1141 
1143 {
1146  maxKinEnergy = e;
1147 }
1148 
1149 //....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
const G4Material * baseMaterial
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:230
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