Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4EmModelManager Class Reference

#include <G4EmModelManager.hh>

Public Member Functions

 G4EmModelManager ()
 
 ~G4EmModelManager ()
 
void Clear ()
 
const G4DataVectorInitialise (const G4ParticleDefinition *part, const G4ParticleDefinition *secPart, G4double minSubRange, G4int verb)
 
void FillDEDXVector (G4PhysicsVector *, const G4MaterialCutsCouple *, G4EmTableType t=fRestricted)
 
void FillLambdaVector (G4PhysicsVector *, const G4MaterialCutsCouple *, G4bool startFromNull=true, G4EmTableType t=fRestricted)
 
G4VEmModelGetModel (G4int, G4bool ver=false)
 
void AddEmModel (G4int, G4VEmModel *, G4VEmFluctuationModel *, const G4Region *)
 
void UpdateEmModel (const G4String &, G4double, G4double)
 
void DumpModelList (G4int verb)
 
G4VEmModelSelectModel (G4double &energy, size_t &index)
 
const G4DataVectorCuts () const
 
const G4DataVectorSubCutoff () const
 
void SetFluoFlag (G4bool val)
 
G4int NumberOfModels () const
 

Detailed Description

Definition at line 142 of file G4EmModelManager.hh.

Constructor & Destructor Documentation

G4EmModelManager::G4EmModelManager ( )

Definition at line 124 of file G4EmModelManager.cc.

References python.hepunit::mm.

124  :
125  nEmModels(0),
126  nRegions(0),
127  particle(0),
128  verboseLevel(0)
129 {
130  maxSubCutInRange = 0.7*mm;
131  models.reserve(4);
132  flucModels.reserve(4);
133  regions.reserve(4);
134  orderOfModels.reserve(4);
135  isUsed.reserve(4);
136  severalModels = true;
137  fluoFlag = false;
138  currRegionModel = 0;
139  currModel = 0;
140  theCuts = 0;
141  theCutsNew = 0;
142  theSubCuts = 0;
143 }
G4EmModelManager::~G4EmModelManager ( )

Definition at line 147 of file G4EmModelManager.cc.

References Clear().

148 {
149  verboseLevel = 0; // no verbosity at destruction
150  Clear();
151  delete theCutsNew;
152  delete theSubCuts;
153 }

Member Function Documentation

void G4EmModelManager::AddEmModel ( G4int  num,
G4VEmModel p,
G4VEmFluctuationModel fm,
const G4Region r 
)

Definition at line 173 of file G4EmModelManager.cc.

References G4VEmModel::DefineForRegion(), G4cout, and G4endl.

Referenced by G4VMultipleScattering::AddEmModel(), G4VEmProcess::AddEmModel(), G4VEnergyLossProcess::AddEmModel(), and G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel().

175 {
176  if(!p) {
177  G4cout << "G4EmModelManager::AddEmModel WARNING: no model defined."
178  << G4endl;
179  return;
180  }
181  models.push_back(p);
182  flucModels.push_back(fm);
183  regions.push_back(r);
184  orderOfModels.push_back(num);
185  isUsed.push_back(0);
186  p->DefineForRegion(r);
187  ++nEmModels;
188 }
virtual void DefineForRegion(const G4Region *)
Definition: G4VEmModel.cc:309
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
void G4EmModelManager::Clear ( )

Definition at line 157 of file G4EmModelManager.cc.

References G4cout, G4endl, and n.

Referenced by Initialise(), and ~G4EmModelManager().

158 {
159  if(1 < verboseLevel) {
160  G4cout << "G4EmModelManager::Clear()" << G4endl;
161  }
162  size_t n = setOfRegionModels.size();
163  if(n > 0) {
164  for(size_t i=0; i<n; ++i) {
165  delete setOfRegionModels[i];
166  setOfRegionModels[i] = 0;
167  }
168  }
169 }
G4GLOB_DLL std::ostream G4cout
const G4int n
#define G4endl
Definition: G4ios.hh:61
const G4DataVector * G4EmModelManager::Cuts ( ) const
inline

Definition at line 245 of file G4EmModelManager.hh.

246 {
247  return theCuts;
248 }
void G4EmModelManager::DumpModelList ( G4int  verb)

