00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 #include "G4VEnergyLoss.hh"
00054 #include "G4PhysicalConstants.hh"
00055 #include "G4SystemOfUnits.hh"
00056 #include "G4EnergyLossMessenger.hh"
00057 #include "G4ProductionCutsTable.hh"
00058
00059
00060
00061 G4bool G4VEnergyLoss::rndmStepFlag = false;
00062 G4bool G4VEnergyLoss::EnlossFlucFlag = true;
00063
00064 G4bool G4VEnergyLoss::subSecFlag = false;
00065 G4double G4VEnergyLoss::MinDeltaCutInRange = 0.100*mm;
00066 G4double* G4VEnergyLoss::MinDeltaEnergy = 0;
00067 G4bool* G4VEnergyLoss::LowerLimitForced = 0;
00068 G4bool G4VEnergyLoss::setMinDeltaCutInRange = false;
00069
00070 G4double G4VEnergyLoss::dRoverRange = 20*perCent;
00071 G4double G4VEnergyLoss::finalRange = 1*mm;
00072 G4double G4VEnergyLoss::finalRangeRequested = -1*mm;
00073 G4double G4VEnergyLoss::c1lim = dRoverRange;
00074 G4double G4VEnergyLoss::c2lim = 2.*(1.-dRoverRange)*finalRange;
00075 G4double G4VEnergyLoss::c3lim = -(1.-dRoverRange)*finalRange*finalRange;
00076
00077 G4EnergyLossMessenger* G4VEnergyLoss::ELossMessenger = 0;
00078
00079
00080
00081 G4VEnergyLoss::G4VEnergyLoss(const G4String& aName , G4ProcessType aType)
00082 : G4VContinuousDiscreteProcess(aName, aType),
00083 lastMaterial(NULL),
00084 nmaxCont1(4),
00085 nmaxCont2(16)
00086 {
00087 if(!ELossMessenger) {
00088 G4cout << "### G4VEnergyLoss class is obsolete "
00089 << "and will be removed for the next release." << G4endl;
00090 ELossMessenger = new G4EnergyLossMessenger();
00091 }
00092
00093 imat = 0;
00094 f1Fluct = f2Fluct = e1Fluct = e2Fluct = rateFluct = ipotFluct = e1LogFluct
00095 = e2LogFluct = ipotLogFluct = taulow = tauhigh = ltaulow = ltauhigh
00096 = ParticleMass = 0.0;
00097 }
00098
00099
00100
00101 G4VEnergyLoss::~G4VEnergyLoss()
00102 {}
00103
00104
00105
00106 void G4VEnergyLoss::SetRndmStep(G4bool value) {rndmStepFlag = value;}
00107
00108
00109
00110 void G4VEnergyLoss::SetEnlossFluc(G4bool value) {EnlossFlucFlag = value;}
00111
00112
00113
00114 void G4VEnergyLoss::SetSubSec(G4bool value) {subSecFlag = value;}
00115
00116
00117
00118 void G4VEnergyLoss::SetMinDeltaCutInRange(G4double value)
00119 {MinDeltaCutInRange = value; setMinDeltaCutInRange = true;}
00120
00121
00122
00123 void G4VEnergyLoss::SetStepFunction(G4double c1, G4double c2)
00124 {
00125 dRoverRange = c1; finalRangeRequested = c2;
00126 c1lim=dRoverRange;
00127 c2lim=2.*(1-dRoverRange)*finalRangeRequested;
00128 c3lim=-(1.-dRoverRange)*finalRangeRequested*finalRangeRequested;
00129 }
00130
00131
00132
00133 G4PhysicsTable* G4VEnergyLoss::BuildRangeTable(
00134 G4PhysicsTable* theDEDXTable,G4PhysicsTable* theRangeTable,
00135 G4double LowestKineticEnergy,G4double HighestKineticEnergy,G4int TotBin)
00136
00137 {
00138 size_t numOfMaterials = theDEDXTable->length();
00139
00140 if(theRangeTable)
00141 { theRangeTable->clearAndDestroy();
00142 delete theRangeTable; }
00143 theRangeTable = new G4PhysicsTable(numOfMaterials);
00144
00145
00146
00147 for (size_t J=0; J<numOfMaterials; J++)
00148 {
00149 G4PhysicsLogVector* aVector;
00150 aVector = new G4PhysicsLogVector(LowestKineticEnergy,
00151 HighestKineticEnergy,TotBin);
00152
00153 BuildRangeVector(theDEDXTable,LowestKineticEnergy,HighestKineticEnergy,
00154 TotBin,J,aVector);
00155 theRangeTable->insert(aVector);
00156 }
00157 return theRangeTable ;
00158 }
00159
00160
00161
00162 void G4VEnergyLoss::BuildRangeVector(G4PhysicsTable* theDEDXTable,
00163 G4double,G4double ,G4int TotBin,
00164 G4int materialIndex,G4PhysicsLogVector* rangeVector)
00165
00166 {
00167 G4int nbin=100,i;
00168 G4bool isOut;
00169
00170 const G4double small = 1.e-6 ;
00171 const G4double masslimit = 0.52*MeV ;
00172
00173 G4double tlim=2.*MeV,t1=0.1*MeV,t2=0.025*MeV ;
00174 G4double tlime=0.2*keV,factor=2.*electron_mass_c2 ;
00175 G4double loss1,loss2,ca,cb,cba ;
00176 G4double losslim,clim ;
00177 G4double taulim,rangelim,ltaulim,
00178 LowEdgeEnergy,tau,Value,tau1,sqtau1 ;
00179 G4double oldValue,tauold ;
00180
00181 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
00182
00183
00184 G4double lossmin = +1.e10 ;
00185 for (G4int i1=0; i1<TotBin; i1++)
00186 {
00187 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i1) ;
00188 Value = physicsVector->GetValue(LowEdgeEnergy,isOut);
00189 if((Value < lossmin)&&(Value > 0.))
00190 lossmin = Value ;
00191 }
00192 for (G4int i2=0; i2<TotBin; i2++)
00193 {
00194 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i2) ;
00195 Value = physicsVector->GetValue(LowEdgeEnergy,isOut);
00196 if(Value < lossmin)
00197 physicsVector->PutValue(i2,small*lossmin) ;
00198 }
00199
00200
00201
00202 if(ParticleMass > masslimit)
00203 {
00204 loss1 = physicsVector->GetValue(t1,isOut);
00205 loss2 = physicsVector->GetValue(t2,isOut);
00206 tau1 = t1/ParticleMass ;
00207 sqtau1 = std::sqrt(tau1) ;
00208 ca = (4.*loss2-loss1)/sqtau1 ;
00209 cb = (2.*loss1-4.*loss2)/tau1 ;
00210 cba = cb/ca ;
00211 taulim = tlim/ParticleMass ;
00212 ltaulim = std::log(taulim) ;
00213
00214
00215 i=-1;
00216 oldValue = 0. ;
00217
00218 do
00219 {
00220 i += 1 ;
00221 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i);
00222 tau = LowEdgeEnergy/ParticleMass;
00223 if ( tau <= tau1 )
00224 {
00225 Value = 2.*ParticleMass*std::log(1.+cba*std::sqrt(tau))/cb ;
00226 }
00227 else
00228 {
00229 Value = 2.*ParticleMass*std::log(1.+cba*sqtau1)/cb ;
00230 if(tau<=taulim)
00231 {
00232 taulow = tau1 ;
00233 tauhigh = tau ;
00234 Value += RangeIntLin(physicsVector,nbin);
00235 }
00236 else
00237 {
00238
00239 taulow = tau1 ;
00240 tauhigh = taulim ;
00241 Value += RangeIntLin(physicsVector,nbin) ;
00242 ltaulow = ltaulim ;
00243 ltauhigh = std::log(tau) ;
00244 Value += RangeIntLog(physicsVector,nbin);
00245 }
00246 }
00247
00248 rangeVector->PutValue(i,Value);
00249 oldValue = Value ;
00250 tauold = tau ;
00251 } while (tau<=taulim) ;
00252 }
00253 else
00254
00255 {
00256 losslim = physicsVector->GetValue(tlime,isOut) ;
00257
00258 taulim = tlime/electron_mass_c2;
00259 clim = losslim;
00260 ltaulim = std::log(taulim);
00261
00262
00263 i=-1;
00264 oldValue = 0.;
00265
00266 do
00267 {
00268 i += 1 ;
00269 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i);
00270 tau = LowEdgeEnergy/electron_mass_c2;
00271 if (tau <= taulim) Value = factor*std::sqrt(tau*taulim)/clim;
00272 else {
00273 rangelim = factor*taulim/losslim ;
00274 ltaulow = std::log(taulim);
00275 ltauhigh = std::log(tau);
00276 Value = rangelim+RangeIntLog(physicsVector,nbin);
00277 }
00278 rangeVector->PutValue(i,Value);
00279 oldValue = Value;
00280 tauold = tau;
00281
00282 } while (tau<=taulim);
00283
00284 }
00285
00286 i += 1 ;
00287 for (G4int j=i; j<TotBin; j++)
00288 {
00289 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(j);
00290 tau = LowEdgeEnergy/ParticleMass;
00291 ltaulow = std::log(tauold);
00292 ltauhigh = std::log(tau);
00293 Value = oldValue+RangeIntLog(physicsVector,nbin);
00294 rangeVector->PutValue(j,Value);
00295 oldValue = Value ;
00296
00297 tauold = tau ;
00298 }
00299 }
00300
00301
00302
00303 G4double G4VEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector,
00304 G4int nbin)
00305
00306 {
00307 G4double dtau,Value,taui,ti,lossi,ci;
00308 G4bool isOut;
00309 dtau = (tauhigh-taulow)/nbin;
00310 Value = 0.;
00311
00312 for (G4int i=0; i<=nbin; i++)
00313 {
00314 taui = taulow + dtau*i ;
00315 ti = ParticleMass*taui;
00316 lossi = physicsVector->GetValue(ti,isOut);
00317 if(i==0)
00318 ci=0.5;
00319 else
00320 {
00321 if(i<nbin)
00322 ci=1.;
00323 else
00324 ci=0.5;
00325 }
00326 Value += ci/lossi;
00327 }
00328 Value *= ParticleMass*dtau;
00329 return Value;
00330 }
00331
00332
00333
00334 G4double G4VEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector,
00335 G4int nbin)
00336
00337 {
00338 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
00339 G4bool isOut;
00340 ltt = ltauhigh-ltaulow;
00341 dltau = ltt/nbin;
00342 Value = 0.;
00343
00344 for (G4int i=0; i<=nbin; i++)
00345 {
00346 ui = ltaulow+dltau*i;
00347 taui = std::exp(ui);
00348 ti = ParticleMass*taui;
00349 lossi = physicsVector->GetValue(ti,isOut);
00350 if(i==0)
00351 ci=0.5;
00352 else
00353 {
00354 if(i<nbin)
00355 ci=1.;
00356 else
00357 ci=0.5;
00358 }
00359 Value += ci*taui/lossi;
00360 }
00361 Value *= ParticleMass*dltau;
00362 return Value;
00363 }
00364
00365
00366
00367 G4PhysicsTable* G4VEnergyLoss::BuildLabTimeTable(G4PhysicsTable* theDEDXTable,
00368 G4PhysicsTable* theLabTimeTable,
00369 G4double LowestKineticEnergy,
00370 G4double HighestKineticEnergy,G4int TotBin)
00371
00372 {
00373 size_t numOfMaterials = theDEDXTable->length();
00374
00375 if(theLabTimeTable)
00376 { theLabTimeTable->clearAndDestroy();
00377 delete theLabTimeTable; }
00378 theLabTimeTable = new G4PhysicsTable(numOfMaterials);
00379
00380
00381 for (size_t J=0; J<numOfMaterials; J++)
00382 {
00383 G4PhysicsLogVector* aVector;
00384
00385 aVector = new G4PhysicsLogVector(LowestKineticEnergy,
00386 HighestKineticEnergy,TotBin);
00387
00388 BuildLabTimeVector(theDEDXTable,
00389 LowestKineticEnergy,HighestKineticEnergy,TotBin,J,aVector);
00390 theLabTimeTable->insert(aVector);
00391
00392
00393 }
00394 return theLabTimeTable ;
00395 }
00396
00397
00398
00399 G4PhysicsTable* G4VEnergyLoss::BuildProperTimeTable(G4PhysicsTable* theDEDXTable,
00400 G4PhysicsTable* theProperTimeTable,
00401 G4double LowestKineticEnergy,
00402 G4double HighestKineticEnergy,G4int TotBin)
00403
00404 {
00405 size_t numOfMaterials = theDEDXTable->length();
00406
00407 if(theProperTimeTable)
00408 { theProperTimeTable->clearAndDestroy();
00409 delete theProperTimeTable; }
00410 theProperTimeTable = new G4PhysicsTable(numOfMaterials);
00411
00412
00413 for (size_t J=0; J<numOfMaterials; J++)
00414 {
00415 G4PhysicsLogVector* aVector;
00416
00417 aVector = new G4PhysicsLogVector(LowestKineticEnergy,
00418 HighestKineticEnergy,TotBin);
00419
00420 BuildProperTimeVector(theDEDXTable,
00421 LowestKineticEnergy,HighestKineticEnergy,TotBin,J,aVector);
00422 theProperTimeTable->insert(aVector);
00423
00424
00425 }
00426 return theProperTimeTable ;
00427 }
00428
00429
00430
00431 void G4VEnergyLoss::BuildLabTimeVector(G4PhysicsTable* theDEDXTable,
00432 G4double,
00433 G4double ,G4int TotBin,
00434 G4int materialIndex, G4PhysicsLogVector* timeVector)
00435
00436 {
00437
00438 G4int nbin=100;
00439 G4bool isOut;
00440 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
00441 G4double losslim,clim,taulim,timelim,
00442 LowEdgeEnergy,tau,Value ;
00443
00444 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
00445
00446
00447 losslim = physicsVector->GetValue(tlim,isOut);
00448 taulim=tlim/ParticleMass ;
00449 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
00450
00451
00452
00453 G4int i=-1;
00454 G4double oldValue = 0. ;
00455 G4double tauold ;
00456 do
00457 {
00458 i += 1 ;
00459 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
00460 tau = LowEdgeEnergy/ParticleMass ;
00461 if ( tau <= taulim )
00462 {
00463 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
00464 }
00465 else
00466 {
00467 timelim=clim ;
00468 ltaulow = std::log(taulim);
00469 ltauhigh = std::log(tau);
00470 Value = timelim+LabTimeIntLog(physicsVector,nbin);
00471 }
00472 timeVector->PutValue(i,Value);
00473 oldValue = Value ;
00474 tauold = tau ;
00475 } while (tau<=taulim) ;
00476 i += 1 ;
00477 for (G4int j=i; j<TotBin; j++)
00478 {
00479 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
00480 tau = LowEdgeEnergy/ParticleMass ;
00481 ltaulow = std::log(tauold);
00482 ltauhigh = std::log(tau);
00483 Value = oldValue+LabTimeIntLog(physicsVector,nbin);
00484 timeVector->PutValue(j,Value);
00485 oldValue = Value ;
00486 tauold = tau ;
00487 }
00488 }
00489
00490
00491
00492 void G4VEnergyLoss::BuildProperTimeVector(G4PhysicsTable* theDEDXTable,
00493 G4double,
00494 G4double ,G4int TotBin,
00495 G4int materialIndex, G4PhysicsLogVector* timeVector)
00496
00497 {
00498 G4int nbin=100;
00499 G4bool isOut;
00500 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
00501 G4double losslim,clim,taulim,timelim,
00502 LowEdgeEnergy,tau,Value ;
00503
00504 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
00505
00506
00507 losslim = physicsVector->GetValue(tlim,isOut);
00508 taulim=tlim/ParticleMass ;
00509 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
00510
00511
00512
00513 G4int i=-1;
00514 G4double oldValue = 0. ;
00515 G4double tauold ;
00516 do
00517 {
00518 i += 1 ;
00519 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
00520 tau = LowEdgeEnergy/ParticleMass ;
00521 if ( tau <= taulim )
00522 {
00523 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
00524 }
00525 else
00526 {
00527 timelim=clim ;
00528 ltaulow = std::log(taulim);
00529 ltauhigh = std::log(tau);
00530 Value = timelim+ProperTimeIntLog(physicsVector,nbin);
00531 }
00532 timeVector->PutValue(i,Value);
00533 oldValue = Value ;
00534 tauold = tau ;
00535 } while (tau<=taulim) ;
00536 i += 1 ;
00537 for (G4int j=i; j<TotBin; j++)
00538 {
00539 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
00540 tau = LowEdgeEnergy/ParticleMass ;
00541 ltaulow = std::log(tauold);
00542 ltauhigh = std::log(tau);
00543 Value = oldValue+ProperTimeIntLog(physicsVector,nbin);
00544 timeVector->PutValue(j,Value);
00545 oldValue = Value ;
00546 tauold = tau ;
00547 }
00548 }
00549
00550
00551
00552 G4double G4VEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector,
00553 G4int nbin)
00554
00555 {
00556 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
00557 G4bool isOut;
00558 ltt = ltauhigh-ltaulow;
00559 dltau = ltt/nbin;
00560 Value = 0.;
00561
00562 for (G4int i=0; i<=nbin; i++)
00563 {
00564 ui = ltaulow+dltau*i;
00565 taui = std::exp(ui);
00566 ti = ParticleMass*taui;
00567 lossi = physicsVector->GetValue(ti,isOut);
00568 if(i==0)
00569 ci=0.5;
00570 else
00571 {
00572 if(i<nbin)
00573 ci=1.;
00574 else
00575 ci=0.5;
00576 }
00577 Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
00578 }
00579 Value *= ParticleMass*dltau/c_light;
00580 return Value;
00581 }
00582
00583
00584
00585 G4double G4VEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector,
00586 G4int nbin)
00587
00588 {
00589 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
00590 G4bool isOut;
00591 ltt = ltauhigh-ltaulow;
00592 dltau = ltt/nbin;
00593 Value = 0.;
00594
00595 for (G4int i=0; i<=nbin; i++)
00596 {
00597 ui = ltaulow+dltau*i;
00598 taui = std::exp(ui);
00599 ti = ParticleMass*taui;
00600 lossi = physicsVector->GetValue(ti,isOut);
00601 if(i==0)
00602 ci=0.5;
00603 else
00604 {
00605 if(i<nbin)
00606 ci=1.;
00607 else
00608 ci=0.5;
00609 }
00610 Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
00611 }
00612 Value *= ParticleMass*dltau/c_light;
00613 return Value;
00614 }
00615
00616
00617
00618 G4PhysicsTable* G4VEnergyLoss::BuildInverseRangeTable(G4PhysicsTable* theRangeTable,
00619 G4PhysicsTable* theRangeCoeffATable,
00620 G4PhysicsTable* theRangeCoeffBTable,
00621 G4PhysicsTable* theRangeCoeffCTable,
00622 G4PhysicsTable* theInverseRangeTable,
00623 G4double LowestKineticEnergy,
00624 G4double HighestKineticEnergy,G4int TotBin)
00625
00626 {
00627 G4double SmallestRange,BiggestRange ;
00628 G4bool isOut ;
00629 size_t numOfMaterials = theRangeTable->length();
00630
00631 if(theInverseRangeTable)
00632 { theInverseRangeTable->clearAndDestroy();
00633 delete theInverseRangeTable; }
00634 theInverseRangeTable = new G4PhysicsTable(numOfMaterials);
00635
00636
00637 for (size_t J=0; J<numOfMaterials; J++)
00638 {
00639 SmallestRange = (*theRangeTable)(J)->
00640 GetValue(LowestKineticEnergy,isOut) ;
00641 BiggestRange = (*theRangeTable)(J)->
00642 GetValue(HighestKineticEnergy,isOut) ;
00643 G4PhysicsLogVector* aVector;
00644 aVector = new G4PhysicsLogVector(SmallestRange,
00645 BiggestRange,TotBin);
00646
00647 InvertRangeVector(theRangeTable,
00648 theRangeCoeffATable,
00649 theRangeCoeffBTable,
00650 theRangeCoeffCTable,
00651 LowestKineticEnergy,HighestKineticEnergy,TotBin,J, aVector);
00652
00653 theInverseRangeTable->insert(aVector);
00654 }
00655 return theInverseRangeTable ;
00656 }
00657
00658
00659
00660 void G4VEnergyLoss::InvertRangeVector(G4PhysicsTable* theRangeTable,
00661 G4PhysicsTable* theRangeCoeffATable,
00662 G4PhysicsTable* theRangeCoeffBTable,
00663 G4PhysicsTable* theRangeCoeffCTable,
00664 G4double LowestKineticEnergy,
00665 G4double HighestKineticEnergy,G4int TotBin,
00666 G4int materialIndex,G4PhysicsLogVector* aVector)
00667
00668 {
00669 G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
00670 G4double RTable = std::exp(std::log(HighestKineticEnergy/LowestKineticEnergy)/TotBin) ;
00671 G4double Tbin = LowestKineticEnergy/RTable ;
00672 G4double rangebin = 0.0 ;
00673 G4int binnumber = -1 ;
00674 G4bool isOut ;
00675
00676
00677 for( G4int i=0; i<TotBin; i++)
00678 {
00679 LowEdgeRange = aVector->GetLowEdgeEnergy(i) ;
00680 if( rangebin < LowEdgeRange )
00681 {
00682 do
00683 {
00684 binnumber += 1 ;
00685 Tbin *= RTable ;
00686 rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
00687 }
00688 while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
00689 }
00690
00691 if(binnumber == 0)
00692 KineticEnergy = LowestKineticEnergy ;
00693 else if(binnumber == TotBin-1)
00694 KineticEnergy = HighestKineticEnergy ;
00695 else
00696 {
00697 A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
00698 B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
00699 C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
00700 if(A==0.)
00701 KineticEnergy = (LowEdgeRange -C )/B ;
00702 else
00703 {
00704 discr = B*B - 4.*A*(C-LowEdgeRange);
00705 discr = discr>0. ? std::sqrt(discr) : 0.;
00706 KineticEnergy = 0.5*(discr-B)/A ;
00707 }
00708 }
00709
00710 aVector->PutValue(i,KineticEnergy) ;
00711 }
00712 }
00713
00714
00715
00716 G4PhysicsTable* G4VEnergyLoss::BuildRangeCoeffATable(G4PhysicsTable* theRangeTable,
00717 G4PhysicsTable* theRangeCoeffATable,
00718 G4double LowestKineticEnergy,
00719 G4double HighestKineticEnergy,G4int TotBin)
00720
00721
00722 {
00723 G4int numOfMaterials = theRangeTable->length();
00724
00725 if(theRangeCoeffATable)
00726 { theRangeCoeffATable->clearAndDestroy();
00727 delete theRangeCoeffATable; }
00728 theRangeCoeffATable = new G4PhysicsTable(numOfMaterials);
00729
00730 G4double RTable = std::exp(std::log(HighestKineticEnergy/LowestKineticEnergy)/TotBin) ;
00731 G4double R2 = RTable*RTable ;
00732 G4double R1 = RTable+1.;
00733 G4double w = R1*(RTable-1.)*(RTable-1.);
00734 G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ;
00735 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
00736 G4bool isOut;
00737
00738
00739 for (G4int J=0; J<numOfMaterials; J++)
00740 {
00741 G4int binmax=TotBin ;
00742 G4PhysicsLinearVector* aVector =
00743 new G4PhysicsLinearVector(0.,binmax, TotBin);
00744 Ti = LowestKineticEnergy ;
00745 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
00746
00747 for ( G4int i=0; i<TotBin; i++)
00748 {
00749 Ri = rangeVector->GetValue(Ti,isOut) ;
00750 if ( i==0 )
00751 Rim = 0. ;
00752 else
00753 {
00754 Tim = Ti/RTable ;
00755 Rim = rangeVector->GetValue(Tim,isOut);
00756 }
00757 if ( i==(TotBin-1))
00758 Rip = Ri ;
00759 else
00760 {
00761 Tip = Ti*RTable ;
00762 Rip = rangeVector->GetValue(Tip,isOut);
00763 }
00764 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ;
00765
00766 aVector->PutValue(i,Value);
00767 Ti = RTable*Ti ;
00768 }
00769
00770 theRangeCoeffATable->insert(aVector);
00771 }
00772 return theRangeCoeffATable ;
00773 }
00774
00775
00776
00777 G4PhysicsTable* G4VEnergyLoss::BuildRangeCoeffBTable(G4PhysicsTable* theRangeTable,
00778 G4PhysicsTable* theRangeCoeffBTable,
00779 G4double LowestKineticEnergy,
00780 G4double HighestKineticEnergy,G4int TotBin)
00781
00782
00783 {
00784 G4int numOfMaterials = theRangeTable->length();
00785
00786 if(theRangeCoeffBTable)
00787 { theRangeCoeffBTable->clearAndDestroy();
00788 delete theRangeCoeffBTable; }
00789 theRangeCoeffBTable = new G4PhysicsTable(numOfMaterials);
00790
00791 G4double RTable = std::exp(std::log(HighestKineticEnergy/LowestKineticEnergy)/TotBin) ;
00792 G4double R2 = RTable*RTable ;
00793 G4double R1 = RTable+1.;
00794 G4double w = R1*(RTable-1.)*(RTable-1.);
00795 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 = -R2*R1/w ;
00796 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
00797 G4bool isOut;
00798
00799
00800 for (G4int J=0; J<numOfMaterials; J++)
00801 {
00802 G4int binmax=TotBin ;
00803 G4PhysicsLinearVector* aVector =
00804 new G4PhysicsLinearVector(0.,binmax, TotBin);
00805 Ti = LowestKineticEnergy ;
00806 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
00807
00808 for ( G4int i=0; i<TotBin; i++)
00809 {
00810 Ri = rangeVector->GetValue(Ti,isOut) ;
00811 if ( i==0 )
00812 Rim = 0. ;
00813 else
00814 {
00815 Tim = Ti/RTable ;
00816 Rim = rangeVector->GetValue(Tim,isOut);
00817 }
00818 if ( i==(TotBin-1))
00819 Rip = Ri ;
00820 else
00821 {
00822 Tip = Ti*RTable ;
00823 Rip = rangeVector->GetValue(Tip,isOut);
00824 }
00825 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
00826
00827 aVector->PutValue(i,Value);
00828 Ti = RTable*Ti ;
00829 }
00830 theRangeCoeffBTable->insert(aVector);
00831 }
00832 return theRangeCoeffBTable ;
00833 }
00834
00835
00836
00837 G4PhysicsTable* G4VEnergyLoss::BuildRangeCoeffCTable(G4PhysicsTable* theRangeTable,
00838 G4PhysicsTable* theRangeCoeffCTable,
00839 G4double LowestKineticEnergy,
00840 G4double HighestKineticEnergy,G4int TotBin)
00841
00842
00843 {
00844 G4int numOfMaterials = theRangeTable->length();
00845
00846 if(theRangeCoeffCTable)
00847 { theRangeCoeffCTable->clearAndDestroy();
00848 delete theRangeCoeffCTable; }
00849 theRangeCoeffCTable = new G4PhysicsTable(numOfMaterials);
00850
00851 G4double RTable = std::exp(std::log(HighestKineticEnergy/LowestKineticEnergy)/TotBin) ;
00852 G4double R2 = RTable*RTable ;
00853 G4double R1 = RTable+1.;
00854 G4double w = R1*(RTable-1.)*(RTable-1.);
00855 G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ;
00856 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
00857 G4bool isOut;
00858
00859
00860 for (G4int J=0; J<numOfMaterials; J++)
00861 {
00862 G4int binmax=TotBin ;
00863 G4PhysicsLinearVector* aVector =
00864 new G4PhysicsLinearVector(0.,binmax, TotBin);
00865 Ti = LowestKineticEnergy ;
00866 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
00867
00868 for ( G4int i=0; i<TotBin; i++)
00869 {
00870 Ri = rangeVector->GetValue(Ti,isOut) ;
00871 if ( i==0 )
00872 Rim = 0. ;
00873 else
00874 {
00875 Tim = Ti/RTable ;
00876 Rim = rangeVector->GetValue(Tim,isOut);
00877 }
00878 if ( i==(TotBin-1))
00879 Rip = Ri ;
00880 else
00881 {
00882 Tip = Ti*RTable ;
00883 Rip = rangeVector->GetValue(Tip,isOut);
00884 }
00885 Value = w1*Rip + w2*Ri + w3*Rim ;
00886
00887 aVector->PutValue(i,Value);
00888 Ti = RTable*Ti ;
00889 }
00890 theRangeCoeffCTable->insert(aVector);
00891 }
00892 return theRangeCoeffCTable ;
00893 }
00894
00895
00896
00897 G4double G4VEnergyLoss::GetLossWithFluct(const G4DynamicParticle* aParticle,
00898 const G4MaterialCutsCouple* couple,
00899 G4double ChargeSquare,
00900 G4double MeanLoss,
00901 G4double step )
00902
00903
00904 {
00905 const G4double minLoss = 1.*eV ;
00906 const G4double probLim = 0.01 ;
00907 const G4double sumaLim = -std::log(probLim) ;
00908 const G4double alim=10.;
00909 const G4double kappa = 10. ;
00910 const G4double factor = twopi_mc2_rcl2 ;
00911 const G4Material* aMaterial = couple->GetMaterial();
00912
00913
00914
00915 if (aMaterial != lastMaterial)
00916 {
00917 lastMaterial = aMaterial;
00918 imat = couple->GetIndex();
00919 f1Fluct = aMaterial->GetIonisation()->GetF1fluct();
00920 f2Fluct = aMaterial->GetIonisation()->GetF2fluct();
00921 e1Fluct = aMaterial->GetIonisation()->GetEnergy1fluct();
00922 e2Fluct = aMaterial->GetIonisation()->GetEnergy2fluct();
00923 e1LogFluct = aMaterial->GetIonisation()->GetLogEnergy1fluct();
00924 e2LogFluct = aMaterial->GetIonisation()->GetLogEnergy2fluct();
00925 rateFluct = aMaterial->GetIonisation()->GetRateionexcfluct();
00926 ipotFluct = aMaterial->GetIonisation()->GetMeanExcitationEnergy();
00927 ipotLogFluct = aMaterial->GetIonisation()->GetLogMeanExcEnergy();
00928 }
00929 G4double threshold,w1,w2,C,
00930 beta2,suma,e0,loss,lossc ,w,electronDensity;
00931 G4double a1,a2,a3;
00932 G4double p1,p2,p3 ;
00933 G4int nb;
00934 G4double Corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
00935 G4double dp3;
00936 G4double siga ;
00937
00938
00939 if(MeanLoss < minLoss) return MeanLoss ;
00940
00941
00942 G4double Tkin = aParticle->GetKineticEnergy();
00943 ParticleMass = aParticle->GetMass() ;
00944
00945 threshold = (*((G4ProductionCutsTable::GetProductionCutsTable())
00946 ->GetEnergyCutsVector(1)))[imat];
00947 G4double rmass = electron_mass_c2/ParticleMass;
00948 G4double tau = Tkin/ParticleMass, tau1 = tau+1., tau2 = tau*(tau+2.);
00949 G4double Tm = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
00950
00951 if(Tm > threshold) Tm = threshold;
00952
00953 beta2 = tau2/(tau1*tau1);
00954
00955
00956 if(MeanLoss >= kappa*Tm || Tm < kappa*ipotFluct)
00957 {
00958 electronDensity = aMaterial->GetElectronDensity() ;
00959 siga = std::sqrt(Tm*(1.0-0.5*beta2)*step*
00960 factor*electronDensity*ChargeSquare/beta2) ;
00961 do {
00962 loss = G4RandGauss::shoot(MeanLoss,siga) ;
00963 } while (loss < 0. || loss > 2.*MeanLoss);
00964 return loss;
00965 }
00966
00967 w1 = Tm/ipotFluct;
00968 w2 = std::log(2.*electron_mass_c2*tau2);
00969
00970 C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
00971
00972 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
00973 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
00974 a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipotFluct*Tm*std::log(w1));
00975 if(a1 < 0.) a1 = 0.;
00976 if(a2 < 0.) a2 = 0.;
00977 if(a3 < 0.) a3 = 0.;
00978
00979 suma = a1+a2+a3;
00980
00981 loss = 0. ;
00982
00983 if(suma < sumaLim)
00984 {
00985 e0 = aMaterial->GetIonisation()->GetEnergy0fluct();
00986
00987 if(Tm == ipotFluct)
00988 {
00989 a3 = MeanLoss/e0;
00990
00991 if(a3>alim)
00992 {
00993 siga=std::sqrt(a3) ;
00994 p3 = std::max(0.,G4RandGauss::shoot(a3,siga)+0.5);
00995 }
00996 else
00997 p3 = G4float(G4Poisson(a3));
00998
00999 loss = p3*e0 ;
01000
01001 if(p3 > 0)
01002 loss += (1.-2.*G4UniformRand())*e0 ;
01003
01004 }
01005 else
01006 {
01007 Tm = Tm-ipotFluct+e0 ;
01008 a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(Tm/e0));
01009
01010 if(a3>alim)
01011 {
01012 siga=std::sqrt(a3) ;
01013 p3 = std::max(0.,G4RandGauss::shoot(a3,siga)+0.5);
01014 }
01015 else
01016 p3 = G4float(G4Poisson(a3));
01017
01018 if(p3 > 0)
01019 {
01020 w = (Tm-e0)/Tm ;
01021 if(p3 > G4float(nmaxCont2))
01022 {
01023 dp3 = p3 ;
01024 Corrfac = dp3/G4float(nmaxCont2) ;
01025 p3 = G4float(nmaxCont2) ;
01026 }
01027 else
01028 Corrfac = 1. ;
01029
01030 for(G4int i=0; i<G4int(p3); i++) loss += 1./(1.-w*G4UniformRand()) ;
01031 loss *= e0*Corrfac ;
01032 }
01033 }
01034 }
01035
01036 else
01037 {
01038
01039 if(a1>alim)
01040 {
01041 siga=std::sqrt(a1) ;
01042 p1 = std::max(0.,G4RandGauss::shoot(a1,siga)+0.5);
01043 }
01044 else
01045 p1 = G4float(G4Poisson(a1));
01046
01047
01048 if(a2>alim)
01049 {
01050 siga=std::sqrt(a2) ;
01051 p2 = std::max(0.,G4RandGauss::shoot(a2,siga)+0.5);
01052 }
01053 else
01054 p2 = G4float(G4Poisson(a2));
01055
01056 loss = p1*e1Fluct+p2*e2Fluct;
01057
01058 if(p2 > 0)
01059 loss += (1.-2.*G4UniformRand())*e2Fluct;
01060 else if (loss>0.)
01061 loss += (1.-2.*G4UniformRand())*e1Fluct;
01062
01063
01064 if(a3 > 0.)
01065 {
01066 if(a3>alim)
01067 {
01068 siga=std::sqrt(a3) ;
01069 p3 = std::max(0.,G4RandGauss::shoot(a3,siga)+0.5);
01070 }
01071 else
01072 p3 = G4float(G4Poisson(a3));
01073
01074 lossc = 0.;
01075 if(p3 > 0)
01076 {
01077 na = 0.;
01078 alfa = 1.;
01079 if (p3 > G4float(nmaxCont2))
01080 {
01081 dp3 = p3;
01082 rfac = dp3/(G4float(nmaxCont2)+dp3);
01083 namean = p3*rfac;
01084 sa = G4float(nmaxCont1)*rfac;
01085 na = G4RandGauss::shoot(namean,sa);
01086 if (na > 0.)
01087 {
01088 alfa = w1*(G4float(nmaxCont2)+p3)/(w1*G4float(nmaxCont2)+p3);
01089 alfa1 = alfa*std::log(alfa)/(alfa-1.);
01090 ea = na*ipotFluct*alfa1;
01091 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
01092 lossc += G4RandGauss::shoot(ea,sea);
01093 }
01094 }
01095
01096 nb = G4int(p3-na);
01097 if (nb > 0)
01098 {
01099 w2 = alfa*ipotFluct;
01100 w = (Tm-w2)/Tm;
01101 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
01102 }
01103 }
01104 loss += lossc;
01105 }
01106 }
01107
01108 if( loss < 0.)
01109 loss = 0.;
01110
01111 return loss ;
01112 }
01113
01114
01115
01116 G4bool G4VEnergyLoss::EqualCutVectors( G4double* vec1, G4double* vec2 )
01117 {
01118 if ( (vec1==0 ) || (vec2==0) ) return false;
01119
01120 G4bool flag = true;
01121
01122 for (size_t j=0; flag && j<G4Material::GetNumberOfMaterials(); j++){
01123 flag = (vec1[j] == vec2[j]);
01124 }
01125
01126 return flag;
01127 }
01128
01129
01130
01131 G4double* G4VEnergyLoss::CopyCutVectors( G4double* dest, G4double* source )
01132 {
01133 if ( dest != 0) delete [] dest;
01134 dest = new G4double [G4Material::GetNumberOfMaterials()];
01135 for (size_t j=0; j<G4Material::GetNumberOfMaterials(); j++){
01136 dest[j] = source[j];
01137 }
01138 return dest;
01139 }
01140
01141
01142
01143 G4bool G4VEnergyLoss::CutsWhereModified()
01144 {
01145 G4bool wasModified = false;
01146 const G4ProductionCutsTable* theCoupleTable=
01147 G4ProductionCutsTable::GetProductionCutsTable();
01148 size_t numOfCouples = theCoupleTable->GetTableSize();
01149
01150 for (size_t j=0; j<numOfCouples; j++){
01151 if (theCoupleTable->GetMaterialCutsCouple(j)->IsRecalcNeeded()) {
01152 wasModified = true;
01153 break;
01154 }
01155 }
01156 return wasModified;
01157 }
01158
01159