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 #include "G4ParticleChangeForLoss.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "G4Track.hh"
00047 #include "G4Step.hh"
00048 #include "G4DynamicParticle.hh"
00049 #include "G4ExceptionSeverity.hh"
00050
00051
00052 G4ParticleChangeForLoss::G4ParticleChangeForLoss()
00053 : G4VParticleChange(), currentTrack(0), proposedKinEnergy(0.),
00054 lowEnergyLimit(1.0*eV), currentCharge(0.)
00055 {
00056 theSteppingControlFlag = NormalCondition;
00057 debugFlag = false;
00058 #ifdef G4VERBOSE
00059 if (verboseLevel>2) {
00060 G4cout << "G4ParticleChangeForLoss::G4ParticleChangeForLoss() " << G4endl;
00061 }
00062 #endif
00063 }
00064
00065 G4ParticleChangeForLoss::~G4ParticleChangeForLoss()
00066 {
00067 #ifdef G4VERBOSE
00068 if (verboseLevel>2) {
00069 G4cout << "G4ParticleChangeForLoss::~G4ParticleChangeForLoss() " << G4endl;
00070 }
00071 #endif
00072 }
00073
00074 G4ParticleChangeForLoss::
00075 G4ParticleChangeForLoss(const G4ParticleChangeForLoss &right)
00076 : G4VParticleChange(right)
00077 {
00078 if (verboseLevel>1) {
00079 G4cout << "G4ParticleChangeForLoss:: copy constructor is called " << G4endl;
00080 }
00081 currentTrack = right.currentTrack;
00082 proposedKinEnergy = right.proposedKinEnergy;
00083 lowEnergyLimit = right.lowEnergyLimit;
00084 currentCharge = right.currentCharge;
00085 proposedMomentumDirection = right.proposedMomentumDirection;
00086 }
00087
00088
00089 G4ParticleChangeForLoss & G4ParticleChangeForLoss::operator=(
00090 const G4ParticleChangeForLoss &right)
00091 {
00092 #ifdef G4VERBOSE
00093 if (verboseLevel>1) {
00094 G4cout << "G4ParticleChangeForLoss:: assignment operator is called " << G4endl;
00095 }
00096 #endif
00097
00098 if (this != &right) {
00099 if (theNumberOfSecondaries>0) {
00100 #ifdef G4VERBOSE
00101 if (verboseLevel>0) {
00102 G4cout << "G4ParticleChangeForLoss: assignment operator Warning ";
00103 G4cout << "theListOfSecondaries is not empty ";
00104 }
00105 #endif
00106 for (G4int index= 0; index<theNumberOfSecondaries; index++){
00107 if ( (*theListOfSecondaries)[index] ) delete (*theListOfSecondaries)[index] ;
00108 }
00109 }
00110 delete theListOfSecondaries;
00111 theListOfSecondaries = new G4TrackFastVector();
00112 theNumberOfSecondaries = right.theNumberOfSecondaries;
00113 for (G4int index = 0; index<theNumberOfSecondaries; index++){
00114 G4Track* newTrack = new G4Track(*((*right.theListOfSecondaries)[index] ));
00115 theListOfSecondaries->SetElement(index, newTrack); }
00116
00117 theStatusChange = right.theStatusChange;
00118 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
00119 theSteppingControlFlag = right.theSteppingControlFlag;
00120 theParentWeight = right.theParentWeight;
00121 isParentWeightProposed = right.isParentWeightProposed;
00122 fSetSecondaryWeightByProcess = right.fSetSecondaryWeightByProcess;
00123
00124 currentTrack = right.currentTrack;
00125 proposedKinEnergy = right.proposedKinEnergy;
00126 currentCharge = right.currentCharge;
00127 proposedMomentumDirection = right.proposedMomentumDirection;
00128 }
00129 return *this;
00130 }
00131
00132
00133
00134
00135
00136 void G4ParticleChangeForLoss::DumpInfo() const
00137 {
00138
00139 G4VParticleChange::DumpInfo();
00140
00141 G4int oldprc = G4cout.precision(3);
00142 G4cout << " Charge (eplus) : "
00143 << std::setw(20) << currentCharge/eplus
00144 << G4endl;
00145 G4cout << " Kinetic Energy (MeV): "
00146 << std::setw(20) << proposedKinEnergy/MeV
00147 << G4endl;
00148 G4cout << " Momentum Direct - x : "
00149 << std::setw(20) << proposedMomentumDirection.x()
00150 << G4endl;
00151 G4cout << " Momentum Direct - y : "
00152 << std::setw(20) << proposedMomentumDirection.y()
00153 << G4endl;
00154 G4cout << " Momentum Direct - z : "
00155 << std::setw(20) << proposedMomentumDirection.z()
00156 << G4endl;
00157 G4cout.precision(oldprc);
00158 }
00159
00160 G4bool G4ParticleChangeForLoss::CheckIt(const G4Track& aTrack)
00161 {
00162 G4bool itsOK = true;
00163 G4bool exitWithError = false;
00164
00165 G4double accuracy;
00166
00167
00168 accuracy = ( proposedKinEnergy - aTrack.GetKineticEnergy())/MeV;
00169 if (accuracy > accuracyForWarning) {
00170 itsOK = false;
00171 exitWithError = (accuracy > accuracyForException);
00172 #ifdef G4VERBOSE
00173 G4cout << "G4ParticleChangeForLoss::CheckIt: ";
00174 G4cout << "KinEnergy become larger than the initial value!"
00175 << " Difference: " << accuracy << "[MeV] " <<G4endl;
00176 G4cout << aTrack.GetDefinition()->GetParticleName()
00177 << " E=" << aTrack.GetKineticEnergy()/MeV
00178 << " pos=" << aTrack.GetPosition().x()/m
00179 << ", " << aTrack.GetPosition().y()/m
00180 << ", " << aTrack.GetPosition().z()/m
00181 <<G4endl;
00182 #endif
00183 }
00184
00185
00186 #ifdef G4VERBOSE
00187 if (!itsOK) DumpInfo();
00188 #endif
00189
00190
00191 if (exitWithError) {
00192 G4Exception("G4ParticleChangeForLoss::CheckIt",
00193 "TRACK004", EventMustBeAborted,
00194 "energy was illegal");
00195 }
00196
00197
00198 if (!itsOK) {
00199 proposedKinEnergy = aTrack.GetKineticEnergy();
00200 }
00201
00202 itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
00203 return itsOK;
00204 }
00205
00206
00207
00208
00209
00210 G4Step* G4ParticleChangeForLoss::UpdateStepForAlongStep(G4Step* pStep)
00211 {
00212 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
00213 G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
00214 G4Track* pTrack = pStep->GetTrack();
00215
00216
00217 G4double preKinEnergy = pPreStepPoint->GetKineticEnergy();
00218 G4double kinEnergy = pPostStepPoint->GetKineticEnergy() +
00219 (proposedKinEnergy - preKinEnergy);
00220
00221
00222 if (kinEnergy < lowEnergyLimit) {
00223 theLocalEnergyDeposit += kinEnergy;
00224 kinEnergy = 0.0;
00225 pPostStepPoint->SetVelocity(0.0);
00226 } else {
00227 pPostStepPoint->SetCharge( currentCharge );
00228
00229 pTrack->SetKineticEnergy(kinEnergy);
00230 pPostStepPoint->SetVelocity(pTrack->CalculateVelocity());
00231 pTrack->SetKineticEnergy(preKinEnergy);
00232 }
00233 pPostStepPoint->SetKineticEnergy( kinEnergy );
00234
00235 if (isParentWeightProposed ){
00236 pPostStepPoint->SetWeight( theParentWeight );
00237 }
00238
00239 pStep->AddTotalEnergyDeposit( theLocalEnergyDeposit );
00240 pStep->AddNonIonizingEnergyDeposit( theNonIonizingEnergyDeposit );
00241 return pStep;
00242 }
00243
00244 G4Step* G4ParticleChangeForLoss::UpdateStepForPostStep(G4Step* pStep)
00245 {
00246 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
00247 G4Track* pTrack = pStep->GetTrack();
00248
00249 pPostStepPoint->SetCharge( currentCharge );
00250 pPostStepPoint->SetMomentumDirection( proposedMomentumDirection );
00251 pPostStepPoint->SetKineticEnergy( proposedKinEnergy );
00252 pTrack->SetKineticEnergy( proposedKinEnergy );
00253 if(proposedKinEnergy > 0.0) {
00254 pPostStepPoint->SetVelocity(pTrack->CalculateVelocity());
00255 } else {
00256 pPostStepPoint->SetVelocity(0.0);
00257 }
00258 pPostStepPoint->SetPolarization( proposedPolarization );
00259
00260 if (isParentWeightProposed ){
00261 pPostStepPoint->SetWeight( theParentWeight );
00262 }
00263
00264 pStep->AddTotalEnergyDeposit( theLocalEnergyDeposit );
00265 pStep->AddNonIonizingEnergyDeposit( theNonIonizingEnergyDeposit );
00266 return pStep;
00267 }
00268