Definition at line 714 of file G4EmModelManager.cc.

References G4InuclParticleNames::an, G4VEmModel::DeexcitationFlag(), G4PhysicsVector::Energy(), G4BestUnit, G4cout, G4endl, G4VEmModel::GetAngularDistribution(), G4VEmModel::GetCrossSectionTable(), G4VEmAngularDistribution::GetName(), G4Region::GetName(), G4VEmModel::GetName(), G4PhysicsVector::GetVectorLength(), G4VEmModel::HighEnergyActivationLimit(), G4VEmModel::LowEnergyActivationLimit(), G4INCL::Math::max(), G4INCL::Math::min(), n, G4InuclParticleNames::nn, and test::v.

Referenced by G4VMultipleScattering::BuildPhysicsTable(), G4VMultipleScattering::PrintInfoDefinition(), and G4VEnergyLossProcess::PrintInfoDefinition().

715 {
716  if(verb == 0) { return; }
717  for(G4int i=0; i<nRegions; ++i) {
718  G4RegionModels* r = setOfRegionModels[i];
719  const G4Region* reg = r->Region();
720  G4int n = r->NumberOfModels();
721  if(n > 0) {
722  G4cout << " ===== EM models for the G4Region " << reg->GetName()
723  << " ======" << G4endl;;
724  for(G4int j=0; j<n; ++j) {
725  G4VEmModel* model = models[r->ModelIndex(j)];
726  G4double emin =
727  std::max(r->LowEdgeEnergy(j),model->LowEnergyActivationLimit());
728  G4double emax =
729  std::min(r->LowEdgeEnergy(j+1),model->HighEnergyActivationLimit());
730  G4cout << std::setw(20);
731  G4cout << model->GetName() << " : Emin= "
732  << std::setw(8) << G4BestUnit(emin,"Energy")
733  << " Emax= "
734  << std::setw(8) << G4BestUnit(emax,"Energy");
735  G4PhysicsTable* table = model->GetCrossSectionTable();
736  if(table) {
737  size_t kk = table->size();
738  for(size_t k=0; k<kk; ++k) {
739  G4PhysicsVector* v = (*table)[k];
740  if(v) {
741  G4int nn = v->GetVectorLength() - 1;
742  G4cout << " Table with " << nn << " bins Emin= "
743  << std::setw(6) << G4BestUnit(v->Energy(0),"Energy")
744  << " Emax= "
745  << std::setw(6) << G4BestUnit(v->Energy(nn),"Energy");
746  break;
747  }
748  }
749  }
751  if(an) { G4cout << " " << an->GetName(); }
752  if(fluoFlag && model->DeexcitationFlag()) {
753  G4cout << " FluoActive";
754  }
755  G4cout << G4endl;
756  }
757  }
758  if(1 == nEmModels) { break; }
759  }
760  if(theCutsNew) {
761  G4cout << " ===== Limit on energy threshold has been applied "
762  << G4endl;
763  }
764 }
G4double LowEnergyActivationLimit() const
Definition: G4VEmModel.hh:613
G4double HighEnergyActivationLimit() const
Definition: G4VEmModel.hh:606
const G4String & GetName() const
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:578
G4PhysicsTable * GetCrossSectionTable()
Definition: G4VEmModel.hh:784
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
const XML_Char XML_Content * model
const G4int n
G4double Energy(size_t index) const
G4bool DeexcitationFlag() const
Definition: G4VEmModel.hh:641
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
const G4String & GetName() const
Definition: G4VEmModel.hh:753
double G4double
Definition: G4Types.hh:76
const G4String & GetName() const
void G4EmModelManager::FillDEDXVector ( G4PhysicsVector aVector,
const G4MaterialCutsCouple couple,
G4EmTableType  t = fRestricted 
)

Definition at line 561 of file G4EmModelManager.cc.

References DBL_MAX, G4PhysicsVector::Energy(), fSubRestricted, fTotal, G4cout, G4endl, G4MaterialCutsCouple::GetIndex(), G4MaterialCutsCouple::GetMaterial(), G4Material::GetName(), G4ParticleDefinition::GetParticleName(), G4PhysicsVector::GetVectorLength(), G4InuclParticleNames::k0, python.hepunit::MeV, python.hepunit::mm, and G4PhysicsVector::PutValue().

