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 // 26 // >> 27 // $Id: G4RDVeLowEnergyLoss.cc 107396 2017-11-10 08:28:08Z gcosmo $ >> 28 // GEANT4 tag $Name: geant4-09-01-ref-00 $ 27 // 29 // 28 // 30 // 29 // ------------------------------------------- 31 // -------------------------------------------------------------- 30 // GEANT 4 class implementation file 32 // GEANT 4 class implementation file 31 // 33 // 32 // History: first implementation, based on ob 34 // History: first implementation, based on object model of 33 // 2nd December 1995, G.Cosmo 35 // 2nd December 1995, G.Cosmo 34 // ------------------------------------------- 36 // -------------------------------------------------------------- 35 // 37 // 36 // Modifications: 38 // Modifications: 37 // 20/09/00 update fluctuations V.Ivanchenko 39 // 20/09/00 update fluctuations V.Ivanchenko 38 // 22/11/00 minor fix in fluctuations V.Ivanch 40 // 22/11/00 minor fix in fluctuations V.Ivanchenko 39 // 10/05/01 V.Ivanchenko Clean up againist Li 41 // 10/05/01 V.Ivanchenko Clean up againist Linux compilation with -Wall 40 // 22/05/01 V.Ivanchenko Update range calcula 42 // 22/05/01 V.Ivanchenko Update range calculation 41 // 23/11/01 V.Ivanchenko Move static member-f 43 // 23/11/01 V.Ivanchenko Move static member-functions from header to source 42 // 22/01/03 V.Ivanchenko Cut per region 44 // 22/01/03 V.Ivanchenko Cut per region 43 // 11/02/03 V.Ivanchenko Add limits to fluctu 45 // 11/02/03 V.Ivanchenko Add limits to fluctuations 44 // 24/04/03 V.Ivanchenko Fix the problem of t 46 // 24/04/03 V.Ivanchenko Fix the problem of table size 45 // 47 // 46 // ------------------------------------------- 48 // -------------------------------------------------------------- 47 49 48 #include "G4RDVeLowEnergyLoss.hh" 50 #include "G4RDVeLowEnergyLoss.hh" 49 #include "G4PhysicalConstants.hh" 51 #include "G4PhysicalConstants.hh" 50 #include "G4SystemOfUnits.hh" 52 #include "G4SystemOfUnits.hh" 51 #include "G4ProductionCutsTable.hh" 53 #include "G4ProductionCutsTable.hh" 52 54 53 G4double G4RDVeLowEnergyLoss::ParticleMass 55 G4double G4RDVeLowEnergyLoss::ParticleMass ; 54 G4double G4RDVeLowEnergyLoss::taulow 56 G4double G4RDVeLowEnergyLoss::taulow ; 55 G4double G4RDVeLowEnergyLoss::tauhigh 57 G4double G4RDVeLowEnergyLoss::tauhigh ; 56 G4double G4RDVeLowEnergyLoss::ltaulow 58 G4double G4RDVeLowEnergyLoss::ltaulow ; 57 G4double G4RDVeLowEnergyLoss::ltauhigh 59 G4double G4RDVeLowEnergyLoss::ltauhigh ; 58 60 59 61 60 G4bool G4RDVeLowEnergyLoss::rndmStepFlag 62 G4bool G4RDVeLowEnergyLoss::rndmStepFlag = false; 61 G4bool G4RDVeLowEnergyLoss::EnlossFlucFl 63 G4bool G4RDVeLowEnergyLoss::EnlossFlucFlag = true; 62 G4double G4RDVeLowEnergyLoss::dRoverRange 64 G4double G4RDVeLowEnergyLoss::dRoverRange = 20*perCent; 63 G4double G4RDVeLowEnergyLoss::finalRange 65 G4double G4RDVeLowEnergyLoss::finalRange = 200*micrometer; 64 G4double G4RDVeLowEnergyLoss::c1lim = dRov 66 G4double G4RDVeLowEnergyLoss::c1lim = dRoverRange ; 65 G4double G4RDVeLowEnergyLoss::c2lim = 2.*( 67 G4double G4RDVeLowEnergyLoss::c2lim = 2.*(1.-dRoverRange)*finalRange ; 66 G4double G4RDVeLowEnergyLoss::c3lim = -(1. 68 G4double G4RDVeLowEnergyLoss::c3lim = -(1.-dRoverRange)*finalRange*finalRange; 67 69 68 70 69 // 71 // 70 72 71 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss() 73 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss() 72 :G4VContinuousDiscreteProce 74 :G4VContinuousDiscreteProcess("No Name Loss Process"), 73 lastMaterial(0), 75 lastMaterial(0), 74 nmaxCont1(4), 76 nmaxCont1(4), 75 nmaxCont2(16) 77 nmaxCont2(16) 76 { 78 { 77 G4Exception("G4RDVeLowEnergyLoss::G4RDVeLowE 79 G4Exception("G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()", "InvalidCall", 78 FatalException, "Default constru 80 FatalException, "Default constructor is called!"); 79 } 81 } 80 82 81 // 83 // 82 84 83 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(const 85 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(const G4String& aName, G4ProcessType aType) 84 : G4VContinuousDiscreteProce 86 : G4VContinuousDiscreteProcess(aName, aType), 85 lastMaterial(0), 87 lastMaterial(0), 86 nmaxCont1(4), 88 nmaxCont1(4), 87 nmaxCont2(16) 89 nmaxCont2(16) 88 { 90 { 89 } 91 } 90 92 91 // 93 // 92 94 93 G4RDVeLowEnergyLoss::~G4RDVeLowEnergyLoss() 95 G4RDVeLowEnergyLoss::~G4RDVeLowEnergyLoss() 94 { 96 { 95 } 97 } 96 98 97 // 99 // 98 100 99 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(G4RDV 101 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(G4RDVeLowEnergyLoss& right) 100 : G4VContinuousDiscreteProce 102 : G4VContinuousDiscreteProcess(right), 101 lastMaterial(0), 103 lastMaterial(0), 102 nmaxCont1(4), 104 nmaxCont1(4), 103 nmaxCont2(16) 105 nmaxCont2(16) 104 { 106 { 105 } 107 } 106 108 107 void G4RDVeLowEnergyLoss::SetRndmStep(G4bool v 109 void G4RDVeLowEnergyLoss::SetRndmStep(G4bool value) 108 { 110 { 109 rndmStepFlag = value; 111 rndmStepFlag = value; 110 } 112 } 111 113 112 void G4RDVeLowEnergyLoss::SetEnlossFluc(G4bool 114 void G4RDVeLowEnergyLoss::SetEnlossFluc(G4bool value) 113 { 115 { 114 EnlossFlucFlag = value; 116 EnlossFlucFlag = value; 115 } 117 } 116 118 117 void G4RDVeLowEnergyLoss::SetStepFunction (G4d 119 void G4RDVeLowEnergyLoss::SetStepFunction (G4double c1, G4double c2) 118 { 120 { 119 dRoverRange = c1; 121 dRoverRange = c1; 120 finalRange = c2; 122 finalRange = c2; 121 c1lim=dRoverRange; 123 c1lim=dRoverRange; 122 c2lim=2.*(1-dRoverRange)*finalRange; 124 c2lim=2.*(1-dRoverRange)*finalRange; 123 c3lim=-(1.-dRoverRange)*finalRange*finalRan 125 c3lim=-(1.-dRoverRange)*finalRange*finalRange; 124 } 126 } 125 127 126 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang 128 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeTable(G4PhysicsTable* theDEDXTable, 127 G4PhysicsTable* theRangeTable, 129 G4PhysicsTable* theRangeTable, 128 G4double lowestKineticEnergy, 130 G4double lowestKineticEnergy, 129 G4double highestKineticEnergy, 131 G4double highestKineticEnergy, 130 G4int TotBin) 132 G4int TotBin) 131 // Build range table from the energy loss tabl 133 // Build range table from the energy loss table 132 { 134 { 133 135 134 G4int numOfCouples = theDEDXTable->length() 136 G4int numOfCouples = theDEDXTable->length(); 135 137 136 if(theRangeTable) 138 if(theRangeTable) 137 { theRangeTable->clearAndDestroy(); 139 { theRangeTable->clearAndDestroy(); 138 delete theRangeTable; } 140 delete theRangeTable; } 139 theRangeTable = new G4PhysicsTable(numOfCou 141 theRangeTable = new G4PhysicsTable(numOfCouples); 140 142 141 // loop for materials 143 // loop for materials 142 144 143 for (G4int J=0; J<numOfCouples; J++) 145 for (G4int J=0; J<numOfCouples; J++) 144 { 146 { 145 G4PhysicsLogVector* aVector; 147 G4PhysicsLogVector* aVector; 146 aVector = new G4PhysicsLogVector(lowestKi 148 aVector = new G4PhysicsLogVector(lowestKineticEnergy, 147 highestKineticEn 149 highestKineticEnergy,TotBin); 148 BuildRangeVector(theDEDXTable,lowestKinet 150 BuildRangeVector(theDEDXTable,lowestKineticEnergy,highestKineticEnergy, 149 TotBin,J,aVector); 151 TotBin,J,aVector); 150 theRangeTable->insert(aVector); 152 theRangeTable->insert(aVector); 151 } 153 } 152 return theRangeTable ; 154 return theRangeTable ; 153 } 155 } 154 156 155 // 157 // 156 158 157 void G4RDVeLowEnergyLoss::BuildRangeVector(G4P 159 void G4RDVeLowEnergyLoss::BuildRangeVector(G4PhysicsTable* theDEDXTable, 158 G4dou 160 G4double lowestKineticEnergy, 159 G4dou 161 G4double, 160 G4int 162 G4int TotBin, 161 G4int 163 G4int materialIndex, 162 G4Phy 164 G4PhysicsLogVector* rangeVector) 163 165 164 // create range vector for a material 166 // create range vector for a material 165 { 167 { 166 G4bool isOut; 168 G4bool isOut; 167 G4PhysicsVector* physicsVector= (*theDEDXTab 169 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex]; 168 G4double energy1 = lowestKineticEnergy; 170 G4double energy1 = lowestKineticEnergy; 169 G4double dedx = physicsVector->GetValue(e 171 G4double dedx = physicsVector->GetValue(energy1,isOut); 170 G4double range = 0.5*energy1/dedx; 172 G4double range = 0.5*energy1/dedx; 171 rangeVector->PutValue(0,range); 173 rangeVector->PutValue(0,range); 172 G4int n = 100; 174 G4int n = 100; 173 G4double del = 1.0/(G4double)n ; 175 G4double del = 1.0/(G4double)n ; 174 176 175 for (G4int j=1; j<TotBin; j++) { 177 for (G4int j=1; j<TotBin; j++) { 176 178 177 G4double energy2 = rangeVector->GetLowEdge 179 G4double energy2 = rangeVector->GetLowEdgeEnergy(j); 178 G4double de = (energy2 - energy1) * del ; 180 G4double de = (energy2 - energy1) * del ; 179 G4double dedx1 = dedx ; 181 G4double dedx1 = dedx ; 180 182 181 for (G4int i=1; i<n; i++) { 183 for (G4int i=1; i<n; i++) { 182 G4double energy = energy1 + i*de ; 184 G4double energy = energy1 + i*de ; 183 G4double dedx2 = physicsVector->GetValu 185 G4double dedx2 = physicsVector->GetValue(energy,isOut); 184 range += 0.5*de*(1.0/dedx1 + 1.0/dedx2) 186 range += 0.5*de*(1.0/dedx1 + 1.0/dedx2); 185 dedx1 = dedx2; 187 dedx1 = dedx2; 186 } 188 } 187 rangeVector->PutValue(j,range); 189 rangeVector->PutValue(j,range); 188 dedx = dedx1 ; 190 dedx = dedx1 ; 189 energy1 = energy2 ; 191 energy1 = energy2 ; 190 } 192 } 191 } 193 } 192 194 193 // 195 // 194 196 195 G4double G4RDVeLowEnergyLoss::RangeIntLin(G4Ph 197 G4double G4RDVeLowEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector, 196 G4int nbin) 198 G4int nbin) 197 // num. integration, linear binning 199 // num. integration, linear binning 198 { 200 { 199 G4double dtau,Value,taui,ti,lossi,ci; 201 G4double dtau,Value,taui,ti,lossi,ci; 200 G4bool isOut; 202 G4bool isOut; 201 dtau = (tauhigh-taulow)/nbin; 203 dtau = (tauhigh-taulow)/nbin; 202 Value = 0.; 204 Value = 0.; 203 205 204 for (G4int i=0; i<=nbin; i++) 206 for (G4int i=0; i<=nbin; i++) 205 { 207 { 206 taui = taulow + dtau*i ; 208 taui = taulow + dtau*i ; 207 ti = ParticleMass*taui; 209 ti = ParticleMass*taui; 208 lossi = physicsVector->GetValue(ti,isOut); 210 lossi = physicsVector->GetValue(ti,isOut); 209 if(i==0) 211 if(i==0) 210 ci=0.5; 212 ci=0.5; 211 else 213 else 212 { 214 { 213 if(i<nbin) 215 if(i<nbin) 214 ci=1.; 216 ci=1.; 215 else 217 else 216 ci=0.5; 218 ci=0.5; 217 } 219 } 218 Value += ci/lossi; 220 Value += ci/lossi; 219 } 221 } 220 Value *= ParticleMass*dtau; 222 Value *= ParticleMass*dtau; 221 return Value; 223 return Value; 222 } 224 } 223 225 224 // 226 // 225 227 226 G4double G4RDVeLowEnergyLoss::RangeIntLog(G4Ph 228 G4double G4RDVeLowEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector, 227 G4int nbin) 229 G4int nbin) 228 // num. integration, logarithmic binning 230 // num. integration, logarithmic binning 229 { 231 { 230 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci 232 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci; 231 G4bool isOut; 233 G4bool isOut; 232 ltt = ltauhigh-ltaulow; 234 ltt = ltauhigh-ltaulow; 233 dltau = ltt/nbin; 235 dltau = ltt/nbin; 234 Value = 0.; 236 Value = 0.; 235 237 236 for (G4int i=0; i<=nbin; i++) 238 for (G4int i=0; i<=nbin; i++) 237 { 239 { 238 ui = ltaulow+dltau*i; 240 ui = ltaulow+dltau*i; 239 taui = std::exp(ui); 241 taui = std::exp(ui); 240 ti = ParticleMass*taui; 242 ti = ParticleMass*taui; 241 lossi = physicsVector->GetValue(ti,isOut); 243 lossi = physicsVector->GetValue(ti,isOut); 242 if(i==0) 244 if(i==0) 243 ci=0.5; 245 ci=0.5; 244 else 246 else 245 { 247 { 246 if(i<nbin) 248 if(i<nbin) 247 ci=1.; 249 ci=1.; 248 else 250 else 249 ci=0.5; 251 ci=0.5; 250 } 252 } 251 Value += ci*taui/lossi; 253 Value += ci*taui/lossi; 252 } 254 } 253 Value *= ParticleMass*dltau; 255 Value *= ParticleMass*dltau; 254 return Value; 256 return Value; 255 } 257 } 256 258 257 259 258 // 260 // 259 261 260 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildLabT 262 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildLabTimeTable(G4PhysicsTable* theDEDXTable, 261 G4PhysicsTable* theLabTimeTab 263 G4PhysicsTable* theLabTimeTable, 262 G4double lowestKineticEnergy, 264 G4double lowestKineticEnergy, 263 G4double highestKineticEnergy 265 G4double highestKineticEnergy,G4int TotBin) 264 266 265 { 267 { 266 268 267 G4int numOfCouples = G4ProductionCutsTable:: 269 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize(); 268 270 269 if(theLabTimeTable) 271 if(theLabTimeTable) 270 { theLabTimeTable->clearAndDestroy(); 272 { theLabTimeTable->clearAndDestroy(); 271 delete theLabTimeTable; } 273 delete theLabTimeTable; } 272 theLabTimeTable = new G4PhysicsTable(numOfCo 274 theLabTimeTable = new G4PhysicsTable(numOfCouples); 273 275 274 276 275 for (G4int J=0; J<numOfCouples; J++) 277 for (G4int J=0; J<numOfCouples; J++) 276 { 278 { 277 G4PhysicsLogVector* aVector; 279 G4PhysicsLogVector* aVector; 278 280 279 aVector = new G4PhysicsLogVector(lowestKin 281 aVector = new G4PhysicsLogVector(lowestKineticEnergy, 280 highestKineticEner 282 highestKineticEnergy,TotBin); 281 283 282 BuildLabTimeVector(theDEDXTable, 284 BuildLabTimeVector(theDEDXTable, 283 lowestKineticEnergy,highestKinet 285 lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector); 284 theLabTimeTable->insert(aVector); 286 theLabTimeTable->insert(aVector); 285 287 286 288 287 } 289 } 288 return theLabTimeTable ; 290 return theLabTimeTable ; 289 } 291 } 290 292 291 // 293 // 292 294 293 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildProp 295 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildProperTimeTable(G4PhysicsTable* theDEDXTable, 294 G4PhysicsTable* theProperTimeTab 296 G4PhysicsTable* theProperTimeTable, 295 G4double lowestKineticEnergy, 297 G4double lowestKineticEnergy, 296 G4double highestKineticEnergy,G4 298 G4double highestKineticEnergy,G4int TotBin) 297 299 298 { 300 { 299 301 300 G4int numOfCouples = G4ProductionCutsTable:: 302 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize(); 301 303 302 if(theProperTimeTable) 304 if(theProperTimeTable) 303 { theProperTimeTable->clearAndDestroy(); 305 { theProperTimeTable->clearAndDestroy(); 304 delete theProperTimeTable; } 306 delete theProperTimeTable; } 305 theProperTimeTable = new G4PhysicsTable(numO 307 theProperTimeTable = new G4PhysicsTable(numOfCouples); 306 308 307 309 308 for (G4int J=0; J<numOfCouples; J++) 310 for (G4int J=0; J<numOfCouples; J++) 309 { 311 { 310 G4PhysicsLogVector* aVector; 312 G4PhysicsLogVector* aVector; 311 313 312 aVector = new G4PhysicsLogVector(lowestKin 314 aVector = new G4PhysicsLogVector(lowestKineticEnergy, 313 highestKineticEner 315 highestKineticEnergy,TotBin); 314 316 315 BuildProperTimeVector(theDEDXTable, 317 BuildProperTimeVector(theDEDXTable, 316 lowestKineticEnergy,highestKinet 318 lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector); 317 theProperTimeTable->insert(aVector); 319 theProperTimeTable->insert(aVector); 318 320 319 321 320 } 322 } 321 return theProperTimeTable ; 323 return theProperTimeTable ; 322 } 324 } 323 325 324 // 326 // 325 327 326 void G4RDVeLowEnergyLoss::BuildLabTimeVector(G 328 void G4RDVeLowEnergyLoss::BuildLabTimeVector(G4PhysicsTable* theDEDXTable, 327 G4double, // lowestKineticEnergy, 329 G4double, // lowestKineticEnergy, 328 G4double, // highestKineticEnergy 330 G4double, // highestKineticEnergy, 329 G4i 331 G4int TotBin, 330 G4int materialIndex, G4PhysicsLog 332 G4int materialIndex, G4PhysicsLogVector* timeVector) 331 // create lab time vector for a material 333 // create lab time vector for a material 332 { 334 { 333 335 334 G4int nbin=100; 336 G4int nbin=100; 335 G4bool isOut; 337 G4bool isOut; 336 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p 338 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ; 337 G4double losslim,clim,taulim,timelim, 339 G4double losslim,clim,taulim,timelim, 338 LowEdgeEnergy,tau,Value ; 340 LowEdgeEnergy,tau,Value ; 339 341 340 G4PhysicsVector* physicsVector= (*theDEDXTab 342 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex]; 341 343 342 // low energy part first... 344 // low energy part first... 343 losslim = physicsVector->GetValue(tlim,isOut 345 losslim = physicsVector->GetValue(tlim,isOut); 344 taulim=tlim/ParticleMass ; 346 taulim=tlim/ParticleMass ; 345 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh 347 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ; 346 348 347 G4int i=-1; 349 G4int i=-1; 348 G4double oldValue = 0. ; 350 G4double oldValue = 0. ; 349 G4double tauold ; 351 G4double tauold ; 350 do 352 do 351 { 353 { 352 i += 1 ; 354 i += 1 ; 353 LowEdgeEnergy = timeVector->GetLowEdgeEner 355 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i); 354 tau = LowEdgeEnergy/ParticleMass ; 356 tau = LowEdgeEnergy/ParticleMass ; 355 if ( tau <= taulim ) 357 if ( tau <= taulim ) 356 { 358 { 357 Value = clim*std::exp(ppar*std::log(tau/ 359 Value = clim*std::exp(ppar*std::log(tau/taulim)) ; 358 } 360 } 359 else 361 else 360 { 362 { 361 timelim=clim ; 363 timelim=clim ; 362 ltaulow = std::log(taulim); 364 ltaulow = std::log(taulim); 363 ltauhigh = std::log(tau); 365 ltauhigh = std::log(tau); 364 Value = timelim+LabTimeIntLog(physicsVec 366 Value = timelim+LabTimeIntLog(physicsVector,nbin); 365 } 367 } 366 timeVector->PutValue(i,Value); 368 timeVector->PutValue(i,Value); 367 oldValue = Value ; 369 oldValue = Value ; 368 tauold = tau ; 370 tauold = tau ; 369 } while (tau<=taulim) ; 371 } while (tau<=taulim) ; 370 i += 1 ; 372 i += 1 ; 371 for (G4int j=i; j<TotBin; j++) 373 for (G4int j=i; j<TotBin; j++) 372 { 374 { 373 LowEdgeEnergy = timeVector->GetLowEdgeEner 375 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j); 374 tau = LowEdgeEnergy/ParticleMass ; 376 tau = LowEdgeEnergy/ParticleMass ; 375 ltaulow = std::log(tauold); 377 ltaulow = std::log(tauold); 376 ltauhigh = std::log(tau); 378 ltauhigh = std::log(tau); 377 Value = oldValue+LabTimeIntLog(physicsVect 379 Value = oldValue+LabTimeIntLog(physicsVector,nbin); 378 timeVector->PutValue(j,Value); 380 timeVector->PutValue(j,Value); 379 oldValue = Value ; 381 oldValue = Value ; 380 tauold = tau ; 382 tauold = tau ; 381 } 383 } 382 } 384 } 383 385 384 // 386 // 385 387 386 void G4RDVeLowEnergyLoss::BuildProperTimeVecto 388 void G4RDVeLowEnergyLoss::BuildProperTimeVector(G4PhysicsTable* theDEDXTable, 387 G4double, // lowestKineticEner 389 G4double, // lowestKineticEnergy, 388 G4double, // highestKineticEne 390 G4double, // highestKineticEnergy, 389 391 G4int TotBin, 390 G4int materialIndex, G4Physics 392 G4int materialIndex, G4PhysicsLogVector* timeVector) 391 // create proper time vector for a material 393 // create proper time vector for a material 392 { 394 { 393 G4int nbin=100; 395 G4int nbin=100; 394 G4bool isOut; 396 G4bool isOut; 395 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p 397 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ; 396 G4double losslim,clim,taulim,timelim, 398 G4double losslim,clim,taulim,timelim, 397 LowEdgeEnergy,tau,Value ; 399 LowEdgeEnergy,tau,Value ; 398 400 399 G4PhysicsVector* physicsVector= (*theDEDXTab 401 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex]; 400 //const G4MaterialTable* theMaterialTable = 402 //const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 401 403 402 // low energy part first... 404 // low energy part first... 403 losslim = physicsVector->GetValue(tlim,isOut 405 losslim = physicsVector->GetValue(tlim,isOut); 404 taulim=tlim/ParticleMass ; 406 taulim=tlim/ParticleMass ; 405 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh 407 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ; 406 408 407 G4int i=-1; 409 G4int i=-1; 408 G4double oldValue = 0. ; 410 G4double oldValue = 0. ; 409 G4double tauold ; 411 G4double tauold ; 410 do 412 do 411 { 413 { 412 i += 1 ; 414 i += 1 ; 413 LowEdgeEnergy = timeVector->GetLowEdgeEner 415 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i); 414 tau = LowEdgeEnergy/ParticleMass ; 416 tau = LowEdgeEnergy/ParticleMass ; 415 if ( tau <= taulim ) 417 if ( tau <= taulim ) 416 { 418 { 417 Value = clim*std::exp(ppar*std::log(tau/ 419 Value = clim*std::exp(ppar*std::log(tau/taulim)) ; 418 } 420 } 419 else 421 else 420 { 422 { 421 timelim=clim ; 423 timelim=clim ; 422 ltaulow = std::log(taulim); 424 ltaulow = std::log(taulim); 423 ltauhigh = std::log(tau); 425 ltauhigh = std::log(tau); 424 Value = timelim+ProperTimeIntLog(physics 426 Value = timelim+ProperTimeIntLog(physicsVector,nbin); 425 } 427 } 426 timeVector->PutValue(i,Value); 428 timeVector->PutValue(i,Value); 427 oldValue = Value ; 429 oldValue = Value ; 428 tauold = tau ; 430 tauold = tau ; 429 } while (tau<=taulim) ; 431 } while (tau<=taulim) ; 430 i += 1 ; 432 i += 1 ; 431 for (G4int j=i; j<TotBin; j++) 433 for (G4int j=i; j<TotBin; j++) 432 { 434 { 433 LowEdgeEnergy = timeVector->GetLowEdgeEner 435 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j); 434 tau = LowEdgeEnergy/ParticleMass ; 436 tau = LowEdgeEnergy/ParticleMass ; 435 ltaulow = std::log(tauold); 437 ltaulow = std::log(tauold); 436 ltauhigh = std::log(tau); 438 ltauhigh = std::log(tau); 437 Value = oldValue+ProperTimeIntLog(physicsV 439 Value = oldValue+ProperTimeIntLog(physicsVector,nbin); 438 timeVector->PutValue(j,Value); 440 timeVector->PutValue(j,Value); 439 oldValue = Value ; 441 oldValue = Value ; 440 tauold = tau ; 442 tauold = tau ; 441 } 443 } 442 } 444 } 443 445 444 // 446 // 445 447 446 G4double G4RDVeLowEnergyLoss::LabTimeIntLog(G4 448 G4double G4RDVeLowEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector, 447 G4int nbin) 449 G4int nbin) 448 // num. integration, logarithmic binning 450 // num. integration, logarithmic binning 449 { 451 { 450 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci 452 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci; 451 G4bool isOut; 453 G4bool isOut; 452 ltt = ltauhigh-ltaulow; 454 ltt = ltauhigh-ltaulow; 453 dltau = ltt/nbin; 455 dltau = ltt/nbin; 454 Value = 0.; 456 Value = 0.; 455 457 456 for (G4int i=0; i<=nbin; i++) 458 for (G4int i=0; i<=nbin; i++) 457 { 459 { 458 ui = ltaulow+dltau*i; 460 ui = ltaulow+dltau*i; 459 taui = std::exp(ui); 461 taui = std::exp(ui); 460 ti = ParticleMass*taui; 462 ti = ParticleMass*taui; 461 lossi = physicsVector->GetValue(ti,isOut); 463 lossi = physicsVector->GetValue(ti,isOut); 462 if(i==0) 464 if(i==0) 463 ci=0.5; 465 ci=0.5; 464 else 466 else 465 { 467 { 466 if(i<nbin) 468 if(i<nbin) 467 ci=1.; 469 ci=1.; 468 else 470 else 469 ci=0.5; 471 ci=0.5; 470 } 472 } 471 Value += ci*taui*(ti+ParticleMass)/(std::s 473 Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi); 472 } 474 } 473 Value *= ParticleMass*dltau/c_light; 475 Value *= ParticleMass*dltau/c_light; 474 return Value; 476 return Value; 475 } 477 } 476 478 477 // 479 // 478 480 479 G4double G4RDVeLowEnergyLoss::ProperTimeIntLog 481 G4double G4RDVeLowEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector, 480 G4int nbin) 482 G4int nbin) 481 // num. integration, logarithmic binning 483 // num. integration, logarithmic binning 482 { 484 { 483 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci 485 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci; 484 G4bool isOut; 486 G4bool isOut; 485 ltt = ltauhigh-ltaulow; 487 ltt = ltauhigh-ltaulow; 486 dltau = ltt/nbin; 488 dltau = ltt/nbin; 487 Value = 0.; 489 Value = 0.; 488 490 489 for (G4int i=0; i<=nbin; i++) 491 for (G4int i=0; i<=nbin; i++) 490 { 492 { 491 ui = ltaulow+dltau*i; 493 ui = ltaulow+dltau*i; 492 taui = std::exp(ui); 494 taui = std::exp(ui); 493 ti = ParticleMass*taui; 495 ti = ParticleMass*taui; 494 lossi = physicsVector->GetValue(ti,isOut); 496 lossi = physicsVector->GetValue(ti,isOut); 495 if(i==0) 497 if(i==0) 496 ci=0.5; 498 ci=0.5; 497 else 499 else 498 { 500 { 499 if(i<nbin) 501 if(i<nbin) 500 ci=1.; 502 ci=1.; 501 else 503 else 502 ci=0.5; 504 ci=0.5; 503 } 505 } 504 Value += ci*taui*ParticleMass/(std::sqrt(t 506 Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi); 505 } 507 } 506 Value *= ParticleMass*dltau/c_light; 508 Value *= ParticleMass*dltau/c_light; 507 return Value; 509 return Value; 508 } 510 } 509 511 510 // 512 // 511 513 512 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildInve 514 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildInverseRangeTable(G4PhysicsTable* theRangeTable, 513 G4PhysicsTable*, 515 G4PhysicsTable*, 514 G4PhysicsTable*, 516 G4PhysicsTable*, 515 G4PhysicsTable*, 517 G4PhysicsTable*, 516 G4PhysicsTable* theInverseRang 518 G4PhysicsTable* theInverseRangeTable, 517 G4double, // lowestKineticEner 519 G4double, // lowestKineticEnergy, 518 G4double, // highestKineticEne 520 G4double, // highestKineticEnergy 519 G4int ) // nbins 521 G4int ) // nbins 520 // Build inverse table of the range table 522 // Build inverse table of the range table 521 { 523 { 522 G4bool b; 524 G4bool b; 523 525 524 G4int numOfCouples = G4ProductionCutsTable:: 526 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize(); 525 527 526 if(theInverseRangeTable) 528 if(theInverseRangeTable) 527 { theInverseRangeTable->clearAndDestroy(); 529 { theInverseRangeTable->clearAndDestroy(); 528 delete theInverseRangeTable; } 530 delete theInverseRangeTable; } 529 theInverseRangeTable = new G4PhysicsTable( 531 theInverseRangeTable = new G4PhysicsTable(numOfCouples); 530 532 531 // loop for materials 533 // loop for materials 532 for (G4int i=0; i<numOfCouples; i++) 534 for (G4int i=0; i<numOfCouples; i++) 533 { 535 { 534 536 535 G4PhysicsVector* pv = (*theRangeTable)[i]; 537 G4PhysicsVector* pv = (*theRangeTable)[i]; 536 size_t nbins = pv->GetVectorLength(); 538 size_t nbins = pv->GetVectorLength(); 537 G4double elow = pv->GetLowEdgeEnergy(0); 539 G4double elow = pv->GetLowEdgeEnergy(0); 538 G4double ehigh = pv->GetLowEdgeEnergy(nbin 540 G4double ehigh = pv->GetLowEdgeEnergy(nbins-1); 539 G4double rlow = pv->GetValue(elow, b); 541 G4double rlow = pv->GetValue(elow, b); 540 G4double rhigh = pv->GetValue(ehigh, b); 542 G4double rhigh = pv->GetValue(ehigh, b); 541 543 542 rhigh *= std::exp(std::log(rhigh/rlow)/((G 544 rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1))); 543 545 544 G4PhysicsLogVector* v = new G4PhysicsLogVe 546 G4PhysicsLogVector* v = new G4PhysicsLogVector(rlow, rhigh, nbins); 545 547 546 v->PutValue(0,elow); 548 v->PutValue(0,elow); 547 G4double energy1 = elow; 549 G4double energy1 = elow; 548 G4double range1 = rlow; 550 G4double range1 = rlow; 549 G4double energy2 = elow; 551 G4double energy2 = elow; 550 G4double range2 = rlow; 552 G4double range2 = rlow; 551 size_t ilow = 0; 553 size_t ilow = 0; 552 size_t ihigh; 554 size_t ihigh; 553 555 554 for (size_t j=1; j<nbins; j++) { 556 for (size_t j=1; j<nbins; j++) { 555 557 556 G4double range = v->GetLowEdgeEnergy(j); 558 G4double range = v->GetLowEdgeEnergy(j); 557 559 558 for (ihigh=ilow+1; ihigh<nbins; ihigh++) 560 for (ihigh=ilow+1; ihigh<nbins; ihigh++) { 559 energy2 = pv->GetLowEdgeEnergy(ihigh); 561 energy2 = pv->GetLowEdgeEnergy(ihigh); 560 range2 = pv->GetValue(energy2, b); 562 range2 = pv->GetValue(energy2, b); 561 if(range2 >= range || ihigh == nbins-1 563 if(range2 >= range || ihigh == nbins-1) { 562 ilow = ihigh - 1; 564 ilow = ihigh - 1; 563 energy1 = pv->GetLowEdgeEnergy(ilow) 565 energy1 = pv->GetLowEdgeEnergy(ilow); 564 range1 = pv->GetValue(energy1, b); 566 range1 = pv->GetValue(energy1, b); 565 break; 567 break; 566 } 568 } 567 } 569 } 568 570 569 G4double e = std::log(energy1) + std::lo 571 G4double e = std::log(energy1) + std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1); 570 572 571 v->PutValue(j,std::exp(e)); 573 v->PutValue(j,std::exp(e)); 572 } 574 } 573 theInverseRangeTable->insert(v); 575 theInverseRangeTable->insert(v); 574 576 575 } 577 } 576 return theInverseRangeTable ; 578 return theInverseRangeTable ; 577 } 579 } 578 580 579 // 581 // 580 582 581 void G4RDVeLowEnergyLoss::InvertRangeVector(G4 583 void G4RDVeLowEnergyLoss::InvertRangeVector(G4PhysicsTable* theRangeTable, 582 G4PhysicsTable* theRangeCoeffATabl 584 G4PhysicsTable* theRangeCoeffATable, 583 G4PhysicsTable* theRangeCoeffBTabl 585 G4PhysicsTable* theRangeCoeffBTable, 584 G4PhysicsTable* theRangeCoeffCTabl 586 G4PhysicsTable* theRangeCoeffCTable, 585 G4double lowestKineticEnergy, 587 G4double lowestKineticEnergy, 586 G4double highestKineticEnergy, G4i 588 G4double highestKineticEnergy, G4int TotBin, 587 G4int materialIndex, G4PhysicsLog 589 G4int materialIndex, G4PhysicsLogVector* aVector) 588 // invert range vector for a material 590 // invert range vector for a material 589 { 591 { 590 G4double LowEdgeRange,A,B,C,discr,KineticEne 592 G4double LowEdgeRange,A,B,C,discr,KineticEnergy ; 591 G4double RTable = std::exp(std::log(highestK 593 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ; 592 G4double Tbin = lowestKineticEnergy/RTable ; 594 G4double Tbin = lowestKineticEnergy/RTable ; 593 G4double rangebin = 0.0 ; 595 G4double rangebin = 0.0 ; 594 G4int binnumber = -1 ; 596 G4int binnumber = -1 ; 595 G4bool isOut ; 597 G4bool isOut ; 596 598 597 //loop for range values 599 //loop for range values 598 for( G4int i=0; i<TotBin; i++) 600 for( G4int i=0; i<TotBin; i++) 599 { 601 { 600 LowEdgeRange = aVector->GetLowEdgeEnergy(i 602 LowEdgeRange = aVector->GetLowEdgeEnergy(i) ; //i.e. GetLowEdgeValue(i) 601 if( rangebin < LowEdgeRange ) 603 if( rangebin < LowEdgeRange ) 602 { 604 { 603 do 605 do 604 { 606 { 605 binnumber += 1 ; 607 binnumber += 1 ; 606 Tbin *= RTable ; 608 Tbin *= RTable ; 607 rangebin = (*theRangeTable)(materialIn 609 rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ; 608 } 610 } 609 while ((rangebin < LowEdgeRange) && (bin 611 while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ; 610 } 612 } 611 613 612 if(binnumber == 0) 614 if(binnumber == 0) 613 KineticEnergy = lowestKineticEnergy ; 615 KineticEnergy = lowestKineticEnergy ; 614 else if(binnumber == TotBin-1) 616 else if(binnumber == TotBin-1) 615 KineticEnergy = highestKineticEnergy ; 617 KineticEnergy = highestKineticEnergy ; 616 else 618 else 617 { 619 { 618 A = (*(*theRangeCoeffATable)(materialInd 620 A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ; 619 B = (*(*theRangeCoeffBTable)(materialInd 621 B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ; 620 C = (*(*theRangeCoeffCTable)(materialInd 622 C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ; 621 if(A==0.) 623 if(A==0.) 622 KineticEnergy = (LowEdgeRange -C )/B 624 KineticEnergy = (LowEdgeRange -C )/B ; 623 else 625 else 624 { 626 { 625 discr = B*B - 4.*A*(C-LowEdgeRange); 627 discr = B*B - 4.*A*(C-LowEdgeRange); 626 discr = discr>0. ? std::sqrt(discr) : 628 discr = discr>0. ? std::sqrt(discr) : 0.; 627 KineticEnergy = 0.5*(discr-B)/A ; 629 KineticEnergy = 0.5*(discr-B)/A ; 628 } 630 } 629 } 631 } 630 632 631 aVector->PutValue(i,KineticEnergy) ; 633 aVector->PutValue(i,KineticEnergy) ; 632 } 634 } 633 } 635 } 634 636 635 // 637 // 636 638 637 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang 639 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffATable(G4PhysicsTable* theRangeTable, 638 G4PhysicsTable* theRangeCoeffAT 640 G4PhysicsTable* theRangeCoeffATable, 639 G4double lowestKineticEnergy, 641 G4double lowestKineticEnergy, 640 G4double highestKineticEnergy, 642 G4double highestKineticEnergy, G4int TotBin) 641 // Build tables of coefficients for the energy 643 // Build tables of coefficients for the energy loss calculation 642 // create table for coefficients "A" 644 // create table for coefficients "A" 643 { 645 { 644 646 645 G4int numOfCouples = G4ProductionCutsTable:: 647 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize(); 646 648 647 if(theRangeCoeffATable) 649 if(theRangeCoeffATable) 648 { theRangeCoeffATable->clearAndDestroy(); 650 { theRangeCoeffATable->clearAndDestroy(); 649 delete theRangeCoeffATable; } 651 delete theRangeCoeffATable; } 650 theRangeCoeffATable = new G4PhysicsTable(num 652 theRangeCoeffATable = new G4PhysicsTable(numOfCouples); 651 653 652 G4double RTable = std::exp(std::log(highestK 654 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ; 653 G4double R2 = RTable*RTable ; 655 G4double R2 = RTable*RTable ; 654 G4double R1 = RTable+1.; 656 G4double R1 = RTable+1.; 655 G4double w = R1*(RTable-1.)*(RTable-1.); 657 G4double w = R1*(RTable-1.)*(RTable-1.); 656 G4double w1 = RTable/w , w2 = -RTable*R1/w , 658 G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ; 657 G4double Ti , Tim , Tip , Ri , Rim , Rip , V 659 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ; 658 G4bool isOut; 660 G4bool isOut; 659 661 660 // loop for materials 662 // loop for materials 661 for (G4int J=0; J<numOfCouples; J++) 663 for (G4int J=0; J<numOfCouples; J++) 662 { 664 { 663 G4int binmax=TotBin ; 665 G4int binmax=TotBin ; 664 G4PhysicsLinearVector* aVector = 666 G4PhysicsLinearVector* aVector = 665 new G4PhysicsLinear 667 new G4PhysicsLinearVector(0.,binmax, TotBin); 666 Ti = lowestKineticEnergy ; 668 Ti = lowestKineticEnergy ; 667 G4PhysicsVector* rangeVector= (*theRangeTa 669 G4PhysicsVector* rangeVector= (*theRangeTable)[J]; 668 670 669 for ( G4int i=0; i<TotBin; i++) 671 for ( G4int i=0; i<TotBin; i++) 670 { 672 { 671 Ri = rangeVector->GetValue(Ti,isOut) ; 673 Ri = rangeVector->GetValue(Ti,isOut) ; 672 if ( i==0 ) 674 if ( i==0 ) 673 Rim = 0. ; 675 Rim = 0. ; 674 else 676 else 675 { 677 { 676 Tim = Ti/RTable ; 678 Tim = Ti/RTable ; 677 Rim = rangeVector->GetValue(Tim,isOut) 679 Rim = rangeVector->GetValue(Tim,isOut); 678 } 680 } 679 if ( i==(TotBin-1)) 681 if ( i==(TotBin-1)) 680 Rip = Ri ; 682 Rip = Ri ; 681 else 683 else 682 { 684 { 683 Tip = Ti*RTable ; 685 Tip = Ti*RTable ; 684 Rip = rangeVector->GetValue(Tip,isOut) 686 Rip = rangeVector->GetValue(Tip,isOut); 685 } 687 } 686 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti 688 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ; 687 689 688 aVector->PutValue(i,Value); 690 aVector->PutValue(i,Value); 689 Ti = RTable*Ti ; 691 Ti = RTable*Ti ; 690 } 692 } 691 693 692 theRangeCoeffATable->insert(aVector); 694 theRangeCoeffATable->insert(aVector); 693 } 695 } 694 return theRangeCoeffATable ; 696 return theRangeCoeffATable ; 695 } 697 } 696 698 697 // 699 // 698 700 699 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang 701 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffBTable(G4PhysicsTable* theRangeTable, 700 G4PhysicsTable* theRangeCoeffBT 702 G4PhysicsTable* theRangeCoeffBTable, 701 G4double lowestKineticEnergy, 703 G4double lowestKineticEnergy, 702 G4double highestKineticEnergy, 704 G4double highestKineticEnergy, G4int TotBin) 703 // Build tables of coefficients for the energy 705 // Build tables of coefficients for the energy loss calculation 704 // create table for coefficients "B" 706 // create table for coefficients "B" 705 { 707 { 706 708 707 G4int numOfCouples = G4ProductionCutsTable:: 709 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize(); 708 710 709 if(theRangeCoeffBTable) 711 if(theRangeCoeffBTable) 710 { theRangeCoeffBTable->clearAndDestroy(); 712 { theRangeCoeffBTable->clearAndDestroy(); 711 delete theRangeCoeffBTable; } 713 delete theRangeCoeffBTable; } 712 theRangeCoeffBTable = new G4PhysicsTable(num 714 theRangeCoeffBTable = new G4PhysicsTable(numOfCouples); 713 715 714 G4double RTable = std::exp(std::log(highestK 716 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ; 715 G4double R2 = RTable*RTable ; 717 G4double R2 = RTable*RTable ; 716 G4double R1 = RTable+1.; 718 G4double R1 = RTable+1.; 717 G4double w = R1*(RTable-1.)*(RTable-1.); 719 G4double w = R1*(RTable-1.)*(RTable-1.); 718 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 720 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 = -R2*R1/w ; 719 G4double Ti , Tim , Tip , Ri , Rim , Rip , V 721 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ; 720 G4bool isOut; 722 G4bool isOut; 721 723 722 // loop for materials 724 // loop for materials 723 for (G4int J=0; J<numOfCouples; J++) 725 for (G4int J=0; J<numOfCouples; J++) 724 { 726 { 725 G4int binmax=TotBin ; 727 G4int binmax=TotBin ; 726 G4PhysicsLinearVector* aVector = 728 G4PhysicsLinearVector* aVector = 727 new G4PhysicsLinearVec 729 new G4PhysicsLinearVector(0.,binmax, TotBin); 728 Ti = lowestKineticEnergy ; 730 Ti = lowestKineticEnergy ; 729 G4PhysicsVector* rangeVector= (*theRangeTa 731 G4PhysicsVector* rangeVector= (*theRangeTable)[J]; 730 732 731 for ( G4int i=0; i<TotBin; i++) 733 for ( G4int i=0; i<TotBin; i++) 732 { 734 { 733 Ri = rangeVector->GetValue(Ti,isOut) ; 735 Ri = rangeVector->GetValue(Ti,isOut) ; 734 if ( i==0 ) 736 if ( i==0 ) 735 Rim = 0. ; 737 Rim = 0. ; 736 else 738 else 737 { 739 { 738 Tim = Ti/RTable ; 740 Tim = Ti/RTable ; 739 Rim = rangeVector->GetValue(Tim,isOut) 741 Rim = rangeVector->GetValue(Tim,isOut); 740 } 742 } 741 if ( i==(TotBin-1)) 743 if ( i==(TotBin-1)) 742 Rip = Ri ; 744 Rip = Ri ; 743 else 745 else 744 { 746 { 745 Tip = Ti*RTable ; 747 Tip = Ti*RTable ; 746 Rip = rangeVector->GetValue(Tip,isOut) 748 Rip = rangeVector->GetValue(Tip,isOut); 747 } 749 } 748 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti; 750 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti; 749 751 750 aVector->PutValue(i,Value); 752 aVector->PutValue(i,Value); 751 Ti = RTable*Ti ; 753 Ti = RTable*Ti ; 752 } 754 } 753 theRangeCoeffBTable->insert(aVector); 755 theRangeCoeffBTable->insert(aVector); 754 } 756 } 755 return theRangeCoeffBTable ; 757 return theRangeCoeffBTable ; 756 } 758 } 757 759 758 // 760 // 759 761 760 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang 762 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffCTable(G4PhysicsTable* theRangeTable, 761 G4PhysicsTable* theRangeCoeffCT 763 G4PhysicsTable* theRangeCoeffCTable, 762 G4double lowestKineticEnergy, 764 G4double lowestKineticEnergy, 763 G4double highestKineticEnergy, 765 G4double highestKineticEnergy, G4int TotBin) 764 // Build tables of coefficients for the energy 766 // Build tables of coefficients for the energy loss calculation 765 // create table for coefficients "C" 767 // create table for coefficients "C" 766 { 768 { 767 769 768 G4int numOfCouples = G4ProductionCutsTable:: 770 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize(); 769 771 770 if(theRangeCoeffCTable) 772 if(theRangeCoeffCTable) 771 { theRangeCoeffCTable->clearAndDestroy(); 773 { theRangeCoeffCTable->clearAndDestroy(); 772 delete theRangeCoeffCTable; } 774 delete theRangeCoeffCTable; } 773 theRangeCoeffCTable = new G4PhysicsTable(num 775 theRangeCoeffCTable = new G4PhysicsTable(numOfCouples); 774 776 775 G4double RTable = std::exp(std::log(highestK 777 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ; 776 G4double R2 = RTable*RTable ; 778 G4double R2 = RTable*RTable ; 777 G4double R1 = RTable+1.; 779 G4double R1 = RTable+1.; 778 G4double w = R1*(RTable-1.)*(RTable-1.); 780 G4double w = R1*(RTable-1.)*(RTable-1.); 779 G4double w1 = 1./w , w2 = -RTable*R1/w , w3 781 G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ; 780 G4double Ti , Tim , Tip , Ri , Rim , Rip , V 782 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ; 781 G4bool isOut; 783 G4bool isOut; 782 784 783 // loop for materials 785 // loop for materials 784 for (G4int J=0; J<numOfCouples; J++) 786 for (G4int J=0; J<numOfCouples; J++) 785 { 787 { 786 G4int binmax=TotBin ; 788 G4int binmax=TotBin ; 787 G4PhysicsLinearVector* aVector = 789 G4PhysicsLinearVector* aVector = 788 new G4PhysicsLinearVecto 790 new G4PhysicsLinearVector(0.,binmax, TotBin); 789 Ti = lowestKineticEnergy ; 791 Ti = lowestKineticEnergy ; 790 G4PhysicsVector* rangeVector= (*theRangeTa 792 G4PhysicsVector* rangeVector= (*theRangeTable)[J]; 791 793 792 for ( G4int i=0; i<TotBin; i++) 794 for ( G4int i=0; i<TotBin; i++) 793 { 795 { 794 Ri = rangeVector->GetValue(Ti,isOut) ; 796 Ri = rangeVector->GetValue(Ti,isOut) ; 795 if ( i==0 ) 797 if ( i==0 ) 796 Rim = 0. ; 798 Rim = 0. ; 797 else 799 else 798 { 800 { 799 Tim = Ti/RTable ; 801 Tim = Ti/RTable ; 800 Rim = rangeVector->GetValue(Tim,isOut) 802 Rim = rangeVector->GetValue(Tim,isOut); 801 } 803 } 802 if ( i==(TotBin-1)) 804 if ( i==(TotBin-1)) 803 Rip = Ri ; 805 Rip = Ri ; 804 else 806 else 805 { 807 { 806 Tip = Ti*RTable ; 808 Tip = Ti*RTable ; 807 Rip = rangeVector->GetValue(Tip,isOut) 809 Rip = rangeVector->GetValue(Tip,isOut); 808 } 810 } 809 Value = w1*Rip + w2*Ri + w3*Rim ; 811 Value = w1*Rip + w2*Ri + w3*Rim ; 810 812 811 aVector->PutValue(i,Value); 813 aVector->PutValue(i,Value); 812 Ti = RTable*Ti ; 814 Ti = RTable*Ti ; 813 } 815 } 814 theRangeCoeffCTable->insert(aVector); 816 theRangeCoeffCTable->insert(aVector); 815 } 817 } 816 return theRangeCoeffCTable ; 818 return theRangeCoeffCTable ; 817 } 819 } 818 820 819 // 821 // 820 822 821 G4double G4RDVeLowEnergyLoss::GetLossWithFluct 823 G4double G4RDVeLowEnergyLoss::GetLossWithFluct(const G4DynamicParticle* aParticle, 822 c 824 const G4MaterialCutsCouple* couple, 823 G4double MeanLoss, 825 G4double MeanLoss, 824 G4double step) 826 G4double step) 825 // calculate actual loss from the mean loss 827 // calculate actual loss from the mean loss 826 // The model used to get the fluctuation is e 828 // The model used to get the fluctuation is essentially the same as in Glandz in Geant3. 827 { 829 { 828 static const G4double minLoss = 1.*eV ; 830 static const G4double minLoss = 1.*eV ; 829 static const G4double probLim = 0.01 ; 831 static const G4double probLim = 0.01 ; 830 static const G4double sumaLim = -std::log(p 832 static const G4double sumaLim = -std::log(probLim) ; 831 static const G4double alim=10.; 833 static const G4double alim=10.; 832 static const G4double kappa = 10. ; 834 static const G4double kappa = 10. ; 833 static const G4double factor = twopi_mc2_rc 835 static const G4double factor = twopi_mc2_rcl2 ; 834 const G4Material* aMaterial = couple->GetMat 836 const G4Material* aMaterial = couple->GetMaterial(); 835 837 836 // check if the material has changed ( cache 838 // check if the material has changed ( cache mechanism) 837 839 838 if (aMaterial != lastMaterial) 840 if (aMaterial != lastMaterial) 839 { 841 { 840 lastMaterial = aMaterial; 842 lastMaterial = aMaterial; 841 imat = couple->GetIndex(); 843 imat = couple->GetIndex(); 842 f1Fluct = aMaterial->GetIonisation( 844 f1Fluct = aMaterial->GetIonisation()->GetF1fluct(); 843 f2Fluct = aMaterial->GetIonisation( 845 f2Fluct = aMaterial->GetIonisation()->GetF2fluct(); 844 e1Fluct = aMaterial->GetIonisation( 846 e1Fluct = aMaterial->GetIonisation()->GetEnergy1fluct(); 845 e2Fluct = aMaterial->GetIonisation( 847 e2Fluct = aMaterial->GetIonisation()->GetEnergy2fluct(); 846 e1LogFluct = aMaterial->GetIonisation( 848 e1LogFluct = aMaterial->GetIonisation()->GetLogEnergy1fluct(); 847 e2LogFluct = aMaterial->GetIonisation( 849 e2LogFluct = aMaterial->GetIonisation()->GetLogEnergy2fluct(); 848 rateFluct = aMaterial->GetIonisation( 850 rateFluct = aMaterial->GetIonisation()->GetRateionexcfluct(); 849 ipotFluct = aMaterial->GetIonisation( 851 ipotFluct = aMaterial->GetIonisation()->GetMeanExcitationEnergy(); 850 ipotLogFluct = aMaterial->GetIonisation( 852 ipotLogFluct = aMaterial->GetIonisation()->GetLogMeanExcEnergy(); 851 } 853 } 852 G4double threshold,w1,w2,C, 854 G4double threshold,w1,w2,C, 853 beta2,suma,e0,loss,lossc,w; 855 beta2,suma,e0,loss,lossc,w; 854 G4double a1,a2,a3; 856 G4double a1,a2,a3; 855 G4int p1,p2,p3; 857 G4int p1,p2,p3; 856 G4int nb; 858 G4int nb; 857 G4double Corrfac, na,alfa,rfac,namean,sa,alf 859 G4double Corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea; 858 // G4double dp1; 860 // G4double dp1; 859 G4double dp3; 861 G4double dp3; 860 G4double siga ; 862 G4double siga ; 861 863 862 // shortcut for very very small loss 864 // shortcut for very very small loss 863 if(MeanLoss < minLoss) return MeanLoss ; 865 if(MeanLoss < minLoss) return MeanLoss ; 864 866 865 // get particle data 867 // get particle data 866 G4double Tkin = aParticle->GetKineticEnerg 868 G4double Tkin = aParticle->GetKineticEnergy(); 867 869 868 // G4cout << "MGP -- Fluc Tkin " << Tkin/ke 870 // G4cout << "MGP -- Fluc Tkin " << Tkin/keV << " keV " << " MeanLoss = " << MeanLoss/keV << G4endl; 869 871 870 threshold = (*((G4ProductionCutsTable::GetPr 872 threshold = (*((G4ProductionCutsTable::GetProductionCutsTable()) 871 ->GetEnergyCutsVector(1)))[ima 873 ->GetEnergyCutsVector(1)))[imat]; 872 G4double rmass = electron_mass_c2/ParticleMa 874 G4double rmass = electron_mass_c2/ParticleMass; 873 G4double tau = Tkin/ParticleMass, tau1 = t 875 G4double tau = Tkin/ParticleMass, tau1 = tau+1., tau2 = tau*(tau+2.); 874 G4double Tm = 2.*electron_mass_c2*tau2/(1 876 G4double Tm = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass); 875 877 876 // G4cout << "MGP Particle mass " << Particl 878 // G4cout << "MGP Particle mass " << ParticleMass/MeV << " Tm " << Tm << G4endl; 877 879 878 if(Tm > threshold) Tm = threshold; 880 if(Tm > threshold) Tm = threshold; 879 beta2 = tau2/(tau1*tau1); 881 beta2 = tau2/(tau1*tau1); 880 882 881 // Gaussian fluctuation ? 883 // Gaussian fluctuation ? 882 if(MeanLoss >= kappa*Tm || MeanLoss <= kappa 884 if(MeanLoss >= kappa*Tm || MeanLoss <= kappa*ipotFluct) 883 { 885 { 884 G4double electronDensity = aMaterial->GetE 886 G4double electronDensity = aMaterial->GetElectronDensity() ; 885 siga = std::sqrt(Tm*(1.0-0.5*beta2)*step* 887 siga = std::sqrt(Tm*(1.0-0.5*beta2)*step* 886 factor*electronDensity/beta2) 888 factor*electronDensity/beta2) ; 887 do { 889 do { 888 loss = G4RandGauss::shoot(MeanLoss,siga) 890 loss = G4RandGauss::shoot(MeanLoss,siga) ; 889 } while (loss < 0. || loss > 2.0*MeanLoss) 891 } while (loss < 0. || loss > 2.0*MeanLoss); 890 return loss ; 892 return loss ; 891 } 893 } 892 894 893 w1 = Tm/ipotFluct; 895 w1 = Tm/ipotFluct; 894 w2 = std::log(2.*electron_mass_c2*tau2); 896 w2 = std::log(2.*electron_mass_c2*tau2); 895 897 896 C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct 898 C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2); 897 899 898 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct 900 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct; 899 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct 901 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct; 900 a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipot 902 a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipotFluct*Tm*std::log(w1)); 901 903 902 suma = a1+a2+a3; 904 suma = a1+a2+a3; 903 905 904 loss = 0. ; 906 loss = 0. ; 905 907 906 if(suma < sumaLim) // very small 908 if(suma < sumaLim) // very small Step 907 { 909 { 908 e0 = aMaterial->GetIonisation()->GetEner 910 e0 = aMaterial->GetIonisation()->GetEnergy0fluct(); 909 // G4cout << "MGP e0 = " << e0/keV << G4 911 // G4cout << "MGP e0 = " << e0/keV << G4endl; 910 912 911 if(Tm == ipotFluct) 913 if(Tm == ipotFluct) 912 { 914 { 913 a3 = MeanLoss/e0; 915 a3 = MeanLoss/e0; 914 916 915 if(a3>alim) 917 if(a3>alim) 916 { 918 { 917 siga=std::sqrt(a3) ; 919 siga=std::sqrt(a3) ; 918 p3 = std::max(0,G4int(G4RandGauss::sho 920 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5)); 919 } 921 } 920 else p3 = G4Poisson(a3); 922 else p3 = G4Poisson(a3); 921 923 922 loss = p3*e0 ; 924 loss = p3*e0 ; 923 925 924 if(p3 > 0) loss += (1.-2.*G4UniformRand()) 926 if(p3 > 0) loss += (1.-2.*G4UniformRand())*e0 ; 925 // G4cout << "MGP very small step " << los 927 // G4cout << "MGP very small step " << loss/keV << G4endl; 926 } 928 } 927 else 929 else 928 { 930 { 929 // G4cout << "MGP old Tm = " << Tm << " 931 // G4cout << "MGP old Tm = " << Tm << " " << ipotFluct << " " << e0 << G4endl; 930 Tm = Tm-ipotFluct+e0 ; 932 Tm = Tm-ipotFluct+e0 ; 931 933 932 // MGP ---- workaround to avoid log argume 934 // MGP ---- workaround to avoid log argument<0, TO BE CHECKED 933 if (Tm <= 0.) 935 if (Tm <= 0.) 934 { 936 { 935 loss = MeanLoss; 937 loss = MeanLoss; 936 p3 = 0; 938 p3 = 0; 937 // G4cout << "MGP correction loss = Me 939 // G4cout << "MGP correction loss = MeanLoss " << loss/keV << G4endl; 938 } 940 } 939 else 941 else 940 { 942 { 941 a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log( 943 a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(Tm/e0)); 942 944 943 // G4cout << "MGP new Tm = " << Tm << 945 // G4cout << "MGP new Tm = " << Tm << " " << ipotFluct << " " << e0 << " a3= " << a3 << G4endl; 944 946 945 if(a3>alim) 947 if(a3>alim) 946 { 948 { 947 siga=std::sqrt(a3) ; 949 siga=std::sqrt(a3) ; 948 p3 = std::max(0,G4int(G4RandGauss::shoot 950 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5)); 949 } 951 } 950 else 952 else 951 p3 = G4Poisson(a3); 953 p3 = G4Poisson(a3); 952 //G4cout << "MGP p3 " << p3 << G4endl; 954 //G4cout << "MGP p3 " << p3 << G4endl; 953 955 954 } 956 } 955 957 956 if(p3 > 0) 958 if(p3 > 0) 957 { 959 { 958 w = (Tm-e0)/Tm ; 960 w = (Tm-e0)/Tm ; 959 if(p3 > nmaxCont2) 961 if(p3 > nmaxCont2) 960 { 962 { 961 // G4cout << "MGP dp3 " << dp3 << " p3 " 963 // G4cout << "MGP dp3 " << dp3 << " p3 " << p3 << " " << nmaxCont2 << G4endl; 962 dp3 = G4double(p3) ; 964 dp3 = G4double(p3) ; 963 Corrfac = dp3/G4double(nmaxCont2) ; 965 Corrfac = dp3/G4double(nmaxCont2) ; 964 p3 = nmaxCont2 ; 966 p3 = nmaxCont2 ; 965 } 967 } 966 else 968 else 967 Corrfac = 1. ; 969 Corrfac = 1. ; 968 970 969 for(G4int i=0; i<p3; i++) loss += 1./( 971 for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ; 970 loss *= e0*Corrfac ; 972 loss *= e0*Corrfac ; 971 // G4cout << "MGP Corrfac = " << Corrf 973 // G4cout << "MGP Corrfac = " << Corrfac << " e0 = " << e0/keV << " loss = " << loss/keV << G4endl; 972 } 974 } 973 } 975 } 974 } 976 } 975 977 976 else // not so 978 else // not so small Step 977 { 979 { 978 // excitation type 1 980 // excitation type 1 979 if(a1>alim) 981 if(a1>alim) 980 { 982 { 981 siga=std::sqrt(a1) ; 983 siga=std::sqrt(a1) ; 982 p1 = std::max(0,int(G4RandGauss::shoot 984 p1 = std::max(0,int(G4RandGauss::shoot(a1,siga)+0.5)); 983 } 985 } 984 else 986 else 985 p1 = G4Poisson(a1); 987 p1 = G4Poisson(a1); 986 988 987 // excitation type 2 989 // excitation type 2 988 if(a2>alim) 990 if(a2>alim) 989 { 991 { 990 siga=std::sqrt(a2) ; 992 siga=std::sqrt(a2) ; 991 p2 = std::max(0,int(G4RandGauss::shoot 993 p2 = std::max(0,int(G4RandGauss::shoot(a2,siga)+0.5)); 992 } 994 } 993 else 995 else 994 p2 = G4Poisson(a2); 996 p2 = G4Poisson(a2); 995 997 996 loss = p1*e1Fluct+p2*e2Fluct; 998 loss = p1*e1Fluct+p2*e2Fluct; 997 999 998 // smearing to avoid unphysical peaks 1000 // smearing to avoid unphysical peaks 999 if(p2 > 0) 1001 if(p2 > 0) 1000 loss += (1.-2.*G4UniformRand())*e2Flu 1002 loss += (1.-2.*G4UniformRand())*e2Fluct; 1001 else if (loss>0.) 1003 else if (loss>0.) 1002 loss += (1.-2.*G4UniformRand())*e1Flu 1004 loss += (1.-2.*G4UniformRand())*e1Fluct; 1003 1005 1004 // ionisation ......................... 1006 // ionisation ....................................... 1005 if(a3 > 0.) 1007 if(a3 > 0.) 1006 { 1008 { 1007 if(a3>alim) 1009 if(a3>alim) 1008 { 1010 { 1009 siga=std::sqrt(a3) ; 1011 siga=std::sqrt(a3) ; 1010 p3 = std::max(0,int(G4RandGauss::shoo 1012 p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5)); 1011 } 1013 } 1012 else 1014 else 1013 p3 = G4Poisson(a3); 1015 p3 = G4Poisson(a3); 1014 1016 1015 lossc = 0.; 1017 lossc = 0.; 1016 if(p3 > 0) 1018 if(p3 > 0) 1017 { 1019 { 1018 na = 0.; 1020 na = 0.; 1019 alfa = 1.; 1021 alfa = 1.; 1020 if (p3 > nmaxCont2) 1022 if (p3 > nmaxCont2) 1021 { 1023 { 1022 dp3 = G4double(p3); 1024 dp3 = G4double(p3); 1023 rfac = dp3/(G4double(nmaxCont2)+d 1025 rfac = dp3/(G4double(nmaxCont2)+dp3); 1024 namean = G4double(p3)*rfac; 1026 namean = G4double(p3)*rfac; 1025 sa = G4double(nmaxCont1)*rfac; 1027 sa = G4double(nmaxCont1)*rfac; 1026 na = G4RandGauss::shoot(namean, 1028 na = G4RandGauss::shoot(namean,sa); 1027 if (na > 0.) 1029 if (na > 0.) 1028 { 1030 { 1029 alfa = w1*G4double(nmaxCont2+p3)/ 1031 alfa = w1*G4double(nmaxCont2+p3)/(w1*G4double(nmaxCont2)+G4double(p3)); 1030 alfa1 = alfa*std::log(alfa)/(alfa- 1032 alfa1 = alfa*std::log(alfa)/(alfa-1.); 1031 ea = na*ipotFluct*alfa1; 1033 ea = na*ipotFluct*alfa1; 1032 sea = ipotFluct*std::sqrt(na*(al 1034 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1)); 1033 lossc += G4RandGauss::shoot(ea,sea) 1035 lossc += G4RandGauss::shoot(ea,sea); 1034 } 1036 } 1035 } 1037 } 1036 1038 1037 nb = G4int(G4double(p3)-na); 1039 nb = G4int(G4double(p3)-na); 1038 if (nb > 0) 1040 if (nb > 0) 1039 { 1041 { 1040 w2 = alfa*ipotFluct; 1042 w2 = alfa*ipotFluct; 1041 w = (Tm-w2)/Tm; 1043 w = (Tm-w2)/Tm; 1042 for (G4int k=0; k<nb; k++) lossc += w2/ 1044 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand()); 1043 } 1045 } 1044 } 1046 } 1045 1047 1046 loss += lossc; 1048 loss += lossc; 1047 } 1049 } 1048 } 1050 } 1049 1051 1050 return loss ; 1052 return loss ; 1051 } 1053 } 1052 1054 1053 // 1055 // 1054 1056 1055 1057