source: trunk/source/processes/cuts/src/G4VRangeToEnergyConverter.cc@ 968

Last change on this file since 968 was 961, checked in by garnier, 17 years ago

update processes
File size: 16.3 KB
Line 
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,v 1.9 2008/03/02 10:52:56 kurasige Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31// --------------------------------------------------------------
32// GEANT 4 class implementation file/ History:
33// 5 Oct. 2002, H.Kuirashige : Structure created based on object model
34// --------------------------------------------------------------
35
36#include "G4VRangeToEnergyConverter.hh"
37#include "G4ParticleTable.hh"
38#include "G4Material.hh"
39#include "G4PhysicsLogVector.hh"
40
41#include "G4ios.hh"
42
43// energy range
44G4double G4VRangeToEnergyConverter::LowestEnergy = 0.99e-3*MeV;
45G4double G4VRangeToEnergyConverter::HighestEnergy = 100.0e6*MeV;
46
47G4VRangeToEnergyConverter::G4VRangeToEnergyConverter():
48 theParticle(0), theLossTable(0), NumberOfElements(0), TotBin(200),
49 verboseLevel(1)
50{
51}
52
53G4VRangeToEnergyConverter::G4VRangeToEnergyConverter(const G4VRangeToEnergyConverter& right)
54{
55 *this = right;
56}
57
58G4VRangeToEnergyConverter & G4VRangeToEnergyConverter::operator=(const G4VRangeToEnergyConverter &right)
59{
60 if (this == &right) return *this;
61 if (theLossTable) {
62 theLossTable->clearAndDestroy();
63 delete theLossTable;
64 theLossTable=0;
65 }
66
67 NumberOfElements = right.NumberOfElements;
68 TotBin = right.TotBin;
69 theParticle = right.theParticle;
70 verboseLevel = right.verboseLevel;
71
72 // create the loss table
73 theLossTable = new G4LossTable();
74 theLossTable->reserve(G4Element::GetNumberOfElements());
75 // fill the loss table
76 for (size_t j=0; j<size_t(NumberOfElements); j++){
77 G4LossVector* aVector= new
78 G4LossVector(LowestEnergy, HighestEnergy, TotBin);
79 for (size_t i=0; i<size_t(TotBin); i++) {
80 G4double Value = (*((*right.theLossTable)[j]))[i];
81 aVector->PutValue(i,Value);
82 }
83 theLossTable->insert(aVector);
84 }
85 return *this;
86}
87
88
89G4VRangeToEnergyConverter::~G4VRangeToEnergyConverter()
90{
91 if (theLossTable) {
92 theLossTable->clearAndDestroy();
93 delete theLossTable;
94 }
95 theLossTable=0;
96}
97
98G4int G4VRangeToEnergyConverter::operator==(const G4VRangeToEnergyConverter &right) const
99{
100 return this == &right;
101}
102
103G4int G4VRangeToEnergyConverter::operator!=(const G4VRangeToEnergyConverter &right) const
104{
105 return this != &right;
106}
107
108
109// **********************************************************************
110// ************************* Convert ***********************************
111// **********************************************************************
112G4double G4VRangeToEnergyConverter::Convert(G4double rangeCut,
113 const G4Material* material)
114{
115 G4double Mass = theParticle->GetPDGMass();
116 G4double theKineticEnergyCuts = 0.;
117
118 // Build the energy loss table
119 BuildLossTable();
120
121 // Build range vector for every material, convert cut into energy-cut,
122 // fill theKineticEnergyCuts and delete the range vector
123 G4double tune = 0.025*mm*g/cm3 ,lowen = 30.*keV ;
124
125 G4int idx = material->GetIndex();
126 G4double density = material->GetDensity() ;
127 if(density > 0.) {
128 G4RangeVector* rangeVector = new G4RangeVector(LowestEnergy, HighestEnergy, TotBin);
129 BuildRangeVector(material, HighestEnergy, Mass, rangeVector);
130 theKineticEnergyCuts = ConvertCutToKineticEnergy(rangeVector, rangeCut, idx);
131
132 if( ((theParticle->GetParticleName()=="e-")||(theParticle->GetParticleName()=="e+"))
133 && (theKineticEnergyCuts < lowen) )
134
135 // corr. should be switched on smoothly
136 { theKineticEnergyCuts /= (1.+(1.-theKineticEnergyCuts/lowen)*
137 tune/(rangeCut*density)); }
138 if(theKineticEnergyCuts < LowestEnergy) {
139 theKineticEnergyCuts = LowestEnergy ;
140 }
141 delete rangeVector;
142 }
143
144 return theKineticEnergyCuts;
145}
146
147// **********************************************************************
148// ************************ SetEnergyRange *****************************
149// **********************************************************************
150void G4VRangeToEnergyConverter::SetEnergyRange(G4double lowedge,
151 G4double highedge)
152{
153 // check LowestEnergy/ HighestEnergy
154 if ( (lowedge<0.0)||(highedge<=lowedge) ){
155 G4cerr << "Error in G4VRangeToEnergyConverter::SetEnergyRange";
156 G4cerr << " : illegal energy range" << "(" << lowedge/GeV;
157 G4cerr << "," << highedge/GeV << ") [GeV]" << G4endl;
158 } else {
159 LowestEnergy = lowedge;
160 HighestEnergy = highedge;
161 }
162}
163
164
165G4double G4VRangeToEnergyConverter::GetLowEdgeEnergy()
166{
167 return LowestEnergy;
168}
169
170
171G4double G4VRangeToEnergyConverter::GetHighEdgeEnergy()
172{
173 return HighestEnergy;
174}
175
176// **********************************************************************
177// ************************ RangeLinSimpson *****************************
178// **********************************************************************
179G4double G4VRangeToEnergyConverter::RangeLinSimpson(
180 G4int numberOfElement,
181 const G4ElementVector* elementVector,
182 const G4double* atomicNumDensityVector,
183 G4double aMass,
184 G4double taulow, G4double tauhigh, G4int nbin)
185{
186 // Simpson numerical integration, linear binning
187 G4double dtau = (tauhigh-taulow)/nbin;
188 G4double Value=0.;
189 for (size_t i=0; i<=size_t(nbin); i++){
190 G4double taui=taulow+dtau*i;
191 G4double ti=aMass*taui;
192 G4double lossi=0.;
193 size_t nEl = (size_t)(numberOfElement);
194 for (size_t j=0; j<nEl; j++) {
195 G4bool isOut;
196 G4int IndEl = (*elementVector)[j]->GetIndex();
197 lossi += atomicNumDensityVector[j]*
198 (*theLossTable)[IndEl]->GetValue(ti,isOut);
199 }
200 if ( i==0 ) {
201 Value += 0.5/lossi;
202 } else {
203 if ( i<size_t(nbin) ) Value += 1./lossi;
204 else Value += 0.5/lossi;
205 }
206 }
207 Value *= aMass*dtau;
208
209 return Value;
210}
211
212
213// **********************************************************************
214// ************************ RangeLogSimpson *****************************
215// **********************************************************************
216G4double G4VRangeToEnergyConverter::RangeLogSimpson(
217 G4int numberOfElement,
218 const G4ElementVector* elementVector,
219 const G4double* atomicNumDensityVector,
220 G4double aMass,
221 G4double ltaulow, G4double ltauhigh,
222 G4int nbin)
223{
224 // Simpson numerical integration, logarithmic binning
225 if(nbin<0) nbin = TotBin;
226 G4double ltt = ltauhigh-ltaulow;
227 G4double dltau = ltt/nbin;
228 G4double Value = 0.;
229 for (size_t i=0; i<=size_t(nbin); i++){
230 G4double ui = ltaulow+dltau*i;
231 G4double taui = std::exp(ui);
232 G4double ti = aMass*taui;
233 G4double lossi = 0.;
234 size_t nEl = (size_t)(numberOfElement);
235
236 for (size_t j=0; j<nEl; j++) {
237 G4bool isOut;
238 G4int IndEl = (*elementVector)[j]->GetIndex();
239 lossi += atomicNumDensityVector[j]*
240 (*theLossTable)[IndEl]->GetValue(ti,isOut);
241 }
242 if ( i==0 ) {
243 Value += 0.5*taui/lossi;
244 } else {
245 if ( i<size_t(nbin) ) Value += taui/lossi;
246 else Value += 0.5*taui/lossi;
247 }
248 }
249 Value *= aMass*dltau;
250
251 return Value;
252}
253
254// **********************************************************************
255// ************************ BuildLossTable ******************************
256// **********************************************************************
257// create Energy Loss Table for charged particles
258// (cross section tabel for neutral )
259void G4VRangeToEnergyConverter::BuildLossTable()
260{
261 // Build dE/dx tables for elements
262 if (size_t(NumberOfElements) != G4Element::GetNumberOfElements()) {
263 if (theLossTable!=0) {
264 theLossTable->clearAndDestroy();
265 delete theLossTable;
266 }
267 theLossTable =0;
268 NumberOfElements = 0;
269 }
270
271 if (NumberOfElements ==0) {
272 NumberOfElements = G4Element::GetNumberOfElements();
273 theLossTable = new G4LossTable();
274 theLossTable->reserve(G4Element::GetNumberOfElements());
275#ifdef G4VERBOSE
276 if (GetVerboseLevel()>3) {
277 G4cout << "G4VRangeToEnergyConverter::BuildLossTable() ";
278 G4cout << "Create theLossTable[" << theLossTable << "]";
279 G4cout << " NumberOfElements=" << NumberOfElements <<G4endl;
280 }
281#endif
282
283
284 // fill the loss table
285 for (size_t j=0; j<size_t(NumberOfElements); j++){
286 G4double Value;
287 G4LossVector* aVector= new
288 G4LossVector(LowestEnergy, HighestEnergy, TotBin);
289 for (size_t i=0; i<size_t(TotBin); i++) {
290 Value = ComputeLoss( (*G4Element::GetElementTable())[j]->GetZ(),
291 aVector->GetLowEdgeEnergy(i)
292 );
293 aVector->PutValue(i,Value);
294 }
295 theLossTable->insert(aVector);
296 }
297 }
298}
299
300// **********************************************************************
301// ************************** ComputeLoss *******************************
302// **********************************************************************
303G4double G4VRangeToEnergyConverter::ComputeLoss(G4double AtomicNumber,
304 G4double KineticEnergy) const
305{
306 // calculate dE/dx
307
308 static G4double Z;
309 static G4double ionpot, tau0, taum, taul, ca, cba, cc;
310
311 G4double z2Particle = theParticle->GetPDGCharge()/eplus;
312 z2Particle *= z2Particle;
313 if (z2Particle < 0.1) return 0.0;
314
315 if( std::abs(AtomicNumber-Z)>0.1 ){
316 // recalculate constants
317 Z = AtomicNumber;
318 G4double Z13 = std::exp(std::log(Z)/3.);
319 tau0 = 0.1*Z13*MeV/proton_mass_c2;
320 taum = 0.035*Z13*MeV/proton_mass_c2;
321 taul = 2.*MeV/proton_mass_c2;
322 ionpot = 1.6e-5*MeV*std::exp(0.9*std::log(Z));
323 cc = (taul+1.)*(taul+1.)*std::log(2.*electron_mass_c2*taul*(taul+2.)/ionpot)/(taul*(taul+2.))-1.;
324 cc = 2.*twopi_mc2_rcl2*Z*cc*std::sqrt(taul);
325 ca = cc/((1.-0.5*std::sqrt(tau0/taum))*tau0);
326 cba = -0.5/std::sqrt(taum);
327 }
328
329 G4double tau = KineticEnergy/theParticle->GetPDGMass();
330 G4double dEdx;
331 if ( tau <= tau0 ) {
332 dEdx = ca*(std::sqrt(tau)+cba*tau);
333 } else {
334 if( tau <= taul ) {
335 dEdx = cc/std::sqrt(tau);
336 } else {
337 dEdx = (tau+1.)*(tau+1.)*
338 std::log(2.*electron_mass_c2*tau*(tau+2.)/ionpot)/(tau*(tau+2.))-1.;
339 dEdx = 2.*twopi_mc2_rcl2*Z*dEdx;
340 }
341 }
342 return dEdx*z2Particle ;
343}
344
345// **********************************************************************
346// ************************ BuildRangeVector ****************************
347// **********************************************************************
348void G4VRangeToEnergyConverter::BuildRangeVector(
349 const G4Material* aMaterial,
350 G4double maxEnergy,
351 G4double aMass,
352 G4RangeVector* rangeVector)
353{
354 // create range vector for a material
355 const G4double tlim=2.*MeV, t1=0.1*MeV, t2=0.025*MeV;
356 const G4int maxnbint=100;
357
358 const G4ElementVector* elementVector = aMaterial->GetElementVector();
359 const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
360
361 G4int NumEl = aMaterial->GetNumberOfElements();
362
363 // calculate parameters of the low energy part first
364 G4double loss1=0.;
365 G4double loss2=0.;
366 size_t i;
367 for (i=0; i<size_t(NumEl); i++) {
368 G4bool isOut;
369 G4int IndEl = (*elementVector)[i]->GetIndex();
370 loss1 += atomicNumDensityVector[i]*
371 (*theLossTable)[IndEl]->GetValue(t1,isOut);
372 loss2 += atomicNumDensityVector[i]*
373 (*theLossTable)[IndEl]->GetValue(t2,isOut);
374 }
375 G4double tau1 = t1/proton_mass_c2;
376 G4double sqtau1 = std::sqrt(tau1);
377 G4double ca = (4.*loss2-loss1)/sqtau1;
378 G4double cb = (2.*loss1-4.*loss2)/tau1;
379 G4double cba = cb/ca;
380 G4double taulim = tlim/proton_mass_c2;
381 G4double taumax = maxEnergy/aMass;
382 G4double ltaumax = std::log(taumax);
383
384 // now we can fill the range vector....
385 G4double rmax = 0.0;
386 for (i=0; i<size_t(TotBin); i++) {
387 G4double LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i);
388 G4double tau = LowEdgeEnergy/aMass;
389 G4double Value;
390
391 if ( tau <= tau1 ){
392 Value =2.*aMass*std::log(1.+cba*std::sqrt(tau))/cb;
393 } else {
394 Value = 2.*aMass*std::log(1.+cba*sqtau1)/cb;
395 if ( tau <= taulim ) {
396 G4int nbin = (G4int)(maxnbint*(tau-tau1)/(taulim-tau1));
397 if ( nbin<1 ) nbin = 1;
398 Value += RangeLinSimpson( NumEl, elementVector,
399 atomicNumDensityVector, aMass,
400 tau1, tau, nbin);
401 } else {
402 Value += RangeLinSimpson( NumEl, elementVector,
403 atomicNumDensityVector, aMass,
404 tau1, taulim, maxnbint);
405 G4double ltaulow = std::log(taulim);
406 G4double ltauhigh = std::log(tau);
407 G4int nbin = (G4int)(maxnbint*(ltauhigh-ltaulow)/(ltaumax-ltaulow));
408 if ( nbin<1 ) nbin = 1;
409 Value += RangeLogSimpson(NumEl, elementVector,
410 atomicNumDensityVector, aMass,
411 ltaulow, ltauhigh, nbin);
412 }
413 }
414 rangeVector->PutValue(i,Value);
415 if (rmax < Value) rmax = Value;
416 }
417}
418
419// **********************************************************************
420// ****************** ConvertCutToKineticEnergy *************************
421// **********************************************************************
422G4double G4VRangeToEnergyConverter::ConvertCutToKineticEnergy(
423 G4RangeVector* rangeVector,
424 G4double theCutInLength,
425 size_t materialIndex
426 ) const
427{
428 const G4double epsilon=0.01;
429
430 // find max. range and the corresponding energy (rmax,Tmax)
431 G4double rmax= -1.e10*mm;
432 G4double Tmax= HighestEnergy;
433 G4double fac = std::exp( std::log(HighestEnergy/LowestEnergy)/TotBin );
434 G4double T=LowestEnergy/fac;
435 G4bool isOut;
436
437 for (size_t ibin=0; ibin<size_t(TotBin); ibin++) {
438 T *= fac;
439 G4double r=rangeVector->GetValue(T,isOut);
440 if ( r>rmax ) {
441 Tmax=T;
442 rmax=r;
443 }
444 }
445
446 // check cut in length is smaller than range max
447 if ( theCutInLength >= rmax ) {
448#ifdef G4VERBOSE
449 if (GetVerboseLevel()>2) {
450 G4cout << "G4VRangeToEnergyConverter::ConvertCutToKineticEnergy ";
451 G4cout << " for " << theParticle->GetParticleName() << G4endl;
452 G4cout << "The cut in range [" << theCutInLength/mm << " (mm)] ";
453 G4cout << " is too big " ;
454 G4cout << " for material idx=" << materialIndex <<G4endl;
455 G4cout << "The cut in energy is set" << DBL_MAX/GeV << "GeV " <<G4endl;
456 }
457#endif
458 return DBL_MAX;
459 }
460
461 // convert range to energy
462 G4double T1 = LowestEnergy;
463 G4double r1 = rangeVector->GetValue(T1,isOut);
464 if ( theCutInLength <= r1 )
465 {
466 return T1;
467 }
468
469 G4double T2 = Tmax ;
470 G4double T3 = std::sqrt(T1*T2);
471 G4double r3 = rangeVector->GetValue(T3,isOut);
472 while ( std::abs(1.-r3/theCutInLength)>epsilon ) {
473 if ( theCutInLength <= r3 ) {
474 T2 = T3;
475 } else {
476 T1 = T3;
477 }
478 T3 = std::sqrt(T1*T2);
479 r3 = rangeVector->GetValue(T3,isOut);
480 }
481
482 return T3;
483}
484
Note: See TracBrowser for help on using the repository browser.