Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VRangeToEnergyConverter.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 //
27 // $Id: G4VRangeToEnergyConverter.cc 93090 2015-10-05 13:14:43Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------
31 // GEANT 4 class implementation file/ History:
32 // 5 Oct. 2002, H.Kuirashige : Structure created based on object model
33 // --------------------------------------------------------------
34 
36 #include "G4ParticleTable.hh"
37 #include "G4Material.hh"
38 #include "G4MaterialTable.hh"
39 #include "G4PhysicsLogVector.hh"
40 
41 #include "G4ios.hh"
42 #include "G4SystemOfUnits.hh"
43 
47 
49  theParticle(0), theLossTable(0), NumberOfElements(0), TotBin(300),
50  verboseLevel(1)
51 {
52  fMaxEnergyCut = 0.;
53 }
54 
55 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter(const G4VRangeToEnergyConverter& right) : theParticle(right.theParticle), theLossTable(0), TotBin(right.TotBin)
56 {
57  fMaxEnergyCut = 0.;
58  if (theLossTable) {
60  delete theLossTable;
61  theLossTable=0;
62  }
63  *this = right;
64 }
65 
67 {
68  if (this == &right) return *this;
69  if (theLossTable) {
71  delete theLossTable;
72  theLossTable=0;
73  }
74 
77  theParticle = right.theParticle;
78  verboseLevel = right.verboseLevel;
79 
80  // create the loss table
81  theLossTable = new G4LossTable();
83  // fill the loss table
84  for (size_t j=0; j<size_t(NumberOfElements); j++){
86  for (size_t i=0; i<=size_t(TotBin); i++) {
87  G4double Value = (*((*right.theLossTable)[j]))[i];
88  aVector->PutValue(i,Value);
89  }
90  theLossTable->insert(aVector);
91  }
92 
93  // clean up range vector store
94  for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
95  delete fRangeVectorStore.at(idx);
96  }
97  fRangeVectorStore.clear();
98 
99  // copy range vector store
100  for (size_t j=0; j<((right.fRangeVectorStore).size()); j++){
101  G4RangeVector* vector = (right.fRangeVectorStore).at(j);
102  G4RangeVector* rangeVector = 0;
103  if (vector !=0 ) {
104  rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin);
106  for (size_t i=0; i<=size_t(TotBin); i++) {
107  G4double Value = (*vector)[i];
108  rangeVector->PutValue(i,Value);
109  }
110  }
111  fRangeVectorStore.push_back(rangeVector);
112  }
113  return *this;
114 }
115 
116 
118 {
119  Reset();
120  // Comment out Reset() for MT application
121 
131 
132 }
133 
135 {
136  return this == &right;
137 }
138 
140 {
141  return this != &right;
142 }
143 
144 
145 // **********************************************************************
146 // ************************* Convert ***********************************
147 // **********************************************************************
149  const G4Material* material)
150 {
151 #ifdef G4VERBOSE
152  if (GetVerboseLevel()>3) {
153  G4cout << "G4VRangeToEnergyConverter::Convert() ";
154  G4cout << "Convert for " << material->GetName()
155  << " with Range Cut " << rangeCut/mm << "[mm]" << G4endl;
156  }
157 #endif
158 
159  G4double theKineticEnergyCuts = 0.;
160 
161  if (fMaxEnergyCut != MaxEnergyCut) {
163  // clear loss table and renge vectors
164  Reset();
165  }
166 
167  // Build the energy loss table
168  BuildLossTable();
169 
170  // Build range vector for every material, convert cut into energy-cut,
171  // fill theKineticEnergyCuts and delete the range vector
172  static const G4double tune = 0.025*mm*g/cm3 ,lowen = 30.*keV ;
173 
174  // check density
175  G4double density = material->GetDensity() ;
176  if(density <= 0.) {
177 #ifdef G4VERBOSE
178  if (GetVerboseLevel()>0) {
179  G4cout << "G4VRangeToEnergyConverter::Convert() ";
180  G4cout << material->GetName() << "has zero density "
181  << "( " << density << ")" << G4endl;
182  }
183 #endif
184  return 0.;
185  }
186 
187  // initialize RangeVectorStore
189  G4int ext_size = table->size() - fRangeVectorStore.size();
190  for (int i=0; i<ext_size; i++) fRangeVectorStore.push_back(0);
191 
192  // Build Range Vector
193  G4int idx = material->GetIndex();
194  G4RangeVector* rangeVector = fRangeVectorStore.at(idx);
195  if (rangeVector == 0) {
196  rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin);
197  BuildRangeVector(material, rangeVector);
198  fRangeVectorStore.at(idx) = rangeVector;
199  }
200 
201  // Convert Range Cut ro Kinetic Energy Cut
202  theKineticEnergyCuts = ConvertCutToKineticEnergy(rangeVector, rangeCut, idx);
203 
204  if( ((theParticle->GetParticleName()=="e-")||(theParticle->GetParticleName()=="e+"))
205  && (theKineticEnergyCuts < lowen) ) {
206  // corr. should be switched on smoothly
207  theKineticEnergyCuts /= (1.+(1.-theKineticEnergyCuts/lowen)*
208  tune/(rangeCut*density));
209  }
210 
211  if(theKineticEnergyCuts < LowestEnergy) {
212  theKineticEnergyCuts = LowestEnergy ;
213  } else if(theKineticEnergyCuts > MaxEnergyCut) {
214  theKineticEnergyCuts = MaxEnergyCut;
215  }
216 
217  return theKineticEnergyCuts;
218 }
219 
220 // **********************************************************************
221 // ************************ SetEnergyRange *****************************
222 // **********************************************************************
224  G4double highedge)
225 {
226  // check LowestEnergy/ HighestEnergy
227  if ( (lowedge<0.0)||(highedge<=lowedge) ){
228 #ifdef G4VERBOSE
229  G4cerr << "Error in G4VRangeToEnergyConverter::SetEnergyRange";
230  G4cerr << " : illegal energy range" << "(" << lowedge/GeV;
231  G4cerr << "," << highedge/GeV << ") [GeV]" << G4endl;
232 #endif
233  G4Exception( "G4VRangeToEnergyConverter::SetEnergyRange()",
234  "ProcCuts101",
235  JustWarning, "Illegal energy range ");
236  } else {
237  LowestEnergy = lowedge;
238  HighestEnergy = highedge;
239  }
240 }
241 
242 
244 {
245  return LowestEnergy;
246 }
247 
248 
250 {
251  return HighestEnergy;
252 }
253 
254 // **********************************************************************
255 // ******************* Get/SetMaxEnergyCut *****************************
256 // **********************************************************************
258 {
259  return MaxEnergyCut;
260 }
261 
263 {
265 }
266 
267 // **********************************************************************
268 // ************************ Reset **************************************
269 // **********************************************************************
271 {
272  // delete loss table
273  if (theLossTable) {
275  delete theLossTable;
276  }
277  theLossTable=0;
279 
280  //clear RangeVectorStore
281  for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
282  delete fRangeVectorStore.at(idx);
283  }
284  fRangeVectorStore.clear();
285 }
286 
287 
288 // **********************************************************************
289 // ************************ BuildLossTable ******************************
290 // **********************************************************************
291 // create Energy Loss Table for charged particles
292 // (cross section tabel for neutral )
294 {
295  if (size_t(NumberOfElements) == G4Element::GetNumberOfElements()) return;
296 
297  // clear Loss table and Range vectors
298  Reset();
299 
300  // Build dE/dx tables for elements
302  theLossTable = new G4LossTable();
304 #ifdef G4VERBOSE
305  if (GetVerboseLevel()>3) {
306  G4cout << "G4VRangeToEnergyConverter::BuildLossTable() ";
307  G4cout << "Create theLossTable[" << theLossTable << "]";
308  G4cout << " NumberOfElements=" << NumberOfElements <<G4endl;
309  }
310 #endif
311 
312 
313  // fill the loss table
314  for (size_t j=0; j<size_t(NumberOfElements); j++){
315  G4double Value;
316  G4LossVector* aVector= 0;
318  for (size_t i=0; i<=size_t(TotBin); i++) {
319  Value = ComputeLoss( (*G4Element::GetElementTable())[j]->GetZ(),
320  aVector->Energy(i)
321  );
322  aVector->PutValue(i,Value);
323  }
324  theLossTable->insert(aVector);
325  }
326 }
327 
328 // **********************************************************************
329 // ************************ BuildRangeVector ****************************
330 // **********************************************************************
332  G4PhysicsLogVector* rangeVector)
333 {
334  // create range vector for a material
335  const G4ElementVector* elementVector = aMaterial->GetElementVector();
336  const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
337  G4int NumEl = aMaterial->GetNumberOfElements();
338 
339  // calculate parameters of the low energy part first
340  size_t i;
341  std::vector<G4double> lossV;
342  for ( size_t ib=0; ib<=size_t(TotBin); ib++) {
343  G4double loss=0.;
344  for (i=0; i<size_t(NumEl); i++) {
345  G4int IndEl = (*elementVector)[i]->GetIndex();
346  loss += atomicNumDensityVector[i]*
347  (*((*theLossTable)[IndEl]))[ib];
348  }
349  lossV.push_back(loss);
350  }
351 
352  // Integrate with Simpson formula with logarithmic binning
353  G4double dltau = 1.0;
354  if (LowestEnergy>0.) {
355  G4double ltt =std::log(MaxEnergyCut/LowestEnergy);
356  dltau = ltt/TotBin;
357  }
358 
359  G4double s0 = 0.;
360  G4double Value;
361  for ( i=0; i<=size_t(TotBin); i++) {
362  G4double t = rangeVector->GetLowEdgeEnergy(i);
363  G4double q = t/lossV[i];
364  if (i==0) s0 += 0.5*q;
365  else s0 += q;
366 
367  if (i==0) {
368  Value = (s0 + 0.5*q)*dltau ;
369  } else {
370  Value = (s0 - 0.5*q)*dltau ;
371  }
372  rangeVector->PutValue(i,Value);
373  }
374 }
375 
376 // **********************************************************************
377 // ****************** ConvertCutToKineticEnergy *************************
378 // **********************************************************************
380  G4RangeVector* rangeVector,
381  G4double theCutInLength,
382 #ifdef G4VERBOSE
383  size_t materialIndex
384 #else
385  size_t
386 #endif
387  ) const
388 {
389  const G4double epsilon=0.01;
390 
391  // find max. range and the corresponding energy (rmax,Tmax)
392  G4double rmax= -1.e10*mm;
393 
394  G4double T1 = LowestEnergy;
395  G4double r1 =(*rangeVector)[0] ;
396 
397  G4double T2 = MaxEnergyCut;
398 
399  // check theCutInLength < r1
400  if ( theCutInLength <= r1 ) { return T1; }
401 
402  // scan range vector to find nearest bin
403  // ( suppose that r(Ti) > r(Tj) if Ti >Tj )
404  for (size_t ibin=0; ibin<=size_t(TotBin); ibin++) {
405  G4double T=rangeVector->GetLowEdgeEnergy(ibin);
406  G4double r=(*rangeVector)[ibin];
407  if ( r>rmax ) rmax=r;
408  if (r <theCutInLength ) {
409  T1 = T;
410  r1 = r;
411  } else if (r >theCutInLength ) {
412  T2 = T;
413  break;
414  }
415  }
416 
417  // check cut in length is smaller than range max
418  if ( theCutInLength >= rmax ) {
419 #ifdef G4VERBOSE
420  if (GetVerboseLevel()>2) {
421  G4cout << "G4VRangeToEnergyConverter::ConvertCutToKineticEnergy ";
422  G4cout << " for " << theParticle->GetParticleName() << G4endl;
423  G4cout << "The cut in range [" << theCutInLength/mm << " (mm)] ";
424  G4cout << " is too big " ;
425  G4cout << " for material idx=" << materialIndex <<G4endl;
426  }
427 #endif
428  return MaxEnergyCut;
429  }
430 
431  // convert range to energy
432  G4double T3 = std::sqrt(T1*T2);
433  G4double r3 = rangeVector->Value(T3);
434  const size_t MAX_LOOP = 1000;
435  for (size_t loop_count=0; loop_count<MAX_LOOP; ++loop_count){
436  if (std::fabs(1.-r3/theCutInLength)<epsilon ) break;
437  if ( theCutInLength <= r3 ) {
438  T2 = T3;
439  } else {
440  T1 = T3;
441  }
442  T3 = std::sqrt(T1*T2);
443  r3 = rangeVector->Value(T3);
444  }
445 
446  return T3;
447 }
448 
static constexpr double mm
Definition: G4SIunits.hh:115
G4int operator!=(const G4VRangeToEnergyConverter &right) const
std::vector< G4Element * > G4ElementVector
void insert(G4PhysicsVector *)
size_t GetIndex() const
Definition: G4Material.hh:262
const G4String & GetName() const
Definition: G4Material.hh:178
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
G4double GetDensity() const
Definition: G4Material.hh:180
virtual void BuildRangeVector(const G4Material *aMaterial, G4RangeVector *rangeVector)
G4double GetLowEdgeEnergy(size_t binNumber) const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
string material
Definition: eplot.py:19
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
virtual G4double ComputeLoss(G4double AtomicNumber, G4double KineticEnergy)=0
G4GLOB_DLL std::ostream G4cout
G4double ConvertCutToKineticEnergy(G4RangeVector *theRangeVector, G4double theCutInLength, size_t materialIndex) const
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
virtual G4double Convert(G4double rangeCut, const G4Material *material)
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4int operator==(const G4VRangeToEnergyConverter &right) const
static void SetMaxEnergyCut(G4double value)
void PutValue(size_t index, G4double theValue)
G4VRangeToEnergyConverter & operator=(const G4VRangeToEnergyConverter &right)
static void SetEnergyRange(G4double lowedge, G4double highedge)
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:216
static constexpr double cm3
Definition: G4SIunits.hh:121
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::vector< G4RangeVector * > fRangeVectorStore
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
static constexpr double keV
Definition: G4SIunits.hh:216
double epsilon(double density, double temperature)
void clearAndDestroy()
const G4ParticleDefinition * theParticle
G4GLOB_DLL std::ostream G4cerr