Referenced by G4VEnergyLossProcess::BuildDEDXTable().

564 {
565  size_t i = couple->GetIndex();
566  G4double cut = (*theCuts)[i];
567  G4double emin = 0.0;
568 
569  if(fTotal == tType) { cut = DBL_MAX; }
570  else if(fSubRestricted == tType) {
571  emin = cut;
572  if(theSubCuts) { emin = (*theSubCuts)[i]; }
573  }
574 
575  if(1 < verboseLevel) {
576  G4cout << "G4EmModelManager::FillDEDXVector() for "
577  << couple->GetMaterial()->GetName()
578  << " cut(MeV)= " << cut
579  << " emin(MeV)= " << emin
580  << " Type " << tType
581  << " for " << particle->GetParticleName()
582  << G4endl;
583  }
584 
585  G4int reg = 0;
586  if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
587  const G4RegionModels* regModels = setOfRegionModels[reg];
588  G4int nmod = regModels->NumberOfModels();
589 
590  // Calculate energy losses vector
591 
592  //G4cout << "nmod= " << nmod << G4endl;
593  size_t totBinsLoss = aVector->GetVectorLength();
594  G4double del = 0.0;
595  G4int k0 = 0;
596 
597  for(size_t j=0; j<totBinsLoss; ++j) {
598 
599  G4double e = aVector->Energy(j);
600 
601  // Choose a model of energy losses
602  G4int k = 0;
603  if (nmod > 1) {
604  k = nmod;
605  do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
606  //G4cout << "k= " << k << G4endl;
607  if(k > 0 && k != k0) {
608  k0 = k;
609  G4double elow = regModels->LowEdgeEnergy(k);
610  G4double dedx1 = ComputeDEDX(models[regModels->ModelIndex(k-1)],
611  couple,elow,cut,emin);
612  G4double dedx2 = ComputeDEDX(models[regModels->ModelIndex(k)],
613  couple,elow,cut,emin);
614  del = 0.0;
615  if(dedx2 > 0.0) { del = (dedx1/dedx2 - 1.0)*elow; }
616  //G4cout << "elow= " << elow
617  // << " dedx1= " << dedx1 << " dedx2= " << dedx2 << G4endl;
618  }
619  }
620  G4double dedx =
621  ComputeDEDX(models[regModels->ModelIndex(k)],couple,e,cut,emin);
622  dedx *= (1.0 + del/e);
623 
624  if(2 < verboseLevel) {
625  G4cout << "Material= " << couple->GetMaterial()->GetName()
626  << " E(MeV)= " << e/MeV
627  << " dEdx(MeV/mm)= " << dedx*mm/MeV
628  << " del= " << del*mm/MeV<< " k= " << k
629  << " modelIdx= " << regModels->ModelIndex(k)
630  << G4endl;
631  }
632  if(dedx < 0.0) { dedx = 0.0; }
633  aVector->PutValue(j, dedx);
634  }
635 }
const G4String & GetName() const
Definition: G4Material.hh:176
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
void PutValue(size_t index, G4double theValue)
G4double Energy(size_t index) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
const G4Material * GetMaterial() const
void G4EmModelManager::FillLambdaVector ( G4PhysicsVector aVector,
const G4MaterialCutsCouple couple,
G4bool  startFromNull = true,
G4EmTableType  t = fRestricted 
)

Definition at line 639 of file G4EmModelManager.cc.

References G4VEmModel::CrossSection(), DBL_MAX, G4PhysicsVector::Energy(), fIsCrossSectionPrim, fSubRestricted, G4cout, G4endl, G4MaterialCutsCouple::GetIndex(), G4MaterialCutsCouple::GetMaterial(), G4PhysicsVector::GetMaxEnergy(), G4Material::GetName(), G4ParticleDefinition::GetParticleName(), G4PhysicsVector::GetVectorLength(), G4InuclParticleNames::k0, python.hepunit::MeV, python.hepunit::mm, and G4PhysicsVector::PutValue().

Referenced by G4VEnergyLossProcess::BuildLambdaTable().

