G4EmUtility Class Reference

#include <G4EmUtility.hh>

Static Public Member Functions
static const G4Region *	FindRegion (const G4String &regionName, const G4int verbose=0)
static const G4Element *	SampleRandomElement (const G4Material *)
static const G4Isotope *	SampleRandomIsotope (const G4Element *)
static std::vector< G4double > *	FindCrossSectionMax (G4PhysicsTable *)
static std::vector< G4double > *	FindCrossSectionMax (G4VDiscreteProcess *, const G4ParticleDefinition *)
static std::vector< G4TwoPeaksXS * > *	FillPeaksStructure (G4PhysicsTable *, G4LossTableBuilder *)
static void	InitialiseElementSelectors (G4VEmModel *, const G4ParticleDefinition *, const G4DataVector &cuts, const G4double emin, const G4double emax)

Detailed Description

Definition at line 51 of file G4EmUtility.hh.

Member Function Documentation

FillPeaksStructure()

std::vector< G4TwoPeaksXS * > * G4EmUtility::FillPeaksStructure ( G4PhysicsTable *  p, G4LossTableBuilder *  bld  )

static

Definition at line 211 of file G4EmUtility.cc.

{
  std::vector<G4TwoPeaksXS*>* ptr = nullptr;
  if(nullptr == p) { return ptr; }
 
  const G4int n = (G4int)p->length();
  ptr = new std::vector<G4TwoPeaksXS*>;
  ptr->resize(n, nullptr);
 
  G4double e, ss, xs, ee;
  G4double e1peak, e1deep, e2peak, e2deep, e3peak;
  G4bool isDeep = false;
 
  // first loop on existing vectors
  for (G4int i=0; i<n; ++i) {
    const G4PhysicsVector* pv = (*p)[i];
    ee = xs = 0.0;
    e1peak = e1deep = e2peak = e2deep = e3peak = DBL_MAX;
    if(nullptr != pv) {
      G4int nb = (G4int)pv->GetVectorLength();
      for (G4int j=0; j<nb; ++j) {
    e = pv->Energy(j);
    ss = (*pv)(j);
    // find out 1st peak
    if(e1peak == DBL_MAX) {
      if(ss >= xs) {
        xs = ss;
        ee = e;
        continue;
      } else {
        e1peak = ee;
      }
    }
    // find out the deep
    if(e1deep == DBL_MAX) {
      if(ss <= xs) {
        xs = ss;
        ee = e;
        continue;
      } else {
        e1deep = ee;
        isDeep = true;
      }
    }
    // find out 2nd peak
    if(e2peak == DBL_MAX) {
      if(ss >= xs) {
        xs = ss;
        ee = e;
        continue;
      } else {
        e2peak = ee;
      }
    }
    if(e2deep == DBL_MAX) {
      if(ss <= xs) {
        xs = ss;
        ee = e;
        continue;
      } else {
        e2deep = ee;
        break;
      }
    }
    // find out 3d peak
    if(e3peak == DBL_MAX) {
      if(ss >= xs) {
        xs = ss;
        ee = e;
        continue;
      } else {
        e3peak = ee;
      }
    }
      }
    }
    G4TwoPeaksXS* x = (*ptr)[i];
    if(nullptr == x) { 
      x = new G4TwoPeaksXS(); 
      (*ptr)[i] = x;
    }
    x->e1peak = e1peak;
    x->e1deep = e1deep;
    x->e2peak = e2peak;
    x->e2deep = e2deep;
    x->e3peak = e3peak;
  }
  // case of no 1st peak in all vectors
  if(!isDeep) {
    delete ptr;
    ptr = nullptr;
    return ptr;
  }
  // check base particles
  if(!bld->GetBaseMaterialFlag()) { return ptr; }
 
  auto theDensityIdx = bld->GetCoupleIndexes();
  // second loop using base materials
  for (G4int i=0; i<n; ++i) {
    const G4PhysicsVector* pv = (*p)[i];
    if (nullptr == pv) {
      G4int j = (*theDensityIdx)[i];
      if(j == i) { continue; }
      G4TwoPeaksXS* x = (*ptr)[i];
      G4TwoPeaksXS* y = (*ptr)[j];
      if(nullptr == x) { 
    x = new G4TwoPeaksXS(); 
    (*ptr)[i] = x;
      }
      x->e1peak = y->e1peak;
      x->e1deep = y->e1deep;
      x->e2peak = y->e2peak;
      x->e2deep = y->e2deep;
      x->e3peak = y->e3peak;
    }
  }
  return ptr;
}

