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