Geant4  10.02.p03
G4LossTableBuilder Class Reference

#include <G4LossTableBuilder.hh>

Collaboration diagram for G4LossTableBuilder:

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)
 

Private Member Functions

void InitialiseCouples ()
 
G4LossTableBuilderoperator= (const G4LossTableBuilder &right)
 
 G4LossTableBuilder (const G4LossTableBuilder &)
 

Private Attributes

G4bool splineFlag
 
G4bool isInitialized
 
std::vector< G4double > * theDensityFactor
 
std::vector< G4int > * theDensityIdx
 
std::vector< G4bool > * theFlag
 

Detailed Description

Definition at line 61 of file G4LossTableBuilder.hh.

Constructor & Destructor Documentation

◆ G4LossTableBuilder() [1/2]

G4LossTableBuilder::G4LossTableBuilder ( )

Definition at line 72 of file G4LossTableBuilder.cc.

73 {
74  splineFlag = true;
75  isInitialized = false;
76 
77  theDensityFactor = new std::vector<G4double>;
78  theDensityIdx = new std::vector<G4int>;
79  theFlag = new std::vector<G4bool>;
80 }
std::vector< G4int > * theDensityIdx
std::vector< G4double > * theDensityFactor
std::vector< G4bool > * theFlag

◆ ~G4LossTableBuilder()

G4LossTableBuilder::~G4LossTableBuilder ( )
virtual

Definition at line 84 of file G4LossTableBuilder.cc.

85 {
86  delete theDensityFactor;
87  delete theDensityIdx;
88  delete theFlag;
89 }
std::vector< G4int > * theDensityIdx
std::vector< G4double > * theDensityFactor
std::vector< G4bool > * theFlag

◆ G4LossTableBuilder() [2/2]

G4LossTableBuilder::G4LossTableBuilder ( const G4LossTableBuilder )
private

Member Function Documentation

◆ BuildDEDXTable()

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

Definition at line 94 of file G4LossTableBuilder.cc.

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

◆ BuildInverseRangeTable()

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

Definition at line 216 of file G4LossTableBuilder.cc.

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

◆ BuildRangeTable()

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

Definition at line 129 of file G4LossTableBuilder.cc.

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

◆ BuildTableForModel()

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

Definition at line 405 of file G4LossTableBuilder.cc.

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

◆ GetCoupleIndexes()

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

Definition at line 123 of file G4LossTableBuilder.hh.

124 {
125  if(theDensityIdx->size() == 0) { InitialiseCouples(); }
126  return theDensityIdx;
127 }
std::vector< G4int > * theDensityIdx
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetDensityFactors()

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

Definition at line 130 of file G4LossTableBuilder.hh.

131 {
132  if(theDensityIdx->size() == 0) { InitialiseCouples(); }
133  return theDensityFactor;
134 }
std::vector< G4int > * theDensityIdx
std::vector< G4double > * theDensityFactor
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetFlag()

G4bool G4LossTableBuilder::GetFlag ( size_t  idx) const
inline

Definition at line 136 of file G4LossTableBuilder.hh.

137 {
138  return (*theFlag)[idx];
139 }
std::vector< G4bool > * theFlag
Here is the caller graph for this function:

◆ InitialiseBaseMaterials()

void G4LossTableBuilder::InitialiseBaseMaterials ( G4PhysicsTable table)

Definition at line 252 of file G4LossTableBuilder.cc.

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

◆ InitialiseCouples()

void G4LossTableBuilder::InitialiseCouples ( )
private

Definition at line 338 of file G4LossTableBuilder.cc.

339 {
340  isInitialized = true;
341 
342  // variable density initialisation for the cas without tables
343  const G4ProductionCutsTable* theCoupleTable=
345  size_t nCouples = theCoupleTable->GetTableSize();
346  //G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() nCouples= "
347  // << nCouples << G4endl;
348 
349  theDensityFactor->resize(nCouples, 1.0);
350  theDensityIdx->resize(nCouples, -1);
351  theFlag->resize(nCouples, true);
352 
353  for(size_t i=0; i<nCouples; ++i) {
354 
355  // base material is needed only for a couple which is not
356  // initialised and for which tables will be computed
357  if ((*theDensityIdx)[i] < 0) {
358  (*theDensityIdx)[i] = i;
359  const G4MaterialCutsCouple* couple =
360  theCoupleTable->GetMaterialCutsCouple(i);
361  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
362  const G4Material* mat = couple->GetMaterial();
363  const G4Material* bmat = mat->GetBaseMaterial();
364 
365  // base material exists - find it and check if it can be reused
366  if(bmat) {
367  for(size_t j=0; j<nCouples; ++j) {
368 
369  if(j == i) { continue; }
370  const G4MaterialCutsCouple* bcouple =
371  theCoupleTable->GetMaterialCutsCouple(j);
372 
373  if(bcouple->GetMaterial() == bmat &&
374  bcouple->GetProductionCuts() == pcuts) {
375 
376  // based couple exist in the same region
377  (*theDensityIdx)[i] = j;
378  (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
379  (*theFlag)[i] = false;
380 
381  // ensure that there will no double initialisation
382  (*theDensityIdx)[j] = j;
383  (*theDensityFactor)[j] = 1.0;
384  (*theFlag)[j] = true;
385  break;
386  }
387  }
388  }
389  }
390  }
391  /*
392  for(size_t i=0; i<nCouples; ++i) {
393  G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i] << " "
394  << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
395  << " DensityFactor= " << (*theDensityFactor)[i]
396  << G4endl;
397  }
398  G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() end" << G4endl;
399  */
400 }
const G4Material * GetMaterial() const
std::vector< G4int > * theDensityIdx
G4ProductionCuts * GetProductionCuts() const
G4double GetDensity() const
Definition: G4Material.hh:180
Float_t mat
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::vector< G4double > * theDensityFactor
std::vector< G4bool > * theFlag
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
Here is the call graph for this function:
Here is the caller graph for this function:

◆ operator=()

G4LossTableBuilder& G4LossTableBuilder::operator= ( const G4LossTableBuilder right)
private

◆ SetInitialisationFlag()

void G4LossTableBuilder::SetInitialisationFlag ( G4bool  flag)
inline

Definition at line 146 of file G4LossTableBuilder.hh.

147 {
148  isInitialized = flag;
149 }
Here is the caller graph for this function:

◆ SetSplineFlag()

void G4LossTableBuilder::SetSplineFlag ( G4bool  flag)
inline

Definition at line 141 of file G4LossTableBuilder.hh.

142 {
143  splineFlag = flag;
144 }
Here is the caller graph for this function:

Member Data Documentation

◆ isInitialized

G4bool G4LossTableBuilder::isInitialized
private

Definition at line 114 of file G4LossTableBuilder.hh.

◆ splineFlag

G4bool G4LossTableBuilder::splineFlag
private

Definition at line 113 of file G4LossTableBuilder.hh.

◆ theDensityFactor

std::vector<G4double>* G4LossTableBuilder::theDensityFactor
private

Definition at line 116 of file G4LossTableBuilder.hh.

◆ theDensityIdx

std::vector<G4int>* G4LossTableBuilder::theDensityIdx
private

Definition at line 117 of file G4LossTableBuilder.hh.

◆ theFlag

std::vector<G4bool>* G4LossTableBuilder::theFlag
private

Definition at line 118 of file G4LossTableBuilder.hh.


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