Referenced by G4VEnergyLossProcess::SetLambdaTable().

FindCrossSectionMax() [1/2]

std::vector< G4double > * G4EmUtility::FindCrossSectionMax ( G4PhysicsTable *  p )

static

Definition at line 104 of file G4EmUtility.cc.

{
  std::vector<G4double>* ptr = nullptr;
  if(nullptr == p) { return ptr; }
 
  const std::size_t n = p->length();
  ptr = new std::vector<G4double>;
  ptr->resize(n, DBL_MAX);
 
  G4bool isPeak = false;
  G4double e, ss, ee, xs;
 
  // first loop on existing vectors
  for (std::size_t i=0; i<n; ++i) {
    const G4PhysicsVector* pv = (*p)[i];
    xs = ee = 0.0;
    if(nullptr != pv) {
      G4int nb = (G4int)pv->GetVectorLength();
      for (G4int j=0; j<nb; ++j) {
    e = pv->Energy(j);
    ss = (*pv)(j);
    if(ss >= xs) {
      xs = ss;
      ee = e;
      continue;
    } else {
      isPeak = true;
      (*ptr)[i] = ee;
      break;
    }
      }
    }
  }
 
  // there is no peak for any material
  if(!isPeak) {
    delete ptr;
    ptr = nullptr;
  }
  return ptr;
}

Referenced by G4EmTableUtil::BuildEmProcess(), and G4VEnergyLossProcess::SetLambdaTable().

FindCrossSectionMax() [2/2]

std::vector< G4double > * G4EmUtility::FindCrossSectionMax ( G4VDiscreteProcess *  p, const G4ParticleDefinition *  part  )

static

Definition at line 149 of file G4EmUtility.cc.

{
  std::vector<G4double>* ptr = nullptr;
  if(nullptr == p || nullptr == part) { return ptr; }
  /*
  G4cout << "G4EmUtility::FindCrossSectionMax for "
     << p->GetProcessName() << " and " << part->GetParticleName() << G4endl;
  */
  G4EmParameters* theParameters = G4EmParameters::Instance();
  G4double tmin = theParameters->MinKinEnergy();
  G4double tmax = theParameters->MaxKinEnergy();
 
  const G4ProductionCutsTable* theCoupleTable=
        G4ProductionCutsTable::GetProductionCutsTable();
  std::size_t n = theCoupleTable->GetTableSize();
  ptr = new std::vector<G4double>;
  ptr->resize(n, DBL_MAX);
 
  G4bool isPeak = false;
  G4double scale = theParameters->NumberOfBinsPerDecade()/g4log10;
 
  G4double e, sig, ee, x, sm, em, emin, emax;
 
  // first loop on existing vectors
  for (std::size_t i=0; i<n; ++i) {
    auto couple = theCoupleTable->GetMaterialCutsCouple((G4int)i);
    emin = std::max(p->MinPrimaryEnergy(part, couple->GetMaterial()), tmin);
    emax = std::max(tmax, 2*emin);
    ee = G4Log(emax/emin);
 
    G4int nbin = G4lrint(ee*scale);
    if(nbin < 4) { nbin = 4; }
    x = G4Exp(ee/nbin);
    sm = 0.0;
    em = 0.0;
    e = emin;
    for(G4int j=0; j<=nbin; ++j) {
      sig = p->GetCrossSection(e, couple);
      if(sig >= sm) {
    em = e;
    sm = sig;
    e = (j+1 < nbin) ? e*x : emax;
      } else {
    isPeak = true;
    (*ptr)[i] = em;
    break;
      }
    }
    //G4cout << i << ".  em=" << em << " sm=" << sm << G4endl;
  }
  // there is no peak for any couple
  if(!isPeak) {
    delete ptr;
    ptr = nullptr;
  }
  return ptr;
}

FindRegion()

const G4Region * G4EmUtility::FindRegion ( const G4String &  regionName, const G4int  verbose = 0  )

static

Definition at line 47 of file G4EmUtility.cc.

{
  const G4Region* reg = nullptr;
  G4RegionStore* regStore = G4RegionStore::GetInstance();
  G4String r = regionName;
  if(r == "") { r = "DefaultRegionForTheWorld"; }
  reg = regStore->GetRegion(r, true); 
  if(nullptr == reg && verbose > 0) {
    G4cout << "### G4EmUtility WARNING: fails to find a region <"
           << r << G4endl;
  } else if(verbose > 1) {
    G4cout << "### G4EmUtility finds out G4Region <" << r << ">" 
           << G4endl;
  } 
  return reg;  
}

