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