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 #include "G4ErrorEnergyLoss.hh"
00028 #include "G4ErrorPropagatorData.hh"
00029 #include "G4EnergyLossForExtrapolator.hh"
00030
00031
00032 G4ErrorEnergyLoss::G4ErrorEnergyLoss(const G4String& processName,
00033 G4ProcessType type)
00034 : G4VContinuousProcess(processName, type)
00035 {
00036 if (verboseLevel>2) {
00037 G4cout << GetProcessName() << " is created " << G4endl;
00038 }
00039
00040 theELossForExtrapolator = new G4EnergyLossForExtrapolator;
00041
00042 theStepLimit = 1.;
00043 }
00044
00045
00046 void G4ErrorEnergyLoss::InstantiateEforExtrapolator()
00047 {}
00048
00049
00050
00051 G4ErrorEnergyLoss::~G4ErrorEnergyLoss()
00052 {
00053 delete theELossForExtrapolator;
00054 }
00055
00056
00057 G4VParticleChange*
00058 G4ErrorEnergyLoss::AlongStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00059 {
00060 aParticleChange.Initialize(aTrack);
00061
00062 G4ErrorPropagatorData* g4edata = G4ErrorPropagatorData::GetErrorPropagatorData();
00063
00064 G4double kinEnergyStart = aTrack.GetKineticEnergy();
00065 G4double step_length = aStep.GetStepLength();
00066
00067 const G4Material* aMaterial = aTrack.GetMaterial();
00068 const G4ParticleDefinition* aParticleDef = aTrack.GetDynamicParticle()->GetDefinition();
00069 G4double kinEnergyEnd = kinEnergyStart;
00070
00071 if( g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropBackwards) ) {
00072 kinEnergyEnd = theELossForExtrapolator->EnergyBeforeStep( kinEnergyStart,
00073 step_length,
00074 aMaterial,
00075 aParticleDef );
00076 G4double kinEnergyHalfStep = kinEnergyStart - (kinEnergyStart-kinEnergyEnd)/2.;
00077
00078 #ifdef G4VERBOSE
00079 if(G4ErrorPropagatorData::verbose() >= 3 )
00080 G4cout << " G4ErrorEnergyLoss FWD end " << kinEnergyEnd
00081 << " halfstep " << kinEnergyHalfStep << G4endl;
00082 #endif
00083
00084
00085 kinEnergyEnd = theELossForExtrapolator->EnergyBeforeStep( kinEnergyHalfStep,
00086 step_length,
00087 aMaterial,
00088 aParticleDef );
00089 kinEnergyEnd = kinEnergyStart - (kinEnergyHalfStep - kinEnergyEnd );
00090 }else if( g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropForwards) ) {
00091 kinEnergyEnd = theELossForExtrapolator->EnergyAfterStep( kinEnergyStart,
00092 step_length,
00093 aMaterial,
00094 aParticleDef );
00095 G4double kinEnergyHalfStep = kinEnergyStart - (kinEnergyStart-kinEnergyEnd)/2.;
00096 #ifdef G4VERBOSE
00097 if(G4ErrorPropagatorData::verbose() >= 3 )
00098 G4cout << " G4ErrorEnergyLoss BCKD end " << kinEnergyEnd
00099 << " halfstep " << kinEnergyHalfStep << G4endl;
00100 #endif
00101
00102
00103 kinEnergyEnd = theELossForExtrapolator->EnergyAfterStep( kinEnergyHalfStep,
00104 step_length,
00105 aMaterial,
00106 aParticleDef );
00107 kinEnergyEnd = kinEnergyStart - (kinEnergyHalfStep - kinEnergyEnd );
00108 }
00109
00110 G4double edepo = kinEnergyEnd - kinEnergyStart;
00111
00112 #ifdef G4VERBOSE
00113 if( G4ErrorPropagatorData::verbose() >= 2 )
00114 G4cout << "AlongStepDoIt Estart= " << kinEnergyStart << " Eend " << kinEnergyEnd
00115 << " Ediff " << kinEnergyStart-kinEnergyEnd << " step= " << step_length
00116 << " mate= " << aMaterial->GetName()
00117 << " particle= " << aParticleDef->GetParticleName() << G4endl;
00118 #endif
00119
00120 aParticleChange.ClearDebugFlag();
00121 aParticleChange.ProposeLocalEnergyDeposit( edepo );
00122 aParticleChange.SetNumberOfSecondaries(0);
00123
00124 aParticleChange.ProposeEnergy( kinEnergyEnd );
00125
00126 return &aParticleChange;
00127 }
00128
00129
00130
00131 G4double G4ErrorEnergyLoss::GetContinuousStepLimit(const G4Track& aTrack,
00132 G4double ,
00133 G4double currentMinimumStep,
00134 G4double& )
00135 {
00136 G4double Step = DBL_MAX;
00137 if( theStepLimit != 1. ) {
00138 G4double kinEnergyStart = aTrack.GetKineticEnergy();
00139 G4double kinEnergyLoss = kinEnergyStart;
00140 const G4Material* aMaterial = aTrack.GetMaterial();
00141 const G4ParticleDefinition* aParticleDef = aTrack.GetDynamicParticle()->GetDefinition();
00142 G4ErrorPropagatorData* g4edata = G4ErrorPropagatorData::GetErrorPropagatorData();
00143 if( g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropBackwards) ) {
00144 kinEnergyLoss = - kinEnergyStart +
00145 theELossForExtrapolator->EnergyBeforeStep( kinEnergyStart, currentMinimumStep,
00146 aMaterial, aParticleDef );
00147 }else if( g4edata->GetMode() == G4ErrorMode(G4ErrorMode_PropForwards) ) {
00148 kinEnergyLoss = kinEnergyStart -
00149 theELossForExtrapolator->EnergyAfterStep( kinEnergyStart, currentMinimumStep,
00150 aMaterial, aParticleDef );
00151 }
00152 #ifdef G4VERBOSE
00153 if(G4ErrorPropagatorData::verbose() >= 3 )
00154 G4cout << " G4ErrorEnergyLoss: currentMinimumStep " <<currentMinimumStep
00155 << " kinEnergyLoss " << kinEnergyLoss
00156 << " kinEnergyStart " << kinEnergyStart << G4endl;
00157 #endif
00158 if( kinEnergyLoss / kinEnergyStart > theStepLimit ) {
00159 Step = theStepLimit / (kinEnergyLoss / kinEnergyStart) * currentMinimumStep;
00160 #ifdef G4VERBOSE
00161 if(G4ErrorPropagatorData::verbose() >= 2 )
00162 G4cout << " G4ErrorEnergyLoss: limiting Step " << Step
00163 << " energy loss fraction " << kinEnergyLoss / kinEnergyStart
00164 << " > " << theStepLimit << G4endl;
00165 #endif
00166 }
00167 }
00168
00169 return Step;
00170
00171 }