643 {
644  size_t i = couple->GetIndex();
645  G4double cut = (*theCuts)[i];
646  G4double tmax = DBL_MAX;
647  if (fSubRestricted == tType) {
648  tmax = cut;
649  if(theSubCuts) { cut = (*theSubCuts)[i]; }
650  }
651 
652  G4int reg = 0;
653  if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; }
654  const G4RegionModels* regModels = setOfRegionModels[reg];
655  G4int nmod = regModels->NumberOfModels();
656  if(1 < verboseLevel) {
657  G4cout << "G4EmModelManager::FillLambdaVector() for "
658  << particle->GetParticleName()
659  << " in " << couple->GetMaterial()->GetName()
660  << " Emin(MeV)= " << aVector->Energy(0)
661  << " Emax(MeV)= " << aVector->GetMaxEnergy()
662  << " Type " << tType
663  << " nmod= " << nmod
664  << G4endl;
665  }
666 
667  // Calculate lambda vector
668  size_t totBinsLambda = aVector->GetVectorLength();
669  G4double del = 0.0;
670  G4int k0 = 0;
671  G4int k = 0;
672  G4VEmModel* mod = models[regModels->ModelIndex(0)];
673  for(size_t j=0; j<totBinsLambda; ++j) {
674 
675  G4double e = aVector->Energy(j);
676 
677  // Choose a model
678  if (nmod > 1) {
679  k = nmod;
680  do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k));
681  if(k > 0 && k != k0) {
682  k0 = k;
683  G4double elow = regModels->LowEdgeEnergy(k);
684  G4VEmModel* mod1 = models[regModels->ModelIndex(k-1)];
685  G4double xs1 = mod1->CrossSection(couple,particle,elow,cut,tmax);
686  mod = models[regModels->ModelIndex(k)];
687  G4double xs2 = mod->CrossSection(couple,particle,elow,cut,tmax);
688  del = 0.0;
689  if(xs2 > 0.0) { del = (xs1/xs2 - 1.0)*elow; }
690  //G4cout << "New model k=" << k << " E(MeV)= " << e/MeV
691  // << " Elow(MeV)= " << elow/MeV << " del= " << del << G4endl;
692  }
693  }
694  G4double cross = mod->CrossSection(couple,particle,e,cut,tmax);
695  cross *= (1.0 + del/e);
696  if(fIsCrossSectionPrim == tType) { cross *= e; }
697 
698  if(j==0 && startFromNull) { cross = 0.0; }
699 
700  if(2 < verboseLevel) {
701  G4cout << "FillLambdaVector: " << j << ". e(MeV)= " << e/MeV
702  << " cross(1/mm)= " << cross*mm
703  << " del= " << del*mm << " k= " << k
704  << " modelIdx= " << regModels->ModelIndex(k)
705  << G4endl;
706  }
707  if(cross < 0.0) { cross = 0.0; }
708  aVector->PutValue(j, cross);
709  }
710 }
G4double GetMaxEnergy() const
const G4String & GetName() const
Definition: G4Material.hh:176
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
G4double CrossSection(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:467
void PutValue(size_t index, G4double theValue)
G4double Energy(size_t index) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
const G4Material * GetMaterial() const
G4VEmModel * G4EmModelManager::GetModel ( G4int  i,
G4bool  ver = false 
)

Definition at line 211 of file G4EmModelManager.cc.

References G4cout, G4endl, and G4ParticleDefinition::GetParticleName().

Referenced by G4VMultipleScattering::GetModelByIndex(), G4VEmProcess::GetModelByIndex(), G4VEnergyLossProcess::GetModelByIndex(), G4VMultipleScattering::PreparePhysicsTable(), G4VEmProcess::PreparePhysicsTable(), G4VEnergyLossProcess::PreparePhysicsTable(), and G4VMultipleScattering::StorePhysicsTable().

212 {
213  G4VEmModel* model = 0;
214  if(i >= 0 && i < nEmModels) { model = models[i]; }
215  else if(verboseLevel > 0 && ver) {
216  G4cout << "G4EmModelManager::GetModel WARNING: "
217  << "index " << i << " is wrong Nmodels= "
218  << nEmModels;
219  if(particle) G4cout << " for " << particle->GetParticleName();
220  G4cout<< G4endl;
221  }
222  return model;
223 }
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
const XML_Char XML_Content * model
#define G4endl
Definition: G4ios.hh:61
const G4DataVector * G4EmModelManager::Initialise ( const G4ParticleDefinition part,
const G4ParticleDefinition secPart,
G4double  minSubRange,
G4int  verb 
)

Definition at line 228 of file G4EmModelManager.cc.

References Clear(), G4ProductionCutsTable::ConvertRangeToEnergy(), DBL_MAX, python.hepunit::eV, FatalException, G4cout, G4endl, G4Exception(), G4Gamma::Gamma(), G4ProductionCutsTable::GetEnergyCutsVector(), G4RegionStore::GetInstance(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4Region::GetName(), G4Material::GetName(), G4VEmModel::GetName(), G4ParticleDefinition::GetParticleName(), G4ProductionCuts::GetProductionCut(), G4MaterialCutsCouple::GetProductionCuts(), G4ProductionCutsTable::GetProductionCutsTable(), G4RegionStore::GetRegion(), G4ProductionCutsTable::GetTableSize(), G4VEmModel::HighEnergyLimit(), G4VEmModel::Initialise(), G4VEmModel::LowEnergyLimit(), eplot::material, G4INCL::Math::max(), python.hepunit::MeV, G4INCL::Math::min(), G4VEmModel::MinEnergyCut(), n, G4InuclParticleNames::nn, and G4Positron::Positron().

Referenced by G4AdjointBremsstrahlungModel::AdjointCrossSection(), G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(), G4VMultipleScattering::PreparePhysicsTable(), G4VEmProcess::PreparePhysicsTable(), and G4VEnergyLossProcess::PreparePhysicsTable().

232 {
233  verboseLevel = val;
234  G4String partname = p->GetParticleName();
235  // if(partname == "proton") verboseLevel = 2;
236  //else verboseLevel = 0;
237  if(1 < verboseLevel) {
238  G4cout << "G4EmModelManager::Initialise() for "
239  << partname << G4endl;
240  }
241  // Are models defined?
242  if(nEmModels < 1) {
244  ed << "No models found out for " << p->GetParticleName()
245  << " !";
246  G4Exception("G4EmModelManager::Initialise","em0002",
247  FatalException, ed);
248  }
249 
250  particle = p;
251  Clear(); // needed if run is not first
252 
253  G4RegionStore* regionStore = G4RegionStore::GetInstance();
254  const G4Region* world =
255  regionStore->GetRegion("DefaultRegionForTheWorld", false);
256 
257  // Identify the list of regions with different set of models
258  nRegions = 1;
259  std::vector<const G4Region*> setr;
260  setr.push_back(world);
261  G4bool isWorld = false;
262 
263  for (G4int ii=0; ii<nEmModels; ++ii) {
264  const G4Region* r = regions[ii];
265  if ( r == 0 || r == world) {
266  isWorld = true;
267  regions[ii] = world;
268  } else {
269  G4bool newRegion = true;
270  if (nRegions>1) {
271  for (G4int j=1; j<nRegions; ++j) {
272  if ( r == setr[j] ) { newRegion = false; }
273  }
274  }
275  if (newRegion) {
276  setr.push_back(r);
277  nRegions++;
278  }
279  }
280  }
281  // Are models defined?
282  if(!isWorld) {
284  ed << "No models defined for the World volume for "
285  << p->GetParticleName() << " !";
286  G4Exception("G4EmModelManager::Initialise","em0002",
287  FatalException, ed);
288  }
289 
290  G4ProductionCutsTable* theCoupleTable=
292  size_t numOfCouples = theCoupleTable->GetTableSize();
293 
294  // prepare vectors, shortcut for the case of only 1 model
295  // or only one region
296  if(nRegions > 1 && nEmModels > 1) {
297  idxOfRegionModels.resize(numOfCouples,0);
298  setOfRegionModels.resize((size_t)nRegions,0);
299  } else {
300  idxOfRegionModels.resize(1,0);
301  setOfRegionModels.resize(1,0);
302  }
303 
304  std::vector<G4int> modelAtRegion(nEmModels);
305  std::vector<G4int> modelOrd(nEmModels);
306  G4DataVector eLow(nEmModels+1);
307  G4DataVector eHigh(nEmModels);
308 
309  if(1 < verboseLevel) {
310  G4cout << " Nregions= " << nRegions
311  << " Nmodels= " << nEmModels << G4endl;
312  }
313 
314  // Order models for regions
315  for (G4int reg=0; reg<nRegions; ++reg) {
316  const G4Region* region = setr[reg];
317  G4int n = 0;
318 
319  for (G4int ii=0; ii<nEmModels; ++ii) {
320 
321  G4VEmModel* model = models[ii];
322  if ( region == regions[ii] ) {
323 
324  G4double tmin = model->LowEnergyLimit();
325  G4double tmax = model->HighEnergyLimit();
326  G4int ord = orderOfModels[ii];
327  G4bool push = true;
328  G4bool insert = false;
329  G4int idx = n;
330 
331  if(1 < verboseLevel) {
332  G4cout << "Model #" << ii
333  << " <" << model->GetName() << "> for region <";
334  if (region) G4cout << region->GetName();
335  G4cout << "> "
336  << " tmin(MeV)= " << tmin/MeV
337  << "; tmax(MeV)= " << tmax/MeV
338  << "; order= " << ord
339  << G4endl;
340  }
341 
342  if(n > 0) {
343 
344  // extend energy range to previous models
345  tmin = std::min(tmin, eHigh[n-1]);
346  tmax = std::max(tmax, eLow[0]);
347  //G4cout << "tmin= " << tmin << " tmax= "
348  // << tmax << " ord= " << ord <<G4endl;
349  // empty energy range
350  if( tmax - tmin <= eV) push = false;
351  // low-energy model
352  else if (tmax == eLow[0]) {
353  push = false;
354  insert = true;
355  idx = 0;
356  // resolve intersections
357  } else if(tmin < eHigh[n-1]) {
358  // compare order
359  for(G4int k=0; k<n; ++k) {
360  // new model has lower application
361  if(ord >= modelOrd[k]) {
362  if(tmin < eHigh[k] && tmin >= eLow[k]) tmin = eHigh[k];
363  if(tmax <= eHigh[k] && tmax > eLow[k]) tmax = eLow[k];
364  if(tmax > eHigh[k] && tmin < eLow[k]) {
365  if(tmax - eHigh[k] > eLow[k] - tmin) tmin = eHigh[k];
366  else tmax = eLow[k];
367  }
368  if( tmax - tmin <= eV) {
369  push = false;
370  break;
371  }
372  }
373  }
374  //G4cout << "tmin= " << tmin << " tmax= "
375  // << tmax << " push= " << push << " idx= " << idx <<G4endl;
376  if(push) {
377  if (tmax == eLow[0]) {
378  push = false;
379  insert = true;
380  idx = 0;
381  // continue resolve intersections
382  } else if(tmin < eHigh[n-1]) {
383  // last energy interval
384  if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) {
385  eHigh[n-1] = tmin;
386  // first energy interval
387  } else if(tmin <= eLow[0] && tmax < eHigh[0]) {
388  eLow[0] = tmax;
389  push = false;
390  insert = true;
391  idx = 0;
392  } else {
393  // find energy interval to replace
394  for(G4int k=0; k<n; ++k) {
395  if(tmin <= eLow[k] && tmax >= eHigh[k]) {
396  push = false;
397  modelAtRegion[k] = ii;
398  modelOrd[k] = ord;
399  isUsed[ii] = 1;
400  }
401  }
402  }
403  }
404  }
405  }
406  }
407  if(insert) {
408  for(G4int k=n-1; k>=idx; --k) {
409  modelAtRegion[k+1] = modelAtRegion[k];
410  modelOrd[k+1] = modelOrd[k];
411  eLow[k+1] = eLow[k];
412  eHigh[k+1] = eHigh[k];
413  }
414  }
415  //G4cout << "push= " << push << " insert= " << insert
416  //<< " idx= " << idx <<G4endl;
417  if (push || insert) {
418  ++n;
419  modelAtRegion[idx] = ii;
420  modelOrd[idx] = ord;
421  eLow[idx] = tmin;
422  eHigh[idx] = tmax;
423  isUsed[ii] = 1;
424  }
425  }
426  }
427  eLow[0] = 0.0;
428  eLow[n] = eHigh[n-1];
429 
430  if(1 < verboseLevel) {
431  G4cout << "New G4RegionModels set with " << n << " models for region <";
432  if (region) { G4cout << region->GetName(); }
433  G4cout << "> Elow(MeV)= ";
434  for(G4int ii=0; ii<=n; ++ii) {G4cout << eLow[ii]/MeV << " ";}
435  G4cout << G4endl;
436  }
437  G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region);
438  setOfRegionModels[reg] = rm;
439  // shortcut
440  if(1 == nEmModels) { break; }
441  }
442 
443  currRegionModel = setOfRegionModels[0];
444  currModel = models[0];
445 
446  // Access to materials and build cuts
447  size_t idx = 1;
448  if(secondaryParticle) {
449  if( secondaryParticle == G4Gamma::Gamma() ) { idx = 0; }
450  else if( secondaryParticle == G4Positron::Positron()) { idx = 2; }
451  }
452 
453  theCuts =
454  static_cast<const G4DataVector*>(theCoupleTable->GetEnergyCutsVector(idx));
455 
456  // for the second run the check on cuts should be repeated
457  if(theCutsNew) { *theCutsNew = *theCuts; }
458 
459  if(minSubRange < 1.0) {
460  if( !theSubCuts ) { theSubCuts = new G4DataVector(); }
461  theSubCuts->resize(numOfCouples,DBL_MAX);
462  }
463 
464  // G4cout << "========Start define cuts" << G4endl;
465  // define cut values
466  for(size_t i=0; i<numOfCouples; ++i) {
467 
468  const G4MaterialCutsCouple* couple =
469  theCoupleTable->GetMaterialCutsCouple(i);
470  const G4Material* material = couple->GetMaterial();
471  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
472 
473  G4int reg = 0;
474  if(nRegions > 1 && nEmModels > 1) {
475  reg = nRegions;
476  do {--reg;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
477  idxOfRegionModels[i] = reg;
478  }
479  if(1 < verboseLevel) {
480  G4cout << "G4EmModelManager::Initialise() for "
481  << material->GetName()
482  << " indexOfCouple= " << i
483  << " indexOfRegion= " << reg
484  << G4endl;
485  }
486 
487  G4double cut = (*theCuts)[i];
488  if(secondaryParticle) {
489 
490  // compute subcut
491  if( cut < DBL_MAX && minSubRange < 1.0) {
492  G4double subcut = minSubRange*cut;
493  G4double rcut = std::min(minSubRange*pcuts->GetProductionCut(idx),
494  maxSubCutInRange);
495  G4double tcutmax =
496  theCoupleTable->ConvertRangeToEnergy(secondaryParticle,
497  material,rcut);
498  if(tcutmax < subcut) { subcut = tcutmax; }
499  (*theSubCuts)[i] = subcut;
500  }
501 
502  // note that idxOfRegionModels[] not always filled
503  G4int inn = 0;
504  G4int nnm = 1;
505  if(nRegions > 1 && nEmModels > 1) {
506  inn = idxOfRegionModels[i];
507  }
508  // check cuts and introduce upper limits
509  //G4cout << "idx= " << i << " cut(keV)= " << cut/keV << G4endl;
510  currRegionModel = setOfRegionModels[inn];
511  nnm = currRegionModel->NumberOfModels();
512 
513  //G4cout << "idx= " << i << " Nmod= " << nnm << G4endl;
514 
515  for(G4int jj=0; jj<nnm; ++jj) {
516  //G4cout << "jj= " << jj << " modidx= "
517  // << currRegionModel->ModelIndex(jj) << G4endl;
518  currModel = models[currRegionModel->ModelIndex(jj)];
519  G4double cutlim = currModel->MinEnergyCut(particle,couple);
520  if(cutlim > cut) {
521  if(!theCutsNew) { theCutsNew = new G4DataVector(*theCuts); }
522  (*theCutsNew)[i] = cutlim;
523  /*
524  G4cout << "### " << partname << " energy loss model in "
525  << material->GetName()
526  << " Cut was changed from " << cut/keV << " keV to "
527  << cutlim/keV << " keV " << " due to "
528  << currModel->GetName() << G4endl;
529  */
530  }
531  }
532  }
533  }
534  if(theCutsNew) { theCuts = theCutsNew; }
535 
536  // initialize models
537  G4int nn = 0;
538  severalModels = true;
539  for(G4int jj=0; jj<nEmModels; ++jj) {
540  if(1 == isUsed[jj]) {
541  ++nn;
542  currModel = models[jj];
543  currModel->Initialise(particle, *theCuts);
544  if(flucModels[jj]) { flucModels[jj]->InitialiseMe(particle); }
545  }
546  }
547  if(1 == nn) { severalModels = false; }
548 
549  if(1 < verboseLevel) {
550  G4cout << "G4EmModelManager for " << partname
551  << " is initialised; nRegions= " << nRegions
552  << " severalModels: " << severalModels
553  << G4endl;
554  }
555 
556  return theCuts;
557 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
const G4String & GetName() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
G4double GetProductionCut(G4int index) const
const char * p
Definition: xmltok.h:285
const G4String & GetName() const
Definition: G4Material.hh:176
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
Definition: G4VEmModel.cc:364
int G4int
Definition: G4Types.hh:78
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)=0
static G4RegionStore * GetInstance()
string material
Definition: eplot.py:19
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
const XML_Char XML_Content * model
const G4int n
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4ProductionCutsTable * GetProductionCutsTable()
static G4Positron * Positron()
Definition: G4Positron.cc:94
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
const G4String & GetName() const
Definition: G4VEmModel.hh:753
G4double ConvertRangeToEnergy(const G4ParticleDefinition *particle, const G4Material *material, G4double range)
double G4double
Definition: G4Types.hh:76
G4ProductionCuts * GetProductionCuts() const
#define DBL_MAX
Definition: templates.hh:83
const G4Material * GetMaterial() const
G4int G4EmModelManager::NumberOfModels ( ) const
inline
G4VEmModel * G4EmModelManager::SelectModel ( G4double energy,
size_t &  index 
)
inline

Definition at line 231 of file G4EmModelManager.hh.

Referenced by G4VMultipleScattering::SelectModel(), G4VEmProcess::SelectModel(), G4VEnergyLossProcess::SelectModel(), G4VEmProcess::SelectModelForMaterial(), and G4VEnergyLossProcess::SelectModelForMaterial().

233 {
234  if(severalModels) {
235  if(nRegions > 1) {
236  currRegionModel = setOfRegionModels[idxOfRegionModels[index]];
237  }
238  currModel = models[currRegionModel->SelectIndex(kinEnergy)];
239  }
240  return currModel;
241 }
void G4EmModelManager::SetFluoFlag ( G4bool  val)
inline

Definition at line 259 of file G4EmModelManager.hh.

Referenced by G4VEmProcess::PreparePhysicsTable().

260 {
261  fluoFlag = val;
262 }
const G4DataVector * G4EmModelManager::SubCutoff ( ) const
inline

Definition at line 252 of file G4EmModelManager.hh.

Referenced by G4VEnergyLossProcess::PreparePhysicsTable().

253 {
254  return theSubCuts;
255 }
void G4EmModelManager::UpdateEmModel ( const G4String nam,
G4double  emin,
G4double  emax 
)

Definition at line 192 of file G4EmModelManager.cc.

References G4cout, and G4endl.

Referenced by G4VEmProcess::UpdateEmModel(), and G4VEnergyLossProcess::UpdateEmModel().

194 {
195  if (nEmModels > 0) {
196  for(G4int i=0; i<nEmModels; ++i) {
197  if(nam == models[i]->GetName()) {
198  models[i]->SetLowEnergyLimit(emin);
199  models[i]->SetHighEnergyLimit(emax);
200  break;
201  }
202  }
203  }
204  G4cout << "G4EmModelManager::UpdateEmModel WARNING: no model <"
205  << nam << "> is found out"
206  << G4endl;
207 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

The documentation for this class was generated from the following files: