G4VeLowEnergyLoss Class Reference

#include <G4VeLowEnergyLoss.hh>

Inheritance diagram for G4VeLowEnergyLoss:

G4VContinuousDiscreteProcess G4VProcess

Public Member Functions

 G4VeLowEnergyLoss (const G4String &, G4ProcessType aType=fNotDefined)
 G4VeLowEnergyLoss (G4VeLowEnergyLoss &)
virtual ~G4VeLowEnergyLoss ()
virtual G4double GetContinuousStepLimit (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety)=0
virtual G4VParticleChangeAlongStepDoIt (const G4Track &track, const G4Step &Step)=0
virtual G4double GetMeanFreePath (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)=0
virtual G4VParticleChangePostStepDoIt (const G4Track &track, const G4Step &Step)=0

Static Public Member Functions

static void SetRndmStep (G4bool value)
static void SetEnlossFluc (G4bool value)
static void SetStepFunction (G4double c1, G4double c2)

Protected Member Functions

G4double GetLossWithFluct (const G4DynamicParticle *aParticle, const G4MaterialCutsCouple *couple, G4double MeanLoss, G4double step)

Static Protected Member Functions

static G4PhysicsTableBuildRangeTable (G4PhysicsTable *theDEDXTable, G4PhysicsTable *theRangeTable, G4double Tmin, G4double Tmax, G4int nbin)
static G4PhysicsTableBuildLabTimeTable (G4PhysicsTable *theDEDXTable, G4PhysicsTable *theLabTimeTable, G4double Tmin, G4double Tmax, G4int nbin)
static G4PhysicsTableBuildProperTimeTable (G4PhysicsTable *theDEDXTable, G4PhysicsTable *ProperTimeTable, G4double Tmin, G4double Tmax, G4int nbin)
static G4PhysicsTableBuildRangeCoeffATable (G4PhysicsTable *theRangeTable, G4PhysicsTable *theCoeffATable, G4double Tmin, G4double Tmax, G4int nbin)
static G4PhysicsTableBuildRangeCoeffBTable (G4PhysicsTable *theRangeTable, G4PhysicsTable *theCoeffBTable, G4double Tmin, G4double Tmax, G4int nbin)
static G4PhysicsTableBuildRangeCoeffCTable (G4PhysicsTable *theRangeTable, G4PhysicsTable *theCoeffCTable, G4double Tmin, G4double Tmax, G4int nbin)
static G4PhysicsTableBuildInverseRangeTable (G4PhysicsTable *theRangeTable, G4PhysicsTable *theRangeCoeffATable, G4PhysicsTable *theRangeCoeffBTable, G4PhysicsTable *theRangeCoeffCTable, G4PhysicsTable *theInverseRangeTable, G4double Tmin, G4double Tmax, G4int nbin)

Protected Attributes

const G4MateriallastMaterial
G4int imat
G4double f1Fluct
G4double f2Fluct
G4double e1Fluct
G4double e2Fluct
G4double rateFluct
G4double ipotFluct
G4double e1LogFluct
G4double e2LogFluct
G4double ipotLogFluct
const G4int nmaxCont1
const G4int nmaxCont2

Static Protected Attributes

static G4double ParticleMass
static G4double taulow
static G4double tauhigh
static G4double ltaulow
static G4double ltauhigh
static G4double dRoverRange = 20*perCent
static G4double finalRange = 200*micrometer
static G4double c1lim = dRoverRange
static G4double c2lim = 2.*(1.-dRoverRange)*finalRange
static G4double c3lim = -(1.-dRoverRange)*finalRange*finalRange
static G4bool rndmStepFlag = false
static G4bool EnlossFlucFlag = true

Detailed Description

Definition at line 77 of file G4VeLowEnergyLoss.hh.


Constructor & Destructor Documentation

G4VeLowEnergyLoss::G4VeLowEnergyLoss ( const G4String ,
G4ProcessType  aType = fNotDefined 
)

Definition at line 87 of file G4VeLowEnergyLoss.cc.

References G4VeLowEnergyLoss().

Referenced by G4VeLowEnergyLoss().

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 {;}

G4VeLowEnergyLoss::G4VeLowEnergyLoss ( G4VeLowEnergyLoss  ) 

Definition at line 105 of file G4VeLowEnergyLoss.cc.

References G4VeLowEnergyLoss().

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 {;}

G4VeLowEnergyLoss::~G4VeLowEnergyLoss (  )  [virtual]

Definition at line 99 of file G4VeLowEnergyLoss.cc.

00100 {
00101 }


Member Function Documentation

virtual G4VParticleChange* G4VeLowEnergyLoss::AlongStepDoIt ( const G4Track track,
const G4Step Step 
) [pure virtual]

Reimplemented from G4VContinuousDiscreteProcess.

G4PhysicsTable * G4VeLowEnergyLoss::BuildInverseRangeTable ( G4PhysicsTable theRangeTable,
G4PhysicsTable theRangeCoeffATable,
G4PhysicsTable theRangeCoeffBTable,
G4PhysicsTable theRangeCoeffCTable,
G4PhysicsTable theInverseRangeTable,
G4double  Tmin,
G4double  Tmax,
G4int  nbin 
) [static, protected]

Definition at line 524 of file G4VeLowEnergyLoss.cc.

References G4PhysicsTable::clearAndDestroy(), G4PhysicsVector::GetLowEdgeEnergy(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4PhysicsTable::insert(), and G4PhysicsVector::PutValue().

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   // loop for materials
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     //rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1)));
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 }

G4PhysicsTable * G4VeLowEnergyLoss::BuildLabTimeTable ( G4PhysicsTable theDEDXTable,
G4PhysicsTable theLabTimeTable,
G4double  Tmin,
G4double  Tmax,
G4int  nbin 
) [static, protected]

Definition at line 268 of file G4VeLowEnergyLoss.cc.

References G4PhysicsTable::clearAndDestroy(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), and G4PhysicsTable::insert().

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 }

G4PhysicsTable * G4VeLowEnergyLoss::BuildProperTimeTable ( G4PhysicsTable theDEDXTable,
G4PhysicsTable ProperTimeTable,
G4double  Tmin,
G4double  Tmax,
G4int  nbin 
) [static, protected]

Definition at line 301 of file G4VeLowEnergyLoss.cc.

References G4PhysicsTable::clearAndDestroy(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), and G4PhysicsTable::insert().

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 }

G4PhysicsTable * G4VeLowEnergyLoss::BuildRangeCoeffATable ( G4PhysicsTable theRangeTable,
G4PhysicsTable theCoeffATable,
G4double  Tmin,
G4double  Tmax,
G4int  nbin 
) [static, protected]

Definition at line 649 of file G4VeLowEnergyLoss.cc.

References G4PhysicsTable::clearAndDestroy(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4PhysicsVector::GetValue(), G4PhysicsTable::insert(), and G4PhysicsVector::PutValue().

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   //  loop for materials
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 }

G4PhysicsTable * G4VeLowEnergyLoss::BuildRangeCoeffBTable ( G4PhysicsTable theRangeTable,
G4PhysicsTable theCoeffBTable,
G4double  Tmin,
G4double  Tmax,
G4int  nbin 
) [static, protected]

Definition at line 711 of file G4VeLowEnergyLoss.cc.

References G4PhysicsTable::clearAndDestroy(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4PhysicsVector::GetValue(), G4PhysicsTable::insert(), and G4PhysicsVector::PutValue().

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   //  loop for materials
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 }

G4PhysicsTable * G4VeLowEnergyLoss::BuildRangeCoeffCTable ( G4PhysicsTable theRangeTable,
G4PhysicsTable theCoeffCTable,
G4double  Tmin,
G4double  Tmax,
G4int  nbin 
) [static, protected]

Definition at line 772 of file G4VeLowEnergyLoss.cc.

References G4PhysicsTable::clearAndDestroy(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), G4PhysicsVector::GetValue(), G4PhysicsTable::insert(), and G4PhysicsVector::PutValue().

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   //  loop for materials
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 }

G4PhysicsTable * G4VeLowEnergyLoss::BuildRangeTable ( G4PhysicsTable theDEDXTable,
G4PhysicsTable theRangeTable,
G4double  Tmin,
G4double  Tmax,
G4int  nbin 
) [static, protected]

Definition at line 134 of file G4VeLowEnergyLoss.cc.

References G4PhysicsTable::clearAndDestroy(), G4PhysicsTable::insert(), and G4PhysicsTable::length().

00140 {
00141 
00142    G4int numOfCouples = theDEDXTable->length();
00143 
00144    if(theRangeTable)
00145    { theRangeTable->clearAndDestroy();
00146      delete theRangeTable; }
00147    theRangeTable = new G4PhysicsTable(numOfCouples);
00148 
00149    // loop for materials
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 }

virtual G4double G4VeLowEnergyLoss::GetContinuousStepLimit ( const G4Track track,
G4double  previousStepSize,
G4double  currentMinimumStep,
G4double currentSafety 
) [pure virtual]

