| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4VRangeToEnergyConverter class implementat <<
27 // 23 //
28 // Author: H.Kurashige, 05 October 2002 - Firs << 24 // $Id: G4VRangeToEnergyConverter.cc,v 1.3 2004/12/02 06:53:56 kurasige Exp $
29 // ------------------------------------------- << 25 // GEANT4 tag $Name: geant4-07-01 $
>> 26 //
>> 27 //
>> 28 // --------------------------------------------------------------
>> 29 // GEANT 4 class implementation file/ History:
>> 30 // 5 Oct. 2002, H.Kuirashige : Structure created based on object model
>> 31 // --------------------------------------------------------------
30 32
31 #include "G4VRangeToEnergyConverter.hh" 33 #include "G4VRangeToEnergyConverter.hh"
32 #include "G4ParticleTable.hh" 34 #include "G4ParticleTable.hh"
33 #include "G4Element.hh" << 35 #include "G4Material.hh"
34 #include "G4SystemOfUnits.hh" << 36 #include "G4PhysicsLogVector.hh"
35 #include "G4Log.hh" <<
36 #include "G4Exp.hh" <<
37 #include "G4AutoLock.hh" <<
38 37
39 namespace << 38 #include "G4ios.hh"
40 { << 39 #include <iomanip>
41 G4Mutex theREMutex = G4MUTEX_INITIALIZER; << 40 #include <strstream>
42 } <<
43 41
44 G4double G4VRangeToEnergyConverter::sEmin = CL << 42 // energy range
45 G4double G4VRangeToEnergyConverter::sEmax = 10 << 43 G4double G4VRangeToEnergyConverter::LowestEnergy = 0.99e-3*MeV;
>> 44 G4double G4VRangeToEnergyConverter::HighestEnergy = 100.0e6*MeV;
46 45
47 std::vector<G4double>* G4VRangeToEnergyConvert << 46 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter():
>> 47 theParticle(0), theLossTable(0), NumberOfElements(0), TotBin(200),
>> 48 verboseLevel(1)
>> 49 {
>> 50 }
48 51
49 G4int G4VRangeToEnergyConverter::sNbinPerDecad << 52 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter(const G4VRangeToEnergyConverter& right)
50 G4int G4VRangeToEnergyConverter::sNbin = 350; << 53 {
>> 54 *this = right;
>> 55 }
51 56
52 // ------------------------------------------- << 57 G4VRangeToEnergyConverter & G4VRangeToEnergyConverter::operator=(const G4VRangeToEnergyConverter &right)
53 G4VRangeToEnergyConverter::G4VRangeToEnergyCon <<
54 { 58 {
55 if(nullptr == sEnergy) << 59 if (this == &right) return *this;
56 { << 60 if (theLossTable) delete theLossTable;
57 G4AutoLock l(&theREMutex); << 61
58 if(nullptr == sEnergy) << 62 NumberOfElements = right.NumberOfElements;
59 { << 63 TotBin = right.TotBin;
60 isFirstInstance = true; << 64 theParticle = right.theParticle;
>> 65 verboseLevel = right.verboseLevel;
>> 66
>> 67 // create the loss table
>> 68 theLossTable = new G4LossTable();
>> 69 theLossTable->reserve(G4Element::GetNumberOfElements());
>> 70 // fill the loss table
>> 71 for (size_t j=0; j<size_t(NumberOfElements); j++){
>> 72 G4LossVector* aVector= new
>> 73 G4LossVector(LowestEnergy, HighestEnergy, TotBin);
>> 74 for (size_t i=0; i<size_t(TotBin); i++) {
>> 75 G4double Value = (*((*right.theLossTable)[j]))[i];
>> 76 aVector->PutValue(i,Value);
61 } 77 }
62 l.unlock(); << 78 theLossTable->insert(aVector);
63 } <<
64 // this method defines lock itself <<
65 if(isFirstInstance) <<
66 { <<
67 FillEnergyVector(CLHEP::keV, 10.0*CLHEP::G <<
68 } 79 }
>> 80 return *this;
69 } 81 }
70 82
71 // ------------------------------------------- << 83
72 G4VRangeToEnergyConverter::~G4VRangeToEnergyCo 84 G4VRangeToEnergyConverter::~G4VRangeToEnergyConverter()
>> 85 {
>> 86 if (theLossTable) delete theLossTable;
>> 87 }
>> 88
>> 89 G4int G4VRangeToEnergyConverter::operator==(const G4VRangeToEnergyConverter &right) const
73 { 90 {
74 if(isFirstInstance) << 91 return this == &right;
75 { <<
76 delete sEnergy; <<
77 sEnergy = nullptr; <<
78 sEmin = CLHEP::keV; <<
79 sEmax = 10.*CLHEP::GeV; <<
80 } <<
81 } 92 }
82 93
83 // ------------------------------------------- << 94 G4int G4VRangeToEnergyConverter::operator!=(const G4VRangeToEnergyConverter &right) const
84 G4double G4VRangeToEnergyConverter::Convert(co <<
85 co <<
86 { 95 {
87 #ifdef G4VERBOSE << 96 return this != &right;
88 if (GetVerboseLevel()>3) << 97 }
89 { <<
90 G4cout << "G4VRangeToEnergyConverter::Conv <<
91 G4cout << "Convert for " << material->GetN <<
92 << " with Range Cut " << rangeCut/mm << " <<
93 } <<
94 #endif <<
95 98
96 G4double cut = 0.0; <<
97 if(fPDG == 22) <<
98 { <<
99 cut = ConvertForGamma(rangeCut, material); <<
100 } <<
101 else <<
102 { <<
103 cut = ConvertForElectron(rangeCut, materia <<
104 99
105 const G4double tune = 0.025*CLHEP::mm*CLHE << 100 // **********************************************************************
106 const G4double lowen = 30.*CLHEP::keV; << 101 // ************************* Convert ***********************************
107 if(cut < lowen) << 102 // **********************************************************************
108 { << 103 G4double G4VRangeToEnergyConverter::Convert(G4double rangeCut,
109 // corr. should be switched on smoothly << 104 const G4Material* material)
110 cut /= (1.+(1.-cut/lowen)*tune/(rangeCut << 105 {
>> 106 //???????????? G4double Charge = theParticle->GetPDGCharge();
>> 107 G4double Mass = theParticle->GetPDGMass();
>> 108 G4double theKineticEnergyCuts = 0.;
>> 109
>> 110 // Build the energy loss table
>> 111 if (theLossTable ==0) BuildLossTable();
>> 112
>> 113 // Build range vector for every material, convert cut into energy-cut,
>> 114 // fill theKineticEnergyCuts and delete the range vector
>> 115 G4double tune = 0.025*mm*g/cm3 ,lowen = 30.*keV ;
>> 116
>> 117 G4int idx = material->GetIndex();
>> 118 G4double density = material->GetDensity() ;
>> 119 if(density > 0.) {
>> 120 G4RangeVector* rangeVector = new G4RangeVector(LowestEnergy, HighestEnergy, TotBin);
>> 121 BuildRangeVector(material, HighestEnergy, Mass, rangeVector);
>> 122 theKineticEnergyCuts = ConvertCutToKineticEnergy(rangeVector, rangeCut, idx);
>> 123
>> 124 if( ((theParticle->GetParticleName()=="e-")||(theParticle->GetParticleName()=="e+"))
>> 125 && (theKineticEnergyCuts < lowen) )
>> 126 { theKineticEnergyCuts /= (1.+tune/(rangeCut*density)); }
>> 127 if(theKineticEnergyCuts < LowestEnergy) {
>> 128 theKineticEnergyCuts = LowestEnergy ;
111 } 129 }
>> 130 delete rangeVector;
112 } 131 }
113 << 132 return theKineticEnergyCuts;
114 cut = std::max(sEmin, std::min(cut, sEmax)); <<
115 return cut; <<
116 } 133 }
117 134
118 // ------------------------------------------- << 135 // **********************************************************************
119 void G4VRangeToEnergyConverter::SetEnergyRange << 136 // ************************ SetEnergyRange *****************************
120 << 137 // **********************************************************************
121 { << 138 void G4VRangeToEnergyConverter::SetEnergyRange(G4double lowedge,
122 G4double ehigh = std::min(10.*CLHEP::GeV, hi << 139 G4double highedge)
123 if(ehigh > lowedge) << 140 {
124 { << 141 // check LowestEnergy/ HighestEnergy
125 FillEnergyVector(lowedge, ehigh); << 142 if ( (lowedge<0.0)||(highedge<=lowedge) ){
126 } << 143 G4cerr << "Error in G4VRangeToEnergyConverter::SetEnergyRange";
>> 144 G4cerr << " : illegal energy range" << "(" << lowedge/GeV;
>> 145 G4cerr << "," << highedge/GeV << ") [GeV]" << G4endl;
>> 146 } else {
>> 147 LowestEnergy = lowedge;
>> 148 HighestEnergy = highedge;
>> 149 }
127 } 150 }
128 151
129 // ------------------------------------------- << 152
130 G4double G4VRangeToEnergyConverter::GetLowEdge 153 G4double G4VRangeToEnergyConverter::GetLowEdgeEnergy()
131 { 154 {
132 return sEmin; << 155 return LowestEnergy;
133 } 156 }
134 157
135 // ------------------------------------------- << 158
136 G4double G4VRangeToEnergyConverter::GetHighEdg 159 G4double G4VRangeToEnergyConverter::GetHighEdgeEnergy()
137 { 160 {
138 return sEmax; << 161 return HighestEnergy;
139 } 162 }
140 163
141 // ------------------------------------------- << 164 // **********************************************************************
>> 165 // ************************ RangeLinSimpson *****************************
>> 166 // **********************************************************************
>> 167 G4double G4VRangeToEnergyConverter::RangeLinSimpson(
>> 168 G4int numberOfElement,
>> 169 const G4ElementVector* elementVector,
>> 170 const G4double* atomicNumDensityVector,
>> 171 G4double aMass,
>> 172 G4double taulow, G4double tauhigh, G4int nbin)
>> 173 {
>> 174 // Simpson numerical integration, linear binning
>> 175 G4double dtau = (tauhigh-taulow)/nbin;
>> 176 G4double Value=0.;
>> 177 for (size_t i=0; i<=size_t(nbin); i++){
>> 178 G4double taui=taulow+dtau*i;
>> 179 G4double ti=aMass*taui;
>> 180 G4double lossi=0.;
>> 181 size_t nEl = (size_t)(numberOfElement);
>> 182 for (size_t j=0; j<nEl; j++) {
>> 183 G4bool isOut;
>> 184 G4int IndEl = (*elementVector)[j]->GetIndex();
>> 185 lossi += atomicNumDensityVector[j]*
>> 186 (*theLossTable)[IndEl]->GetValue(ti,isOut);
>> 187 }
>> 188 if ( i==0 ) {
>> 189 Value += 0.5/lossi;
>> 190 } else {
>> 191 if ( i<size_t(nbin) ) Value += 1./lossi;
>> 192 else Value += 0.5/lossi;
>> 193 }
>> 194 }
>> 195 Value *= aMass*dtau;
142 196
143 G4double G4VRangeToEnergyConverter::GetMaxEner << 197 return Value;
144 { <<
145 return sEmax; <<
146 } 198 }
147 199
148 // ------------------------------------------- << 200
149 void G4VRangeToEnergyConverter::SetMaxEnergyCu << 201 // **********************************************************************
150 { << 202 // ************************ RangeLogSimpson *****************************
151 G4double ehigh = std::min(10.*CLHEP::GeV, va << 203 // **********************************************************************
152 if(ehigh > sEmin) << 204 G4double G4VRangeToEnergyConverter::RangeLogSimpson(
153 { << 205 G4int numberOfElement,
154 FillEnergyVector(sEmin, ehigh); << 206 const G4ElementVector* elementVector,
>> 207 const G4double* atomicNumDensityVector,
>> 208 G4double aMass,
>> 209 G4double ltaulow, G4double ltauhigh,
>> 210 G4int nbin)
>> 211 {
>> 212 // Simpson numerical integration, logarithmic binning
>> 213 if(nbin<0) nbin = TotBin;
>> 214 G4double ltt = ltauhigh-ltaulow;
>> 215 G4double dltau = ltt/nbin;
>> 216 G4double Value = 0.;
>> 217 for (size_t i=0; i<=size_t(nbin); i++){
>> 218 G4double ui = ltaulow+dltau*i;
>> 219 G4double taui = std::exp(ui);
>> 220 G4double ti = aMass*taui;
>> 221 G4double lossi = 0.;
>> 222 size_t nEl = (size_t)(numberOfElement);
>> 223
>> 224 for (size_t j=0; j<nEl; j++) {
>> 225 G4bool isOut;
>> 226 G4int IndEl = (*elementVector)[j]->GetIndex();
>> 227 lossi += atomicNumDensityVector[j]*
>> 228 (*theLossTable)[IndEl]->GetValue(ti,isOut);
>> 229 }
>> 230 if ( i==0 ) {
>> 231 Value += 0.5*taui/lossi;
>> 232 } else {
>> 233 if ( i<size_t(nbin) ) Value += taui/lossi;
>> 234 else Value += 0.5*taui/lossi;
>> 235 }
155 } 236 }
>> 237 Value *= aMass*dltau;
>> 238
>> 239 return Value;
156 } 240 }
157 241
158 // ------------------------------------------- << 242 // **********************************************************************
159 void G4VRangeToEnergyConverter::FillEnergyVect << 243 // ************************ BuildLossTable ******************************
160 << 244 // **********************************************************************
161 { << 245 // create Energy Loss Table for charged particles
162 if(emin != sEmin || emax != sEmax || nullptr << 246 // (cross section tabel for neutral )
163 { << 247 void G4VRangeToEnergyConverter::BuildLossTable()
164 sEmin = emin; << 248 {
165 sEmax = emax; << 249 // Build dE/dx tables for elements
166 sNbin = sNbinPerDecade*G4lrint(std::log10( << 250 if (size_t(NumberOfElements) != G4Element::GetNumberOfElements()) {
167 if(nullptr == sEnergy) { sEnergy = new std << 251 if (theLossTable!=0) delete theLossTable;
168 sEnergy->resize(sNbin + 1); << 252 theLossTable =0;
169 (*sEnergy)[0] = emin; << 253 NumberOfElements = 0;
170 (*sEnergy)[sNbin] = emax; << 254 }
171 G4double fact = G4Log(emax/emin)/sNbin; << 255
172 for(G4int i=1; i<sNbin; ++i) { (*sEnergy)[ << 256 if (NumberOfElements ==0) {
>> 257 NumberOfElements = G4Element::GetNumberOfElements();
>> 258 theLossTable = new G4LossTable();
>> 259 theLossTable->reserve(G4Element::GetNumberOfElements());
>> 260 #ifdef G4VERBOSE
>> 261 if (GetVerboseLevel()>2) {
>> 262 G4cout << "G4VRangeToEnergyConverter::BuildLossTable() ";
>> 263 G4cout << "Create theLossTable[" << theLossTable << "]";
>> 264 G4cout << " NumberOfElements=" << NumberOfElements <<G4endl;
>> 265 }
>> 266 #endif
>> 267 }
>> 268
>> 269 // fill the loss table
>> 270 for (size_t j=0; j<size_t(NumberOfElements); j++){
>> 271 G4double Value;
>> 272 G4LossVector* aVector= new
>> 273 G4LossVector(LowestEnergy, HighestEnergy, TotBin);
>> 274 for (size_t i=0; i<size_t(TotBin); i++) {
>> 275 Value = ComputeLoss( (*G4Element::GetElementTable())[j]->GetZ(),
>> 276 aVector->GetLowEdgeEnergy(i)
>> 277 );
>> 278 aVector->PutValue(i,Value);
>> 279 }
>> 280 theLossTable->insert(aVector);
173 } 281 }
174 } 282 }
175 283
176 // ------------------------------------------- << 284 // **********************************************************************
177 G4double << 285 // ************************** ComputeLoss *******************************
178 G4VRangeToEnergyConverter::ConvertForGamma(con << 286 // **********************************************************************
179 con << 287 G4double G4VRangeToEnergyConverter::ComputeLoss(G4double AtomicNumber,
180 { << 288 G4double KineticEnergy) const
181 const G4ElementVector* elm = material->GetEl << 289 {
182 const G4double* dens = material->GetAtomicNu << 290 // calculate dE/dx
183 << 291
184 // fill absorption length vector << 292 static G4double Z;
185 G4int nelm = (G4int)material->GetNumberOfEle << 293 static G4double ionpot, tau0, taum, taul, ca, cba, cc;
186 G4double range1 = 0.0; << 294
187 G4double range2 = 0.0; << 295 G4double z2Particle = theParticle->GetPDGCharge()/eplus;
188 G4double e1 = 0.0; << 296 z2Particle *= z2Particle;
189 G4double e2 = 0.0; << 297 if (z2Particle < 0.1) return 0.0;
190 for (G4int i=0; i<sNbin; ++i) << 298
191 { << 299 if( std::abs(AtomicNumber-Z)>0.1 ){
192 e2 = (*sEnergy)[i]; << 300 // recalculate constants
193 G4double sig = 0.; << 301 Z = AtomicNumber;
194 << 302 G4double Z13 = std::exp(std::log(Z)/3.);
195 for (G4int j=0; j<nelm; ++j) << 303 tau0 = 0.1*Z13*MeV/proton_mass_c2;
196 { << 304 taum = 0.035*Z13*MeV/proton_mass_c2;
197 sig += dens[j]*ComputeValue((*elm)[j]->G << 305 taul = 2.*MeV/proton_mass_c2;
198 } << 306 ionpot = 1.6e-5*MeV*std::exp(0.9*std::log(Z));
199 range2 = (sig > 0.0) ? 5./sig : DBL_MAX; << 307 cc = (taul+1.)*(taul+1.)*std::log(2.*electron_mass_c2*taul*(taul+2.)/ionpot)/(taul*(taul+2.))-1.;
200 if(i == 0 || range2 < rangeCut) << 308 cc = 2.*twopi_mc2_rcl2*Z*cc*std::sqrt(taul);
201 { << 309 ca = cc/((1.-0.5*std::sqrt(tau0/taum))*tau0);
202 e1 = e2; << 310 cba = -0.5/std::sqrt(taum);
203 range1 = range2; << 311 }
204 } << 312
205 else << 313 G4double tau = KineticEnergy/theParticle->GetPDGMass();
206 { << 314 G4double dEdx;
207 break; << 315 if ( tau <= tau0 ) {
208 } << 316 dEdx = ca*(std::sqrt(tau)+cba*tau);
209 } << 317 } else {
210 return LiniearInterpolation(e1, e2, range1, << 318 if( tau <= taul ) {
211 } << 319 dEdx = cc/std::sqrt(tau);
212 << 320 } else {
213 // ------------------------------------------- << 321 dEdx = (tau+1.)*(tau+1.)*
214 G4double << 322 std::log(2.*electron_mass_c2*tau*(tau+2.)/ionpot)/(tau*(tau+2.))-1.;
215 G4VRangeToEnergyConverter::ConvertForElectron( << 323 dEdx = 2.*twopi_mc2_rcl2*Z*dEdx;
216 <<
217 { <<
218 const G4ElementVector* elm = material->GetEl <<
219 const G4double* dens = material->GetAtomicNu <<
220 <<
221 // fill absorption length vector <<
222 G4int nelm = (G4int)material->GetNumberOfEle <<
223 G4double dedx1 = 0.0; <<
224 G4double dedx2 = 0.0; <<
225 G4double range1 = 0.0; <<
226 G4double range2 = 0.0; <<
227 G4double e1 = 0.0; <<
228 G4double e2 = 0.0; <<
229 G4double range = 0.; <<
230 for (G4int i=0; i<sNbin; ++i) <<
231 { <<
232 e2 = (*sEnergy)[i]; <<
233 dedx2 = 0.0; <<
234 for (G4int j=0; j<nelm; ++j) <<
235 { <<
236 dedx2 += dens[j]*ComputeValue((*elm)[j]- <<
237 } 324 }
238 range += (dedx1 + dedx2 > 0.0) ? 2*(e2 - e << 325 }
239 range2 = range; << 326 return dEdx*z2Particle ;
240 if(range2 < rangeCut) << 327 }
241 { << 328
242 e1 = e2; << 329 // **********************************************************************
243 dedx1 = dedx2; << 330 // ************************ BuildRangeVector ****************************
244 range1 = range2; << 331 // **********************************************************************
>> 332 void G4VRangeToEnergyConverter::BuildRangeVector(
>> 333 const G4Material* aMaterial,
>> 334 G4double maxEnergy,
>> 335 G4double aMass,
>> 336 G4RangeVector* rangeVector)
>> 337 {
>> 338 // create range vector for a material
>> 339 const G4double tlim=2.*MeV, t1=0.1*MeV, t2=0.025*MeV;
>> 340 const G4int maxnbint=100;
>> 341
>> 342 const G4ElementVector* elementVector = aMaterial->GetElementVector();
>> 343 const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
>> 344
>> 345 G4int NumEl = aMaterial->GetNumberOfElements();
>> 346
>> 347 // calculate parameters of the low energy part first
>> 348 G4double loss1=0.;
>> 349 G4double loss2=0.;
>> 350 size_t i;
>> 351 for (i=0; i<size_t(NumEl); i++) {
>> 352 G4bool isOut;
>> 353 G4int IndEl = (*elementVector)[i]->GetIndex();
>> 354 loss1 += atomicNumDensityVector[i]*
>> 355 (*theLossTable)[IndEl]->GetValue(t1,isOut);
>> 356 loss2 += atomicNumDensityVector[i]*
>> 357 (*theLossTable)[IndEl]->GetValue(t2,isOut);
>> 358 }
>> 359 G4double tau1 = t1/proton_mass_c2;
>> 360 G4double sqtau1 = std::sqrt(tau1);
>> 361 G4double ca = (4.*loss2-loss1)/sqtau1;
>> 362 G4double cb = (2.*loss1-4.*loss2)/tau1;
>> 363 G4double cba = cb/ca;
>> 364 G4double taulim = tlim/proton_mass_c2;
>> 365 G4double taumax = maxEnergy/aMass;
>> 366 G4double ltaumax = std::log(taumax);
>> 367
>> 368 // now we can fill the range vector....
>> 369 G4double rmax = 0.0;
>> 370 for (i=0; i<size_t(TotBin); i++) {
>> 371 G4double LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i);
>> 372 G4double tau = LowEdgeEnergy/aMass;
>> 373 G4double Value;
>> 374
>> 375 if ( tau <= tau1 ){
>> 376 Value =2.*aMass*std::log(1.+cba*std::sqrt(tau))/cb;
>> 377 } else {
>> 378 Value = 2.*aMass*std::log(1.+cba*sqtau1)/cb;
>> 379 if ( tau <= taulim ) {
>> 380 G4int nbin = (G4int)(maxnbint*(tau-tau1)/(taulim-tau1));
>> 381 if ( nbin<1 ) nbin = 1;
>> 382 Value += RangeLinSimpson( NumEl, elementVector,
>> 383 atomicNumDensityVector, aMass,
>> 384 tau1, tau, nbin);
>> 385 } else {
>> 386 Value += RangeLinSimpson( NumEl, elementVector,
>> 387 atomicNumDensityVector, aMass,
>> 388 tau1, taulim, maxnbint);
>> 389 G4double ltaulow = std::log(taulim);
>> 390 G4double ltauhigh = std::log(tau);
>> 391 G4int nbin = (G4int)(maxnbint*(ltauhigh-ltaulow)/(ltaumax-ltaulow));
>> 392 if ( nbin<1 ) nbin = 1;
>> 393 Value += RangeLogSimpson(NumEl, elementVector,
>> 394 atomicNumDensityVector, aMass,
>> 395 ltaulow, ltauhigh, nbin);
>> 396 }
>> 397 }
>> 398 rangeVector->PutValue(i,Value);
>> 399 if (rmax < Value) rmax = Value;
>> 400 }
>> 401 }
>> 402
>> 403 // **********************************************************************
>> 404 // ****************** ConvertCutToKineticEnergy *************************
>> 405 // **********************************************************************
>> 406 G4double G4VRangeToEnergyConverter::ConvertCutToKineticEnergy(
>> 407 G4RangeVector* rangeVector,
>> 408 G4double theCutInLength,
>> 409 size_t materialIndex
>> 410 ) const
>> 411 {
>> 412 const G4double epsilon=0.01;
>> 413
>> 414 // find max. range and the corresponding energy (rmax,Tmax)
>> 415 G4double rmax= -1.e10*mm;
>> 416 G4double Tmax= HighestEnergy;
>> 417 G4double fac = std::exp( std::log(HighestEnergy/LowestEnergy)/TotBin );
>> 418 G4double T=LowestEnergy/fac;
>> 419 G4bool isOut;
>> 420
>> 421 for (size_t ibin=0; ibin<size_t(TotBin); ibin++) {
>> 422 T *= fac;
>> 423 G4double r=rangeVector->GetValue(T,isOut);
>> 424 if ( r>rmax ) {
>> 425 Tmax=T;
>> 426 rmax=r;
>> 427 }
>> 428 }
>> 429
>> 430 // check cut in length is smaller than range max
>> 431 if ( theCutInLength >= rmax ) {
>> 432 #ifdef G4VERBOSE
>> 433 if (GetVerboseLevel()>0) {
>> 434 G4cout << "G4VRangeToEnergyConverter::ConvertCutToKineticEnergy ";
>> 435 G4cout << " for " << theParticle->GetParticleName() << G4endl;
>> 436 G4cout << "The cut in range [" << theCutInLength/mm << " (mm)] ";
>> 437 G4cout << " is too big " ;
>> 438 G4cout << " for material idx=" << materialIndex <<G4endl;
>> 439 G4cout << "The cut in energy is set" << DBL_MAX/GeV << "GeV " <<G4endl;
245 } 440 }
246 else << 441 #endif
247 { << 442 return DBL_MAX;
248 break; << 443 }
>> 444
>> 445 // convert range to energy
>> 446 G4double T1 = LowestEnergy;
>> 447 G4double r1 = rangeVector->GetValue(T1,isOut);
>> 448 if ( theCutInLength <= r1 )
>> 449 {
>> 450 return T1;
>> 451 }
>> 452
>> 453 G4double T2 = Tmax ;
>> 454 G4double T3 = std::sqrt(T1*T2);
>> 455 G4double r3 = rangeVector->GetValue(T3,isOut);
>> 456 while ( std::abs(1.-r3/theCutInLength)>epsilon ) {
>> 457 if ( theCutInLength <= r3 ) {
>> 458 T2 = T3;
>> 459 } else {
>> 460 T1 = T3;
249 } 461 }
>> 462 T3 = std::sqrt(T1*T2);
>> 463 r3 = rangeVector->GetValue(T3,isOut);
250 } 464 }
251 return LiniearInterpolation(e1, e2, range1, << 465
>> 466 return T3;
252 } 467 }
253 468
254 // ------------------------------------------- <<
255 469