| 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 93090 2015-10-05 13:14:43Z 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()
>> 118 {
>> 119 Reset();
>> 120 // Comment out Reset() for MT application
>> 121
>> 122 /////MA // delete loss table without deleteing vectors
>> 123 /////MA if (theLossTable) {
>> 124 /////MA delete theLossTable;
>> 125 /////MA }
>> 126 /////MA theLossTable=0;
>> 127 /////MA NumberOfElements=0;
>> 128 /////MA
>> 129 /////MA //clear RangeVectorStore without deleteing vectors
>> 130 /////MA fRangeVectorStore.clear();
>> 131
>> 132 }
>> 133
>> 134 G4int G4VRangeToEnergyConverter::operator==(const G4VRangeToEnergyConverter &right) const
73 { 135 {
74 if(isFirstInstance) << 136 return this == &right;
75 { << 137 }
76 delete sEnergy; << 138
77 sEnergy = nullptr; << 139 G4int G4VRangeToEnergyConverter::operator!=(const G4VRangeToEnergyConverter &right) const
78 sEmin = CLHEP::keV; << 140 {
79 sEmax = 10.*CLHEP::GeV; << 141 return this != &right;
80 } <<
81 } 142 }
82 143
83 // ------------------------------------------- << 144
84 G4double G4VRangeToEnergyConverter::Convert(co << 145 // **********************************************************************
85 co << 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 const size_t MAX_LOOP = 1000;
>> 435 for (size_t loop_count=0; loop_count<MAX_LOOP; ++loop_count){
>> 436 if (std::fabs(1.-r3/theCutInLength)<epsilon ) break;
>> 437 if ( theCutInLength <= r3 ) {
>> 438 T2 = T3;
>> 439 } else {
>> 440 T1 = T3;
249 } 441 }
>> 442 T3 = std::sqrt(T1*T2);
>> 443 r3 = rangeVector->Value(T3);
250 } 444 }
251 return LiniearInterpolation(e1, e2, range1, << 445
>> 446 return T3;
252 } 447 }
253 448
254 // ------------------------------------------- <<
255 449