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