Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4LossTableBuilder Class Reference

#include <G4LossTableBuilder.hh>

Public Member Functions

 G4LossTableBuilder ()
 
virtual ~G4LossTableBuilder ()
 
void BuildDEDXTable (G4PhysicsTable *dedxTable, const std::vector< G4PhysicsTable * > &)
 
void BuildRangeTable (const G4PhysicsTable *dedxTable, G4PhysicsTable *rangeTable, G4bool isIonisation=false)
 
void BuildInverseRangeTable (const G4PhysicsTable *rangeTable, G4PhysicsTable *invRangeTable, G4bool isIonisation=false)
 
G4PhysicsTableBuildTableForModel (G4PhysicsTable *table, G4VEmModel *model, const G4ParticleDefinition *, G4double emin, G4double emax, G4bool spline)
 
void InitialiseBaseMaterials (G4PhysicsTable *table)
 
const std::vector< G4int > * GetCoupleIndexes ()
 
const std::vector< G4double > * GetDensityFactors ()
 
G4bool GetFlag (size_t idx) const
 
void SetSplineFlag (G4bool flag)
 
void SetInitialisationFlag (G4bool flag)
 

Detailed Description

Definition at line 62 of file G4LossTableBuilder.hh.

Constructor & Destructor Documentation

G4LossTableBuilder::G4LossTableBuilder ( )

Definition at line 73 of file G4LossTableBuilder.cc.

74 {
75  theParameters = G4EmParameters::Instance();
76  splineFlag = true;
77  isInitialized = false;
78 
79  theDensityFactor = new std::vector<G4double>;
80  theDensityIdx = new std::vector<G4int>;
81  theFlag = new std::vector<G4bool>;
82 }
static G4EmParameters * Instance()

Here is the call graph for this function:

G4LossTableBuilder::~G4LossTableBuilder ( )
virtual

Definition at line 86 of file G4LossTableBuilder.cc.

87 {
88  delete theDensityFactor;
89  delete theDensityIdx;
90  delete theFlag;
91 }

Member Function Documentation

void G4LossTableBuilder::BuildDEDXTable ( G4PhysicsTable dedxTable,
const std::vector< G4PhysicsTable * > &  list 
)

Definition at line 96 of file G4LossTableBuilder.cc.

98 {
99  size_t n_processes = list.size();
100  //G4cout << "Nproc= " << n_processes << " Ncoup= "
101  //<< dedxTable->size() << G4endl;
102  if(1 >= n_processes) { return; }
103 
104  size_t nCouples = dedxTable->size();
105  if(0 >= nCouples) { return; }
106 
107  for (size_t i=0; i<nCouples; ++i) {
108  // if ((*theFlag)[i]) {
109  G4PhysicsLogVector* pv0 =
110  static_cast<G4PhysicsLogVector*>((*(list[0]))[i]);
111  if(pv0) {
112  size_t npoints = pv0->GetVectorLength();
113  G4PhysicsLogVector* pv = new G4PhysicsLogVector(*pv0);
114  pv->SetSpline(splineFlag);
115  for (size_t j=0; j<npoints; ++j) {
116  G4double dedx = 0.0;
117  for (size_t k=0; k<n_processes; ++k) {
118  G4PhysicsVector* pv1 = (*(list[k]))[i];
119  dedx += (*pv1)[j];
120  }
121  pv->PutValue(j, dedx);
122  }
123  if(splineFlag) { pv->FillSecondDerivatives(); }
124  G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
125  }
126  }
127 }
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
size_t GetVectorLength() const
void FillSecondDerivatives()
void SetSpline(G4bool)
void PutValue(size_t index, G4double theValue)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4LossTableBuilder::BuildInverseRangeTable ( const G4PhysicsTable rangeTable,
G4PhysicsTable invRangeTable,
G4bool  isIonisation = false 
)

Definition at line 218 of file G4LossTableBuilder.cc.