Referenced by G4EmConfigurator::AddModels(), G4EmDNAPhysicsActivator::ConstructProcess(), G4EmCalculator::FindRegion(), and G4EmConfigurator::PrepareModels().

InitialiseElementSelectors()

void G4EmUtility::InitialiseElementSelectors ( G4VEmModel *  mod, const G4ParticleDefinition *  part, const G4DataVector &  cuts, const G4double  emin, const G4double  emax  )

static

Definition at line 332 of file G4EmUtility.cc.

{
  // using spline for element selectors should be investigated in details
  // because small number of points may provide biased results
  // large number of points requires significant increase of memory
  G4bool spline = false;
 
  G4int nbinsPerDec = G4EmParameters::Instance()->NumberOfBinsPerDecade();
 
  G4ProductionCutsTable* theCoupleTable=
    G4ProductionCutsTable::GetProductionCutsTable();
  std::size_t numOfCouples = theCoupleTable->GetTableSize();
 
  // prepare vector
  auto elmSelectors = mod->GetElementSelectors();
  if(nullptr == elmSelectors) {
    elmSelectors = new std::vector<G4EmElementSelector*>;
  }
  std::size_t nSelectors = elmSelectors->size();
  if(numOfCouples > nSelectors) { 
    for(std::size_t i=nSelectors; i<numOfCouples; ++i) { 
      elmSelectors->push_back(nullptr); 
    }
    nSelectors = numOfCouples;
  }
 
  // initialise vector
  for(std::size_t i=0; i<numOfCouples; ++i) {
 
    // no need in element selectors for infinite cuts
    if(cuts[i] == DBL_MAX) { continue; }
   
    auto couple = theCoupleTable->GetMaterialCutsCouple((G4int)i); 
    auto mat = couple->GetMaterial();
    mod->SetCurrentCouple(couple);
 
    // selector already exist then delete
    delete (*elmSelectors)[i];
 
    G4double emin = std::max(elow, mod->MinPrimaryEnergy(mat, part, cuts[i]));
    G4double emax = std::max(ehigh, 10*emin);
    static const G4double invlog106 = 1.0/(6*G4Log(10.));
    G4int nbins = G4lrint(nbinsPerDec*G4Log(emax/emin)*invlog106);
    nbins = std::max(nbins, 3);
 
    (*elmSelectors)[i] = new G4EmElementSelector(mod,mat,nbins,
                         emin,emax,spline);
    ((*elmSelectors)[i])->Initialise(part, cuts[i]);
    /*      
      G4cout << "G4VEmModel::InitialiseElmSelectors i= " << i 
             << "  "  << part->GetParticleName() 
             << " for " << mod->GetName() << "  cut= " << cuts[i] 
             << "  " << (*elmSelectors)[i] << G4endl;      
      ((*elmSelectors)[i])->Dump(part);
    */
  }
  mod->SetElementSelectors(elmSelectors);
}

Referenced by G4VEmModel::InitialiseElementSelectors().

SampleRandomElement()

const G4Element * G4EmUtility::SampleRandomElement ( const G4Material *  mat )

static

Definition at line 66 of file G4EmUtility.cc.

{
  const G4Element* elm = mat->GetElement(0);
  std::size_t nElements = mat->GetNumberOfElements();
  if(1 < nElements) {
    G4double x = mat->GetTotNbOfElectPerVolume()*G4UniformRand();
    const G4double* y = mat->GetVecNbOfAtomsPerVolume();
    for(std::size_t i=0; i<nElements; ++i) {
      elm = mat->GetElement((G4int)i);
      x -= y[i]*elm->GetZ();
      if(x <= 0.0) { break; }
    }
  }
  return elm;
}

Referenced by G4VEmModel::GetCurrentElement().

SampleRandomIsotope()

const G4Isotope * G4EmUtility::SampleRandomIsotope ( const G4Element *  elm )

static

Definition at line 84 of file G4EmUtility.cc.

{
  const std::size_t ni = elm->GetNumberOfIsotopes();
  const G4Isotope* iso = elm->GetIsotope(0);
  if(ni > 1) {
    const G4double* ab = elm->GetRelativeAbundanceVector();
    G4double x = G4UniformRand();
    for(std::size_t idx=0; idx<ni; ++idx) {
      x -= ab[idx];
      if (x <= 0.0) { 
    iso = elm->GetIsotope((G4int)idx);
    break; 
      }
    }
  }
  return iso;
}

Referenced by G4VEmModel::GetCurrentIsotope().

---

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

• G4EmUtility.hh
• G4EmUtility.cc