Geant4  10.01.p02
G4LossTableBuilder.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: G4LossTableBuilder.cc 81058 2014-05-20 09:04:28Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4LossTableBuilder
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 03.01.2002
38 //
39 // Modifications:
40 //
41 // 23-01-03 V.Ivanchenko Cut per region
42 // 21-07-04 V.Ivanchenko Fix problem of range for dedx=0
43 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
44 // 07-12-04 Fix of BuildDEDX table (V.Ivanchenko)
45 // 27-03-06 Add bool options isIonisation (V.Ivanchenko)
46 // 16-01-07 Fill new (not old) DEDX table (V.Ivanchenko)
47 // 12-02-07 Use G4LPhysicsFreeVector for the inverse range table (V.Ivanchenko)
48 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
49 //
50 // Class Description:
51 //
52 // -------------------------------------------------------------------
53 //
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 
57 #include "G4LossTableBuilder.hh"
58 #include "G4SystemOfUnits.hh"
59 #include "G4PhysicsTable.hh"
60 #include "G4PhysicsLogVector.hh"
61 #include "G4PhysicsTableHelper.hh"
62 #include "G4LPhysicsFreeVector.hh"
63 #include "G4ProductionCutsTable.hh"
64 #include "G4MaterialCutsCouple.hh"
65 #include "G4Material.hh"
66 #include "G4VEmModel.hh"
67 #include "G4ParticleDefinition.hh"
68 #include "G4LossTableManager.hh"
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 
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 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
85 {
86  delete theDensityFactor;
87  delete theDensityIdx;
88  delete theFlag;
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
93 void
95  const std::vector<G4PhysicsTable*>& list)
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  /*
113  G4PhysicsLogVector* pv = new G4PhysicsLogVector(pv0->Energy(0),
114  pv0->GetMaxEnergy(),
115  npoints-1);
116  */
117  pv->SetSpline(splineFlag);
118  for (size_t j=0; j<npoints; ++j) {
119  G4double dedx = 0.0;
120  for (size_t k=0; k<n_processes; ++k) {
121  G4PhysicsVector* pv1 = (*(list[k]))[i];
122  dedx += (*pv1)[j];
123  }
124  pv->PutValue(j, dedx);
125  }
126  if(splineFlag) { pv->FillSecondDerivatives(); }
127  G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
128  }
129  }
130 }
131 
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133 
135  G4PhysicsTable* rangeTable,
136  G4bool isIonisation)
137 // Build range table from the energy loss table
138 {
139  size_t nCouples = dedxTable->size();
140  if(0 >= nCouples) { return; }
141 
142  size_t n = 100;
143  G4double del = 1.0/(G4double)n;
144 
145  for (size_t i=0; i<nCouples; ++i) {
146  if(isIonisation) {
147  if( !(*theFlag)[i] ) { continue; }
148  }
149  G4PhysicsLogVector* pv = static_cast<G4PhysicsLogVector*>((*dedxTable)[i]);
150  size_t npoints = pv->GetVectorLength();
151  size_t bin0 = 0;
152  G4double elow = pv->Energy(0);
153  G4double ehigh = pv->Energy(npoints-1);
154  G4double dedx1 = (*pv)[0];
155 
156  //G4cout << "i= " << i << "npoints= " << npoints << " dedx1= "
157  //<< dedx1 << G4endl;
158 
159  // protection for specific cases dedx=0
160  if(dedx1 == 0.0) {
161  for (size_t k=1; k<npoints; ++k) {
162  bin0++;
163  elow = pv->Energy(k);
164  dedx1 = (*pv)[k];
165  if(dedx1 > 0.0) { break; }
166  }
167  npoints -= bin0;
168  }
169  //G4cout<<"New Range vector" << G4endl;
170  //G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh
171  // <<" bin0= " << bin0 <<G4endl;
172 
173  // initialisation of a new vector
174  if(npoints < 2) { npoints = 2; }
175 
176  delete (*rangeTable)[i];
178  if(0 == bin0) { v = new G4PhysicsLogVector(*pv); }
179  else { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); }
180 
181  // dedx is exact zero cannot build range table
182  if(2 == npoints) {
183  v->PutValue(0,1000.);
184  v->PutValue(1,2000.);
185  G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
186  return;
187  }
188  v->SetSpline(splineFlag);
189 
190  // assumed dedx proportional to beta
191  G4double energy1 = v->Energy(0);
192  G4double range = 2.*energy1/dedx1;
193  //G4cout << "range0= " << range << G4endl;
194  v->PutValue(0,range);
195 
196  for (size_t j=1; j<npoints; ++j) {
197 
198  G4double energy2 = v->Energy(j);
199  G4double de = (energy2 - energy1) * del;
200  G4double energy = energy2 + de*0.5;
201  G4double sum = 0.0;
202  //G4cout << "j= " << j << " e1= " << energy1 << " e2= " << energy2
203  // << " n= " << n << G4endl;
204  for (size_t k=0; k<n; ++k) {
205  energy -= de;
206  dedx1 = pv->Value(energy);
207  if(dedx1 > 0.0) { sum += de/dedx1; }
208  }
209  range += sum;
210  v->PutValue(j,range);
211  energy1 = energy2;
212  }
213  if(splineFlag) { v->FillSecondDerivatives(); }
214  G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
215  }
216 }
217 
218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
219 
220 void
222  G4PhysicsTable* invRangeTable,
223  G4bool isIonisation)
224 // Build inverse range table from the energy loss table
225 {
226  size_t nCouples = rangeTable->size();
227  if(0 >= nCouples) { return; }
228 
229  for (size_t i=0; i<nCouples; ++i) {
230 
231  if(isIonisation) {
232  if( !(*theFlag)[i] ) { continue; }
233  }
234  G4PhysicsVector* pv = (*rangeTable)[i];
235  size_t npoints = pv->GetVectorLength();
236  G4double rlow = (*pv)[0];
237  G4double rhigh = (*pv)[npoints-1];
238 
239  delete (*invRangeTable)[i];
240  G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh);
241  v->SetSpline(splineFlag);
242 
243  for (size_t j=0; j<npoints; ++j) {
244  G4double e = pv->Energy(j);
245  G4double r = (*pv)[j];
246  v->PutValues(j,r,e);
247  }
248  if(splineFlag) { v->FillSecondDerivatives(); }
249 
250  G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);
251  }
252 }
253 
254 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255 
256 void
258 {
259  size_t nCouples = table->size();
260  size_t nFlags = theFlag->size();
261 
262  if(nCouples == nFlags && isInitialized) { return; }
263 
264  isInitialized = true;
265 
266  //G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials Ncouples= "
267  // << nCouples << " FlagSize= " << nFlags << G4endl;
268 
269  // variable density check
270  const G4ProductionCutsTable* theCoupleTable=
272 
273  /*
274  for(size_t i=0; i<nFlags; ++i) {
275  G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i]
276  << " tableFlag= " << table->GetFlag(i) << " "
277  << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
278  << G4endl;
279  }
280  */
281 
282  // expand vectors
283  if(nFlags < nCouples) {
284  for(size_t i=nFlags; i<nCouples; ++i) {
285  theDensityFactor->push_back(1.0);
286  }
287  for(size_t i=nFlags; i<nCouples; ++i) { theDensityIdx->push_back(-1); }
288  for(size_t i=nFlags; i<nCouples; ++i) { theFlag->push_back(true); }
289  }
290  for(size_t i=0; i<nCouples; ++i) {
291 
292  // base material is needed only for a couple which is not
293  // initialised and for which tables will be computed
294  (*theFlag)[i] = table->GetFlag(i);
295  if ((*theDensityIdx)[i] < 0) {
296  (*theDensityIdx)[i] = i;
297  const G4MaterialCutsCouple* couple =
298  theCoupleTable->GetMaterialCutsCouple(i);
299  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
300  const G4Material* mat = couple->GetMaterial();
301  const G4Material* bmat = mat->GetBaseMaterial();
302 
303  // base material exists - find it and check if it can be reused
304  if(bmat) {
305  for(size_t j=0; j<nCouples; ++j) {
306 
307  if(j == i) { continue; }
308  const G4MaterialCutsCouple* bcouple =
309  theCoupleTable->GetMaterialCutsCouple(j);
310 
311  if(bcouple->GetMaterial() == bmat &&
312  bcouple->GetProductionCuts() == pcuts) {
313 
314  // based couple exist in the same region
315  (*theDensityIdx)[i] = j;
316  (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
317  (*theFlag)[i] = false;
318 
319  // ensure that there will no double initialisation
320  (*theDensityIdx)[j] = j;
321  (*theDensityFactor)[j] = 1.0;
322  (*theFlag)[j] = true;
323  break;
324  }
325  }
326  }
327  }
328  }
329  /*
330  for(size_t i=0; i<nCouples; ++i) {
331  G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i]
332  << " TableFlag= " << table->GetFlag(i) << " "
333  << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
334  << G4endl;
335  }
336  G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials end"
337  << G4endl;
338  */
339 }
340 
341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
342 
344 {
345  isInitialized = true;
346 
347  // variable density initialisation for the cas without tables
348  const G4ProductionCutsTable* theCoupleTable=
350  size_t nCouples = theCoupleTable->GetTableSize();
351  //G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() nCouples= "
352  // << nCouples << G4endl;
353 
354  theDensityFactor->resize(nCouples, 1.0);
355  theDensityIdx->resize(nCouples, -1);
356  theFlag->resize(nCouples, true);
357 
358  for(size_t i=0; i<nCouples; ++i) {
359 
360  // base material is needed only for a couple which is not
361  // initialised and for which tables will be computed
362  if ((*theDensityIdx)[i] < 0) {
363  (*theDensityIdx)[i] = i;
364  const G4MaterialCutsCouple* couple =
365  theCoupleTable->GetMaterialCutsCouple(i);
366  const G4ProductionCuts* pcuts = couple->GetProductionCuts();
367  const G4Material* mat = couple->GetMaterial();
368  const G4Material* bmat = mat->GetBaseMaterial();
369 
370  // base material exists - find it and check if it can be reused
371  if(bmat) {
372  for(size_t j=0; j<nCouples; ++j) {
373 
374  if(j == i) { continue; }
375  const G4MaterialCutsCouple* bcouple =
376  theCoupleTable->GetMaterialCutsCouple(j);
377 
378  if(bcouple->GetMaterial() == bmat &&
379  bcouple->GetProductionCuts() == pcuts) {
380 
381  // based couple exist in the same region
382  (*theDensityIdx)[i] = j;
383  (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
384  (*theFlag)[i] = false;
385 
386  // ensure that there will no double initialisation
387  (*theDensityIdx)[j] = j;
388  (*theDensityFactor)[j] = 1.0;
389  (*theFlag)[j] = true;
390  break;
391  }
392  }
393  }
394  }
395  }
396  /*
397  for(size_t i=0; i<nCouples; ++i) {
398  G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i] << " "
399  << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
400  << " DensityFactor= " << (*theDensityFactor)[i]
401  << G4endl;
402  }
403  G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() end" << G4endl;
404  */
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408 
411  G4VEmModel* model,
412  const G4ParticleDefinition* part,
413  G4double emin, G4double emax,
414  G4bool spline)
415 {
416  // check input
418  if(!table) { return table; }
419  if(emin >= emax) {
420  table->clearAndDestroy();
421  delete table;
422  table = 0;
423  return table;
424  }
426 
428 
429  // Access to materials
430  const G4ProductionCutsTable* theCoupleTable=
432  size_t numOfCouples = theCoupleTable->GetTableSize();
433 
434  G4PhysicsLogVector* aVector = 0;
435  //G4PhysicsLogVector* bVector = 0;
436 
437  for(size_t i=0; i<numOfCouples; ++i) {
438 
439  //G4cout<< "i= " << i << " Flag= " << GetFlag(i) << G4endl;
440 
441  if (GetFlag(i)) {
442 
443  // create physics vector and fill it
444  const G4MaterialCutsCouple* couple =
445  theCoupleTable->GetMaterialCutsCouple(i);
446  delete (*table)[i];
447 
448  // if start from zero then change the scale
449 
450  const G4Material* mat = couple->GetMaterial();
451 
452  G4double tmin = std::max(emin,model->MinPrimaryEnergy(mat,part));
453  if(0.0 >= tmin) { tmin = eV; }
454  G4int n = nbins;
455 
456  if(tmin >= emax) {
457  aVector = 0;
458  } else {
459  n *= G4int(std::log10(emax/tmin) + 0.5);
460  if(n < 3) { n = 3; }
461  aVector = new G4PhysicsLogVector(tmin, emax, n);
462  }
463 
464  if(aVector) {
465  aVector->SetSpline(spline);
466  for(G4int j=0; j<=n; ++j) {
467  aVector->PutValue(j, model->Value(couple, part,
468  aVector->Energy(j)));
469  }
470  if(spline) { aVector->FillSecondDerivatives(); }
471  }
472  G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
473  }
474  }
475  /*
476  G4cout << "G4LossTableBuilder::BuildTableForModel done for "
477  << part->GetParticleName() << " and "<< model->GetName()
478  << " " << table << G4endl;
479  */
480  return table;
481 }
482 
483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
484 
void BuildRangeTable(const G4PhysicsTable *dedxTable, G4PhysicsTable *rangeTable, G4bool isIonisation=false)
void PutValues(size_t binNumber, G4double binValue, G4double dataValue)
static void SetPhysicsVector(G4PhysicsTable *physTable, size_t idx, G4PhysicsVector *vec)
static G4LossTableManager * Instance()
std::vector< G4int > * theDensityIdx
G4PhysicsTable * BuildTableForModel(G4PhysicsTable *table, G4VEmModel *model, const G4ParticleDefinition *, G4double emin, G4double emax, G4bool spline)
virtual G4double MinPrimaryEnergy(const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
Definition: G4VEmModel.cc:370
G4double GetDensity() const
Definition: G4Material.hh:180
G4bool GetFlag(size_t idx) const
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
void FillSecondDerivatives()
static G4PhysicsTable * PreparePhysicsTable(G4PhysicsTable *physTable)
void SetSpline(G4bool)
bool G4bool
Definition: G4Types.hh:79
void PutValue(size_t index, G4double theValue)
const G4int n
void BuildInverseRangeTable(const G4PhysicsTable *rangeTable, G4PhysicsTable *invRangeTable, G4bool isIonisation=false)
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
void BuildDEDXTable(G4PhysicsTable *dedxTable, const std::vector< G4PhysicsTable * > &)
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:194
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::vector< G4double > * theDensityFactor
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void InitialiseBaseMaterials(G4PhysicsTable *table)
G4int GetNumberOfBinsPerDecade() const
G4double energy(const ThreeVector &p, const G4double m)
std::vector< G4bool > * theFlag
const G4Material * GetBaseMaterial() const
Definition: G4Material.hh:233
G4bool GetFlag(size_t i) const
double G4double
Definition: G4Types.hh:76
G4ProductionCuts * GetProductionCuts() const
void clearAndDestroy()
virtual G4double Value(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
Definition: G4VEmModel.cc:361
const G4Material * GetMaterial() const