222 {
223  size_t nCouples = rangeTable->size();
224  if(0 >= nCouples) { return; }
225 
226  for (size_t i=0; i<nCouples; ++i) {
227 
228  if(isIonisation) {
229  if( !(*theFlag)[i] ) { continue; }
230  }
231  G4PhysicsVector* pv = (*rangeTable)[i];
232  size_t npoints = pv->GetVectorLength();
233  G4double rlow = (*pv)[0];
234  G4double rhigh = (*pv)[npoints-1];
235 
236  delete (*invRangeTable)[i];
237  G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh);
238  v->SetSpline(splineFlag);
239 
240  for (size_t j=0; j<npoints; ++j) {
241  G4double e = pv->Energy(j);
242  G4double r = (*pv)[j];
243  v->PutValues(j,r,e);
244  }
245  if(splineFlag) { v->FillSecondDerivatives(); }
246 
247  G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);
248  }
249 }
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
size_t GetVectorLength() const
void FillSecondDerivatives()
void SetSpline(G4bool)
G4double Energy(size_t index) const
void PutValues(size_t index, G4double e, G4double dataValue)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4LossTableBuilder::BuildRangeTable ( const G4PhysicsTable dedxTable,
G4PhysicsTable rangeTable,
G4bool  isIonisation = false 
)

Definition at line 131 of file G4LossTableBuilder.cc.

135 {
136  size_t nCouples = dedxTable->size();
137  if(0 >= nCouples) { return; }
138 
139  size_t n = 100;
140  G4double del = 1.0/(G4double)n;
141 
142  for (size_t i=0; i<nCouples; ++i) {
143  if(isIonisation) {
144  if( !(*theFlag)[i] ) { continue; }
145  }
146  G4PhysicsLogVector* pv = static_cast<G4PhysicsLogVector*>((*dedxTable)[i]);
147  size_t npoints = pv->GetVectorLength();
148  size_t bin0 = 0;
149  G4double elow = pv->Energy(0);
150  G4double ehigh = pv->Energy(npoints-1);
151  G4double dedx1 = (*pv)[0];
152 
153  //G4cout << "i= " << i << "npoints= " << npoints << " dedx1= "
154  //<< dedx1 << G4endl;
155 
156  // protection for specific cases dedx=0
157  if(dedx1 == 0.0) {
158  for (size_t k=1; k<npoints; ++k) {
159  bin0++;
160  elow = pv->Energy(k);
161  dedx1 = (*pv)[k];
162  if(dedx1 > 0.0) { break; }
163  }
164  npoints -= bin0;
165  }
166  //G4cout<<"New Range vector" << G4endl;
167  //G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh
168  // <<" bin0= " << bin0 <<G4endl;
169 
170  // initialisation of a new vector
171  if(npoints < 2) { npoints = 2; }
172 
173  delete (*rangeTable)[i];
175  if(0 == bin0) { v = new G4PhysicsLogVector(*pv); }
176  else { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); }
177 
178  // dedx is exact zero cannot build range table
179  if(2 == npoints) {
180  v->PutValue(0,1000.);
181  v->PutValue(1,2000.);
182  G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
183  return;
184  }
185  v->SetSpline(splineFlag);
186 
187  // assumed dedx proportional to beta
188  G4double energy1 = v->Energy(0);
189  G4double range = 2.*energy1/dedx1;
190  //G4cout << "range0= " << range << G4endl;
191  v->PutValue(0,range);
192 
193  for (size_t j=1; j<npoints; ++j) {
194 
195  G4double energy2 = v->Energy(j);
196  G4double de = (energy2 - energy1) * del;
197  G4double energy = energy2 + de*0.5;
198  G4double sum = 0.0;
199  //G4cout << "j= " << j << " e1= " << energy1 << " e2= " << energy2
200  // << " n= " << n << G4endl;
201  for (size_t k=0; k<n; ++k) {
202  energy -= de;
203  dedx1 = pv->Value(energy);
204  if(dedx1 > 0.0) { sum += de/dedx1; }
205  }
206  range += sum;
207  v->PutValue(j,range);
208  energy1 = energy2;
209  }
210  if(splineFlag) { v->FillSecondDerivatives(); }
211  G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
212  }
213 }
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
size_t GetVectorLength() const
void FillSecondDerivatives()
void SetSpline(G4bool)
void PutValue(size_t index, G4double theValue)
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double energy(const ThreeVector &p, const G4double m)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4PhysicsTable * G4LossTableBuilder::BuildTableForModel ( G4PhysicsTable table,
G4VEmModel model,
const G4ParticleDefinition part,
G4double  emin,
G4double  emax,
G4bool  spline 
)