Implements G4VContinuousDiscreteProcess.

G4double G4VeLowEnergyLoss::GetLossWithFluct ( const G4DynamicParticle aParticle,
const G4MaterialCutsCouple couple,
G4double  MeanLoss,
G4double  step 
) [protected]

Definition at line 833 of file G4VeLowEnergyLoss.cc.

References e1Fluct, e1LogFluct, e2Fluct, e2LogFluct, f1Fluct, f2Fluct, G4Poisson(), G4UniformRand, G4Material::GetElectronDensity(), G4IonisParamMat::GetEnergy0fluct(), G4IonisParamMat::GetEnergy1fluct(), G4IonisParamMat::GetEnergy2fluct(), G4IonisParamMat::GetF1fluct(), G4IonisParamMat::GetF2fluct(), G4MaterialCutsCouple::GetIndex(), G4Material::GetIonisation(), G4DynamicParticle::GetKineticEnergy(), G4IonisParamMat::GetLogEnergy1fluct(), G4IonisParamMat::GetLogEnergy2fluct(), G4IonisParamMat::GetLogMeanExcEnergy(), G4MaterialCutsCouple::GetMaterial(), G4IonisParamMat::GetMeanExcitationEnergy(), G4ProductionCutsTable::GetProductionCutsTable(), G4IonisParamMat::GetRateionexcfluct(), imat, ipotFluct, ipotLogFluct, lastMaterial, nmaxCont1, nmaxCont2, ParticleMass, and rateFluct.

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   // check if the material has changed ( cache mechanism)
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   //  G4double dp1;
00871   G4double dp3;
00872   G4double siga ;
00873 
00874   // shortcut for very very small loss
00875   if(MeanLoss < minLoss) return MeanLoss ;
00876 
00877   // get particle data
00878   G4double Tkin   = aParticle->GetKineticEnergy();
00879 
00880   //  G4cout << "MGP -- Fluc Tkin " << Tkin/keV << " keV "  << " MeanLoss = " << MeanLoss/keV << G4endl;
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   // G4cout << "MGP Particle mass " << ParticleMass/MeV << " Tm " << Tm << G4endl;
00889 
00890   if(Tm > threshold) Tm = threshold;
00891   beta2 = tau2/(tau1*tau1);
00892 
00893   // Gaussian fluctuation ?
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)             // very small Step
00919     {
00920       e0 = aMaterial->GetIonisation()->GetEnergy0fluct();
00921       // G4cout << "MGP e0 = " << e0/keV << G4endl;
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           // G4cout << "MGP very small step " << loss/keV << G4endl;
00938         }
00939       else
00940         {
00941           //      G4cout << "MGP old Tm = " << Tm << " " << ipotFluct << " " << e0 << G4endl;
00942           Tm = Tm-ipotFluct+e0 ;
00943 
00944           // MGP ---- workaround to avoid log argument<0, TO BE CHECKED
00945           if (Tm <= 0.)
00946             {
00947               loss = MeanLoss;
00948               p3 = 0;
00949               // G4cout << "MGP correction loss = MeanLoss " << loss/keV << G4endl;
00950             }
00951           else
00952             {
00953               a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(Tm/e0));
00954 
00955               // G4cout << "MGP new Tm = " << Tm << " " << ipotFluct << " " << e0 << " a3= " << a3 << G4endl;
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               //G4cout << "MGP p3 " << p3 << G4endl;
00965 
00966             }
00967               
00968           if(p3 > 0)
00969             {
00970               w = (Tm-e0)/Tm ;
00971               if(p3 > nmaxCont2)
00972                 {
00973                   // G4cout << "MGP dp3 " << dp3 << " p3 " << p3 << " " << nmaxCont2 << G4endl;
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               // G4cout << "MGP Corrfac = " << Corrfac << " e0 = " << e0/keV << " loss = " << loss/keV << G4endl;
00984             }        
00985         }
00986     }
00987 
00988   else                              // not so small Step
00989     {
00990       // excitation type 1
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       // excitation type 2
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       // smearing to avoid unphysical peaks
01011       if(p2 > 0)
01012         loss += (1.-2.*G4UniformRand())*e2Fluct;
01013       else if (loss>0.)
01014         loss += (1.-2.*G4UniformRand())*e1Fluct;   
01015       
01016       // ionisation .......................................
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 }

virtual G4double G4VeLowEnergyLoss::GetMeanFreePath ( const G4Track track,
G4double  previousStepSize,
G4ForceCondition condition 
) [pure virtual]

