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 #include "G4HadronStoppingProcess.hh"
00051 #include "G4HadronicProcessStore.hh"
00052 #include "G4HadronicProcessType.hh"
00053 #include "G4EmCaptureCascade.hh"
00054 #include "G4Nucleus.hh"
00055 #include "G4HadFinalState.hh"
00056 #include "G4HadProjectile.hh"
00057 #include "G4HadSecondary.hh"
00058 #include "G4Material.hh"
00059 #include "G4Element.hh"
00060
00061
00062
00063 G4HadronStoppingProcess::G4HadronStoppingProcess(const G4String& name)
00064 : G4HadronicProcess(name, fHadronAtRest)
00065 {
00066
00067 enableAtRestDoIt = true;
00068 enablePostStepDoIt = false;
00069
00070 fElementSelector = new G4ElementSelector();
00071 fEmCascade = new G4EmCaptureCascade();
00072 fBoundDecay = 0;
00073 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
00074 }
00075
00076
00077
00078 G4HadronStoppingProcess::~G4HadronStoppingProcess()
00079 {
00080 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
00081 delete fElementSelector;
00082
00083 }
00084
00085
00086
00087 G4bool G4HadronStoppingProcess::IsApplicable(const G4ParticleDefinition& p)
00088 {
00089 return (p.GetPDGCharge() < 0.0);
00090 }
00091
00092
00093
00094 void G4HadronStoppingProcess::PreparePhysicsTable(const G4ParticleDefinition& p)
00095 {
00096 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00097 }
00098
00099
00100
00101 void G4HadronStoppingProcess::BuildPhysicsTable(const G4ParticleDefinition& p)
00102 {
00103 G4HadronicProcessStore::Instance()->PrintInfo(&p);
00104 }
00105
00106
00107
00108 G4double G4HadronStoppingProcess::AtRestGetPhysicalInteractionLength(
00109 const G4Track&, G4ForceCondition* condition)
00110 {
00111 *condition = NotForced;
00112 return 0.0;
00113 }
00114
00115
00116
00117 G4double G4HadronStoppingProcess::PostStepGetPhysicalInteractionLength(
00118 const G4Track&, G4double, G4ForceCondition* condition)
00119 {
00120 *condition = NotForced;
00121 return DBL_MAX;
00122 }
00123
00124
00125
00126 G4VParticleChange* G4HadronStoppingProcess::AtRestDoIt(const G4Track& track,
00127 const G4Step&)
00128 {
00129
00130 theTotalResult->Initialize(track);
00131
00132 G4Nucleus* nucleus = GetTargetNucleusPointer();
00133 G4Element* elm = fElementSelector->SelectZandA(track, nucleus);
00134
00135 G4HadFinalState* result = 0;
00136 thePro.Initialise(track);
00137 G4double time0 = track.GetGlobalTime();
00138 G4bool nuclearCapture = true;
00139
00140
00141
00142
00143
00144 result = fEmCascade->ApplyYourself(thePro, *nucleus);
00145 G4double ebound = result->GetLocalEnergyDeposit();
00146 G4double edep = 0.0;
00147 G4int nSecondaries = result->GetNumberOfSecondaries();
00148
00149
00150
00151
00152
00153
00154 thePro.SetBoundEnergy(ebound);
00155 if(fBoundDecay) {
00156 G4HadFinalState* resultDecay =
00157 fBoundDecay->ApplyYourself(thePro, *nucleus);
00158 G4int n = resultDecay->GetNumberOfSecondaries();
00159 if(0 < n) {
00160 nSecondaries += n;
00161 result->AddSecondaries(resultDecay);
00162 }
00163 if(resultDecay->GetStatusChange() == stopAndKill) {
00164 nuclearCapture = false;
00165 }
00166 resultDecay->Clear();
00167 }
00168
00169 if(nuclearCapture) {
00170
00171
00172 G4HadronicInteraction* model = 0;
00173 try {
00174 model = ChooseHadronicInteraction(0.0, track.GetMaterial(), elm);
00175 }
00176 catch(G4HadronicException & aE) {
00177 G4ExceptionDescription ed;
00178 ed << "Target element "<<elm->GetName()<<" Z= "
00179 << nucleus->GetZ_asInt() << " A= "
00180 << nucleus->GetA_asInt() << G4endl;
00181 DumpState(track,"ChooseHadronicInteraction",ed);
00182 ed << " No HadronicInteraction found out" << G4endl;
00183 G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had005",
00184 FatalException, ed);
00185 }
00186
00187 G4HadFinalState* resultNuc = 0;
00188 G4int reentryCount = 0;
00189 do {
00190
00191
00192
00193 try {
00194 resultNuc = model->ApplyYourself(thePro, *nucleus);
00195 ++reentryCount;
00196 }
00197 catch(G4HadronicException aR) {
00198 G4ExceptionDescription ed;
00199 ed << "Call for " << model->GetModelName() << G4endl;
00200 ed << "Target element "<<elm->GetName()<<" Z= "
00201 << nucleus->GetZ_asInt()
00202 << " A= " << nucleus->GetA_asInt() << G4endl;
00203 DumpState(track,"ApplyYourself",ed);
00204 ed << " ApplyYourself failed" << G4endl;
00205 G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
00206 FatalException, ed);
00207 }
00208
00209
00210 resultNuc = CheckResult(thePro, *nucleus, resultNuc);
00211
00212 if(reentryCount>100) {
00213 G4ExceptionDescription ed;
00214 ed << "Call for " << model->GetModelName() << G4endl;
00215 ed << "Target element "<<elm->GetName()<<" Z= "
00216 << nucleus->GetZ_asInt()
00217 << " A= " << nucleus->GetA_asInt() << G4endl;
00218 DumpState(track,"ApplyYourself",ed);
00219 ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
00220 G4Exception("G4HadronStoppingProcess::AtRestDoIt", "had006",
00221 FatalException, ed);
00222 }
00223 }
00224 while(!resultNuc);
00225
00226 edep = resultNuc->GetLocalEnergyDeposit();
00227 nSecondaries += resultNuc->GetNumberOfSecondaries();
00228 result->AddSecondaries(resultNuc);
00229 resultNuc->Clear();
00230 }
00231
00232
00233
00234 theTotalResult->ProposeTrackStatus(fStopAndKill);
00235 theTotalResult->ProposeLocalEnergyDeposit(edep);
00236 theTotalResult->SetNumberOfSecondaries(nSecondaries);
00237 G4double w = track.GetWeight();
00238 theTotalResult->ProposeWeight(w);
00239 for(G4int i=0; i<nSecondaries; ++i) {
00240 G4HadSecondary* sec = result->GetSecondary(i);
00241
00242
00243 G4double time = sec->GetTime();
00244 if(time < 0.0) { time = 0.0; }
00245 time += time0;
00246
00247
00248 G4Track* t = new G4Track(sec->GetParticle(),
00249 time,
00250 track.GetPosition());
00251 t->SetWeight(w*sec->GetWeight());
00252 t->SetTouchableHandle(track.GetTouchableHandle());
00253 theTotalResult->AddSecondary(t);
00254 }
00255 result->Clear();
00256
00257 if (epReportLevel != 0) {
00258 CheckEnergyMomentumConservation(track, *nucleus);
00259 }
00260 return theTotalResult;
00261 }
00262
00263
00264
00265 void G4HadronStoppingProcess::ProcessDescription(std::ostream& outFile) const
00266 {
00267 outFile << "Base process for negatively charged particle capture at rest.\n";
00268 }
00269
00270