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