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

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 13.4 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.15 2009/09/14 07:27:46 kurasige Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
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 "G4MaterialTable.hh"
40#include "G4PhysicsLogVector.hh"
41
42#include "G4ios.hh"
43
44// energy range
45G4double G4VRangeToEnergyConverter::LowestEnergy = 0.99e-3*MeV;
46G4double G4VRangeToEnergyConverter::HighestEnergy = 100.0e6*MeV;
47
48// max energy cut
49G4double G4VRangeToEnergyConverter::MaxEnergyCut = 10.0*GeV;
50
51G4VRangeToEnergyConverter::G4VRangeToEnergyConverter():
52 theParticle(0), theLossTable(0), NumberOfElements(0), TotBin(300),
53 verboseLevel(1)
54{
55 fMaxEnergyCut = 0.;
56}
57
58G4VRangeToEnergyConverter::G4VRangeToEnergyConverter(const G4VRangeToEnergyConverter& right) : TotBin(right.TotBin)
59{
60 *this = right;
61}
62
63G4VRangeToEnergyConverter & G4VRangeToEnergyConverter::operator=(const G4VRangeToEnergyConverter &right)
64{
65 if (this == &right) return *this;
66 if (theLossTable) {
67 theLossTable->clearAndDestroy();
68 delete theLossTable;
69 theLossTable=0;
70 }
71
72 NumberOfElements = right.NumberOfElements;
73 //TotBin = right.TotBin;
74 theParticle = right.theParticle;
75 verboseLevel = right.verboseLevel;
76
77 // create the loss table
78 theLossTable = new G4LossTable();
79 theLossTable->reserve(G4Element::GetNumberOfElements());
80 // fill the loss table
81 for (size_t j=0; j<size_t(NumberOfElements); j++){
82 G4LossVector* aVector= new
83 G4LossVector(LowestEnergy, MaxEnergyCut, TotBin);
84 for (size_t i=0; i<size_t(TotBin); i++) {
85 G4double Value = (*((*right.theLossTable)[j]))[i];
86 aVector->PutValue(i,Value);
87 }
88 theLossTable->insert(aVector);
89 }
90
91 // clean up range vector store
92 for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
93 delete fRangeVectorStore.at(idx);
94 }
95 fRangeVectorStore.clear();
96
97 // copy range vector store
98 for (size_t j=0; j<((right.fRangeVectorStore).size()); j++){
99 G4RangeVector* vector = (right.fRangeVectorStore).at(j);
100 G4RangeVector* rangeVector = 0;
101 if (vector !=0 ) {
102 rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin);
103 for (size_t i=0; i<size_t(TotBin); i++) {
104 G4double Value = (*vector)[i];
105 rangeVector->PutValue(i,Value);
106 }
107 }
108 fRangeVectorStore.push_back(rangeVector);
109 }
110 return *this;
111}
112
113
114G4VRangeToEnergyConverter::~G4VRangeToEnergyConverter()
115{
116 Reset();
117}
118
119G4int G4VRangeToEnergyConverter::operator==(const G4VRangeToEnergyConverter &right) const
120{
121 return this == &right;
122}
123
124G4int G4VRangeToEnergyConverter::operator!=(const G4VRangeToEnergyConverter &right) const
125{
126 return this != &right;
127}
128
129
130// **********************************************************************
131// ************************* Convert ***********************************
132// **********************************************************************
133G4double G4VRangeToEnergyConverter::Convert(G4double rangeCut,
134 const G4Material* material)
135{
136#ifdef G4VERBOSE
137 if (GetVerboseLevel()>3) {
138 G4cout << "G4VRangeToEnergyConverter::Convert() ";
139 G4cout << "Convert for " << material->GetName()
140 << " with Range Cut " << rangeCut/mm << "[mm]" << G4endl;
141 }
142#endif
143
144 G4double theKineticEnergyCuts = 0.;
145
146 if (fMaxEnergyCut != MaxEnergyCut) {
147 fMaxEnergyCut = MaxEnergyCut;
148 // clear loss table and renge vectors
149 Reset();
150 }
151
152 // Build the energy loss table
153 BuildLossTable();
154
155 // Build range vector for every material, convert cut into energy-cut,
156 // fill theKineticEnergyCuts and delete the range vector
157 G4double tune = 0.025*mm*g/cm3 ,lowen = 30.*keV ;
158
159 // check density
160 G4double density = material->GetDensity() ;
161 if(density <= 0.) {
162 #ifdef G4VERBOSE
163 if (GetVerboseLevel()>0) {
164 G4cout << "G4VRangeToEnergyConverter::Convert() ";
165 G4cout << material->GetName() << "has zero density "
166 << "( " << density << ")" << G4endl;
167 }
168#endif
169 return 0.;
170 }
171
172 // initialize RangeVectorStore
173 const G4MaterialTable* table = G4Material::GetMaterialTable();
174 G4int ext_size = table->size() - fRangeVectorStore.size();
175 for (int i=0; i<ext_size; i++) fRangeVectorStore.push_back(0);
176
177 // Build Range Vector
178 G4int idx = material->GetIndex();
179 G4RangeVector* rangeVector = fRangeVectorStore.at(idx);
180 if (rangeVector == 0) {
181 rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin);
182 BuildRangeVector(material, rangeVector);
183 fRangeVectorStore.at(idx) = rangeVector;
184 }
185
186 // Convert Range Cut ro Kinetic Energy Cut
187 theKineticEnergyCuts = ConvertCutToKineticEnergy(rangeVector, rangeCut, idx);
188
189 if( ((theParticle->GetParticleName()=="e-")||(theParticle->GetParticleName()=="e+"))
190 && (theKineticEnergyCuts < lowen) ) {
191 // corr. should be switched on smoothly
192 theKineticEnergyCuts /= (1.+(1.-theKineticEnergyCuts/lowen)*
193 tune/(rangeCut*density));
194 }
195
196 if(theKineticEnergyCuts < LowestEnergy) {
197 theKineticEnergyCuts = LowestEnergy ;
198 } else if(theKineticEnergyCuts > MaxEnergyCut) {
199 theKineticEnergyCuts = MaxEnergyCut;
200 }
201
202 return theKineticEnergyCuts;
203}
204
205// **********************************************************************
206// ************************ SetEnergyRange *****************************
207// **********************************************************************
208void G4VRangeToEnergyConverter::SetEnergyRange(G4double lowedge,
209 G4double highedge)
210{
211 // check LowestEnergy/ HighestEnergy
212 if ( (lowedge<0.0)||(highedge<=lowedge) ){
213 G4cerr << "Error in G4VRangeToEnergyConverter::SetEnergyRange";
214 G4cerr << " : illegal energy range" << "(" << lowedge/GeV;
215 G4cerr << "," << highedge/GeV << ") [GeV]" << G4endl;
216 } else {
217 LowestEnergy = lowedge;
218 HighestEnergy = highedge;
219 }
220}
221
222
223G4double G4VRangeToEnergyConverter::GetLowEdgeEnergy()
224{
225 return LowestEnergy;
226}
227
228
229G4double G4VRangeToEnergyConverter::GetHighEdgeEnergy()
230{
231 return HighestEnergy;
232}
233
234// **********************************************************************
235// ******************* Get/SetMaxEnergyCut *****************************
236// **********************************************************************
237G4double G4VRangeToEnergyConverter::GetMaxEnergyCut()
238{
239 return MaxEnergyCut;
240}
241
242void G4VRangeToEnergyConverter::SetMaxEnergyCut(G4double value)
243{
244 MaxEnergyCut = value;
245}
246
247// **********************************************************************
248// ************************ Reset **************************************
249// **********************************************************************
250void G4VRangeToEnergyConverter::Reset()
251{
252 // delete loss table
253 if (theLossTable) {
254 theLossTable->clearAndDestroy();
255 delete theLossTable;
256 }
257 theLossTable=0;
258 NumberOfElements=0;
259
260 //clear RangeVectorStore
261 for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
262 delete fRangeVectorStore.at(idx);
263 }
264 fRangeVectorStore.clear();
265}
266
267
268// **********************************************************************
269// ************************ BuildLossTable ******************************
270// **********************************************************************
271// create Energy Loss Table for charged particles
272// (cross section tabel for neutral )
273void G4VRangeToEnergyConverter::BuildLossTable()
274{
275 if (size_t(NumberOfElements) == G4Element::GetNumberOfElements()) return;
276
277 // clear Loss table and Range vectors
278 Reset();
279
280 // Build dE/dx tables for elements
281 NumberOfElements = G4Element::GetNumberOfElements();
282 theLossTable = new G4LossTable();
283 theLossTable->reserve(G4Element::GetNumberOfElements());
284#ifdef G4VERBOSE
285 if (GetVerboseLevel()>3) {
286 G4cout << "G4VRangeToEnergyConverter::BuildLossTable() ";
287 G4cout << "Create theLossTable[" << theLossTable << "]";
288 G4cout << " NumberOfElements=" << NumberOfElements <<G4endl;
289 }
290#endif
291
292
293 // fill the loss table
294 for (size_t j=0; j<size_t(NumberOfElements); j++){
295 G4double Value;
296 G4LossVector* aVector= 0;
297 aVector= new G4LossVector(LowestEnergy, MaxEnergyCut, TotBin);
298 for (size_t i=0; i<size_t(TotBin); i++) {
299 Value = ComputeLoss( (*G4Element::GetElementTable())[j]->GetZ(),
300 aVector->GetLowEdgeEnergy(i)
301 );
302 aVector->PutValue(i,Value);
303 }
304 theLossTable->insert(aVector);
305 }
306}
307
308// **********************************************************************
309// ************************ BuildRangeVector ****************************
310// **********************************************************************
311void G4VRangeToEnergyConverter::BuildRangeVector(const G4Material* aMaterial,
312 G4PhysicsLogVector* rangeVector)
313{
314 // create range vector for a material
315 const G4ElementVector* elementVector = aMaterial->GetElementVector();
316 const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
317 G4int NumEl = aMaterial->GetNumberOfElements();
318
319 // calculate parameters of the low energy part first
320 size_t i;
321 std::vector<G4double> lossV;
322 for ( size_t ib=0; ib<size_t(TotBin); ib++) {
323 G4double loss=0.;
324 for (i=0; i<size_t(NumEl); i++) {
325 G4int IndEl = (*elementVector)[i]->GetIndex();
326 loss += atomicNumDensityVector[i]*
327 (*((*theLossTable)[IndEl]))[ib];
328 }
329 lossV.push_back(loss);
330 }
331
332 // Integrate with Simpson formula with logarithmic binning
333 G4double ltt = std::log(MaxEnergyCut/LowestEnergy);
334 G4double dltau = ltt/TotBin;
335
336 G4double s0 = 0.;
337 G4double Value;
338 for ( i=0; i<size_t(TotBin); i++) {
339 G4double t = rangeVector->GetLowEdgeEnergy(i);
340 G4double s = t/lossV[i];
341 if (i==0) s0 += 0.5*s;
342 else s0 += s;
343
344 if (i==0) {
345 Value = (s0 + 0.5*s)*dltau ;
346 } else {
347 Value = (s0 - 0.5*s)*dltau ;
348 }
349 rangeVector->PutValue(i,Value);
350 }
351}
352
353// **********************************************************************
354// ****************** ConvertCutToKineticEnergy *************************
355// **********************************************************************
356G4double G4VRangeToEnergyConverter::ConvertCutToKineticEnergy(
357 G4RangeVector* rangeVector,
358 G4double theCutInLength,
359 size_t materialIndex
360 ) const
361{
362 const G4double epsilon=0.01;
363
364 // find max. range and the corresponding energy (rmax,Tmax)
365 G4double rmax= -1.e10*mm;
366
367 G4double T1 = LowestEnergy;
368 G4double r1 =(*rangeVector)[0] ;
369
370 G4double T2 = MaxEnergyCut;
371
372 // check theCutInLength < r1
373 if ( theCutInLength <= r1 ) { return T1; }
374
375 // scan range vector to find nearest bin
376 // ( suppose that r(Ti) > r(Tj) if Ti >Tj )
377 for (size_t ibin=0; ibin<size_t(TotBin); ibin++) {
378 G4double T=rangeVector->GetLowEdgeEnergy(ibin);
379 G4double r=(*rangeVector)[ibin];
380 if ( r>rmax ) rmax=r;
381 if (r <theCutInLength ) {
382 T1 = T;
383 r1 = r;
384 } else if (r >theCutInLength ) {
385 T2 = T;
386 break;
387 }
388 }
389
390 // check cut in length is smaller than range max
391 if ( theCutInLength >= rmax ) {
392#ifdef G4VERBOSE
393 if (GetVerboseLevel()>2) {
394 G4cout << "G4VRangeToEnergyConverter::ConvertCutToKineticEnergy ";
395 G4cout << " for " << theParticle->GetParticleName() << G4endl;
396 G4cout << "The cut in range [" << theCutInLength/mm << " (mm)] ";
397 G4cout << " is too big " ;
398 G4cout << " for material idx=" << materialIndex <<G4endl;
399 }
400#endif
401 return MaxEnergyCut;
402 }
403
404 // convert range to energy
405 G4double T3 = std::sqrt(T1*T2);
406 G4double r3 = rangeVector->Value(T3);
407 while ( std::fabs(1.-r3/theCutInLength)>epsilon ) {
408 if ( theCutInLength <= r3 ) {
409 T2 = T3;
410 } else {
411 T1 = T3;
412 }
413 T3 = std::sqrt(T1*T2);
414 r3 = rangeVector->Value(T3);
415 }
416
417 return T3;
418}
419
Note: See TracBrowser for help on using the repository browser.