Definition at line 407 of file G4LossTableBuilder.cc.

412 {
413  // check input
415  if(!table) { return table; }
416  if(emin >= emax) {
417  table->clearAndDestroy();
418  delete table;
419  table = 0;
420  return table;
421  }
423 
424  G4int nbins = theParameters->NumberOfBinsPerDecade();
425 
426  // Access to materials
427  const G4ProductionCutsTable* theCoupleTable=
429  size_t numOfCouples = theCoupleTable->GetTableSize();
430 
431  G4PhysicsLogVector* aVector = 0;
432 
433  for(size_t i=0; i<numOfCouples; ++i) {
434 
435  //G4cout<< "i= " << i << " Flag= " << GetFlag(i) << G4endl;
436 
437  if (GetFlag(i)) {
438 
439  // create physics vector and fill it
440  const G4MaterialCutsCouple* couple =
441  theCoupleTable->GetMaterialCutsCouple(i);
442  delete (*table)[i];
443 
444  // if start from zero then change the scale
445 
446  const G4Material* mat = couple->GetMaterial();
447 
448  G4double tmin = std::max(emin,model->MinPrimaryEnergy(mat,part));
449  if(0.0 >= tmin) { tmin = eV; }
450  G4int n = nbins;
451 
452  if(tmin >= emax) {
453  aVector = 0;
454  } else {
455  n *= G4int(std::log10(emax/tmin) + 0.5);
456  if(n < 3) { n = 3; }
457  aVector = new G4PhysicsLogVector(tmin, emax, n);
458  }
459 
460  if(aVector) {
461  aVector->SetSpline(spline);
462  for(G4int j=0; j<=n; ++j) {
463  aVector->PutValue(j, model->Value(couple, part,
464  aVector->Energy(j)));
465  }
466  if(spline) { aVector->FillSecondDerivatives(); }
467  }
468  G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
469  }
470  }
471  /*
472  G4cout << "G4LossTableBuilder::BuildTableForModel done for "
473  << part->GetParticleName() << " and "<< model->GetName()
474  << " " << table << G4endl;
475  */
476  return table;
477 }
G4int NumberOfBinsPerDecade() const
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
Definition: G4VEmModel.cc:386
G4bool GetFlag(size_t idx) const
int G4int
Definition: G4Types.hh:78
void FillSecondDerivatives()
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
void SetSpline(G4bool)
void PutValue(size_t index, G4double theValue)
static constexpr double eV
Definition: G4SIunits.hh:215
G4double Energy(size_t index) const
static const G4double emax
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void InitialiseBaseMaterials(G4PhysicsTable *table)
double G4double
Definition: G4Types.hh:76
void clearAndDestroy()
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:377
const G4Material * GetMaterial() const

Here is the call graph for this function:

Here is the caller graph for this function:

const std::vector< G4int > * G4LossTableBuilder::GetCoupleIndexes ( )
inline

Definition at line 126 of file G4LossTableBuilder.hh.

127 {
128  if(theDensityIdx->size() == 0) { InitialiseCouples(); }
129  return theDensityIdx;
130 }

Here is the caller graph for this function:

const std::vector< G4double > * G4LossTableBuilder::GetDensityFactors ( )
inline

Definition at line 133 of file G4LossTableBuilder.hh.

