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