Implements G4VContinuousDiscreteProcess.

virtual G4VParticleChange* G4VeLowEnergyLoss::PostStepDoIt ( const G4Track track,
const G4Step Step 
) [pure virtual]

Reimplemented from G4VContinuousDiscreteProcess.

void G4VeLowEnergyLoss::SetEnlossFluc ( G4bool  value  )  [static]

Definition at line 120 of file G4VeLowEnergyLoss.cc.

References EnlossFlucFlag.

00121 {
00122    EnlossFlucFlag = value;
00123 }

void G4VeLowEnergyLoss::SetRndmStep ( G4bool  value  )  [static]

Definition at line 115 of file G4VeLowEnergyLoss.cc.

References rndmStepFlag.

00116 {
00117    rndmStepFlag   = value;
00118 }

void G4VeLowEnergyLoss::SetStepFunction ( G4double  c1,
G4double  c2 
) [static]

Definition at line 125 of file G4VeLowEnergyLoss.cc.

References c1lim, c2lim, c3lim, dRoverRange, and finalRange.

00126 {
00127    dRoverRange = c1;
00128    finalRange = c2;
00129    c1lim=dRoverRange;
00130    c2lim=2.*(1-dRoverRange)*finalRange;
00131    c3lim=-(1.-dRoverRange)*finalRange*finalRange;
00132 }


Field Documentation

G4double G4VeLowEnergyLoss::c1lim = dRoverRange [static, protected]

Definition at line 237 of file G4VeLowEnergyLoss.hh.

Referenced by SetStepFunction().

G4double G4VeLowEnergyLoss::c2lim = 2.*(1.-dRoverRange)*finalRange [static, protected]

Definition at line 237 of file G4VeLowEnergyLoss.hh.

Referenced by SetStepFunction().

G4double G4VeLowEnergyLoss::c3lim = -(1.-dRoverRange)*finalRange*finalRange [static, protected]

Definition at line 237 of file G4VeLowEnergyLoss.hh.

Referenced by SetStepFunction().

G4double G4VeLowEnergyLoss::dRoverRange = 20*perCent [static, protected]

Definition at line 234 of file G4VeLowEnergyLoss.hh.

Referenced by SetStepFunction().

G4double G4VeLowEnergyLoss::e1Fluct [protected]

Definition at line 125 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::e1LogFluct [protected]

Definition at line 126 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::e2Fluct [protected]

Definition at line 125 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::e2LogFluct [protected]

Definition at line 126 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4bool G4VeLowEnergyLoss::EnlossFlucFlag = true [static, protected]

Definition at line 240 of file G4VeLowEnergyLoss.hh.

Referenced by SetEnlossFluc().

G4double G4VeLowEnergyLoss::f1Fluct [protected]

Definition at line 125 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::f2Fluct [protected]

Definition at line 125 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::finalRange = 200*micrometer [static, protected]

Definition at line 236 of file G4VeLowEnergyLoss.hh.

Referenced by SetStepFunction().

G4int G4VeLowEnergyLoss::imat [protected]

Definition at line 124 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::ipotFluct [protected]

Definition at line 125 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::ipotLogFluct [protected]

Definition at line 126 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

const G4Material* G4VeLowEnergyLoss::lastMaterial [protected]

Definition at line 123 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::ltauhigh [static, protected]

Definition at line 231 of file G4VeLowEnergyLoss.hh.

G4double G4VeLowEnergyLoss::ltaulow [static, protected]

Definition at line 231 of file G4VeLowEnergyLoss.hh.

const G4int G4VeLowEnergyLoss::nmaxCont1 [protected]

Definition at line 128 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

const G4int G4VeLowEnergyLoss::nmaxCont2 [protected]

Definition at line 128 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::ParticleMass [static, protected]

Definition at line 231 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4double G4VeLowEnergyLoss::rateFluct [protected]

Definition at line 125 of file G4VeLowEnergyLoss.hh.

Referenced by GetLossWithFluct().

G4bool G4VeLowEnergyLoss::rndmStepFlag = false [static, protected]

Definition at line 239 of file G4VeLowEnergyLoss.hh.

Referenced by SetRndmStep().

G4double G4VeLowEnergyLoss::tauhigh [static, protected]

Definition at line 231 of file G4VeLowEnergyLoss.hh.

G4double G4VeLowEnergyLoss::taulow [static, protected]

Definition at line 231 of file G4VeLowEnergyLoss.hh.


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:41 2013 for Geant4 by  doxygen 1.4.7