134 {
135  if(theDensityIdx->size() == 0) { InitialiseCouples(); }
136  return theDensityFactor;
137 }

Here is the caller graph for this function:

G4bool G4LossTableBuilder::GetFlag ( size_t  idx) const
inline

Definition at line 139 of file G4LossTableBuilder.hh.

140 {
141  return (*theFlag)[idx];
142 }

Here is the caller graph for this function:

void G4LossTableBuilder::InitialiseBaseMaterials ( G4PhysicsTable table)

Definition at line 254 of file G4LossTableBuilder.cc.

255 {
256  size_t nCouples = table->size();
257  size_t nFlags = theFlag->size();
258 
259  if(nCouples == nFlags && isInitialized) { return; }
260 
261  isInitialized = true;
262 
263  //G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials Ncouples= "
264  // << nCouples << " FlagSize= " << nFlags << G4endl;
265 
266  // variable density check
267  const G4ProductionCutsTable* theCoupleTable=
269 
270  /*
271  for(size_t i=0; i<nFlags; ++i) {
272  G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i]
273  << " tableFlag= " << table->GetFlag(i) << " "
274  << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
275  << G4endl;
276  }
277  */
278 
279  // expand vectors
280  if(nFlags < nCouples) {
281  for(size_t i=nFlags; i<nCouples; ++i) {
282  theDensityFactor->push_back(1.0);
283  }
284  for(size_t i=nFlags; i<nCouples; ++i) { theDensityIdx->push_back(-1); }
285  for(size_t i=nFlags; i<nCouples; ++i) { theFlag->push_back(true); }
286  }
287  for(size_t i=0; i<nCouples; ++i) {
288 
289  // base material is needed only for a couple which is not
290  // initialised and for which tables will be computed
291  (*theFlag)[i] = table->GetFlag(i);
292  if ((*theDensityIdx)[i] < 0) {
293  (*theDensityIdx)[i] = i;
294  const G4MaterialCutsCouple* couple =
295  theCoupleTable->GetMaterialCutsCouple(i);
296  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
297  const G4Material* mat = couple->GetMaterial();
298  const G4Material* bmat = mat->GetBaseMaterial();
299 
300  // base material exists - find it and check if it can be reused
301  if(bmat) {
302  for(size_t j=0; j<nCouples; ++j) {
303 
304  if(j == i) { continue; }
305  const G4MaterialCutsCouple* bcouple =
306  theCoupleTable->GetMaterialCutsCouple(j);
307 
308  if(bcouple->GetMaterial() == bmat &&
309  bcouple->GetProductionCuts() == pcuts) {
310 
311  // based couple exist in the same region
312  (*theDensityIdx)[i] = j;
313  (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
314  (*theFlag)[i] = false;
315 
316  // ensure that there will no double initialisation
317  (*theDensityIdx)[j] = j;
318  (*theDensityFactor)[j] = 1.0;
319  (*theFlag)[j] = true;
320  break;
321  }
322  }
323  }
324  }
325  }
326  /*
327  for(size_t i=0; i<nCouples; ++i) {
328  G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i]
329  << " TableFlag= " << table->GetFlag(i) << " "
330  << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
331  << G4endl;
332  }
333  G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials end"
334  << G4endl;
335  */
336 }
G4double GetDensity() const
Definition: G4Material.hh:180
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
G4bool GetFlag(size_t i) const
G4ProductionCuts * GetProductionCuts() const
const G4Material * GetMaterial() const

Here is the call graph for this function:

Here is the caller graph for this function:

void G4LossTableBuilder::SetInitialisationFlag ( G4bool  flag)
inline

Definition at line 149 of file G4LossTableBuilder.hh.

150 {
151  isInitialized = flag;
152 }
void G4LossTableBuilder::SetSplineFlag ( G4bool  flag)
inline

Definition at line 144 of file G4LossTableBuilder.hh.

145 {
146  splineFlag = flag;
147 }

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