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 #include <iostream>
00036 #include <typeinfo>
00037
00038 #include "G4HadronElasticProcess.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4Nucleus.hh"
00041 #include "G4ProcessManager.hh"
00042 #include "G4CrossSectionDataStore.hh"
00043 #include "G4HadronElasticDataSet.hh"
00044 #include "G4ProductionCutsTable.hh"
00045 #include "G4HadronicException.hh"
00046 #include "G4HadronicDeprecate.hh"
00047
00048 G4HadronElasticProcess::G4HadronElasticProcess(const G4String& pName)
00049 : G4HadronicProcess(pName, fHadronElastic), isInitialised(false)
00050 {
00051 AddDataSet(new G4HadronElasticDataSet);
00052 lowestEnergy = 1.*keV;
00053 }
00054
00055 G4HadronElasticProcess::~G4HadronElasticProcess()
00056 {}
00057
00058 void G4HadronElasticProcess::Description() const
00059 {
00060 char* dirName = getenv("G4PhysListDocDir");
00061 if (dirName) {
00062 std::ofstream outFile;
00063 G4String outFileName = GetProcessName() + ".html";
00064 G4String pathName = G4String(dirName) + "/" + outFileName;
00065 outFile.open(pathName);
00066 outFile << "<html>\n";
00067 outFile << "<head>\n";
00068
00069 outFile << "<title>Description of G4HadronElasticProcess</title>\n";
00070 outFile << "</head>\n";
00071 outFile << "<body>\n";
00072
00073 outFile << "G4HadronElasticProcess handles the elastic scattering of\n"
00074 << "hadrons by invoking one or more hadronic models and one or\n"
00075 << "more hadronic cross sections.\n";
00076
00077 outFile << "</body>\n";
00078 outFile << "</html>\n";
00079 outFile.close();
00080 }
00081 }
00082
00083 G4VParticleChange*
00084 G4HadronElasticProcess::PostStepDoIt(const G4Track& track,
00085 const G4Step& )
00086 {
00087 theTotalResult->Clear();
00088 theTotalResult->Initialize(track);
00089 G4double weight = track.GetWeight();
00090 theTotalResult->ProposeWeight(weight);
00091
00092
00093 ClearNumberOfInteractionLengthLeft();
00094
00095 G4double kineticEnergy = track.GetKineticEnergy();
00096 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
00097 const G4ParticleDefinition* part = dynParticle->GetDefinition();
00098
00099
00100
00101 if (kineticEnergy <= lowestEnergy) return theTotalResult;
00102
00103 G4Material* material = track.GetMaterial();
00104 G4Nucleus* targNucleus = GetTargetNucleusPointer();
00105
00106
00107 G4Element* elm = 0;
00108 try
00109 {
00110 elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,
00111 *targNucleus);
00112 }
00113 catch(G4HadronicException & aR)
00114 {
00115 G4ExceptionDescription ed;
00116 DumpState(track,"SampleZandA",ed);
00117 ed << " PostStepDoIt failed on element selection" << G4endl;
00118 G4Exception("G4HadronElasticProcess::PostStepDoIt", "had003",
00119 FatalException, ed);
00120 }
00121 G4HadronicInteraction* hadi = 0;
00122 try
00123 {
00124 hadi = ChooseHadronicInteraction( kineticEnergy, material, elm );
00125 }
00126 catch(G4HadronicException & aE)
00127 {
00128 G4ExceptionDescription ed;
00129 ed << "Target element "<< elm->GetName()<<" Z= "
00130 << targNucleus->GetZ_asInt() << " A= "
00131 << targNucleus->GetA_asInt() << G4endl;
00132 DumpState(track,"ChooseHadronicInteraction",ed);
00133 ed << " No HadronicInteraction found out" << G4endl;
00134 G4Exception("G4HadronElasticProcess::PostStepDoIt", "had005",
00135 FatalException, ed);
00136 }
00137
00138 size_t idx = track.GetMaterialCutsCouple()->GetIndex();
00139 G4double tcut = (*(G4ProductionCutsTable::GetProductionCutsTable()
00140 ->GetEnergyCutsVector(3)))[idx];
00141 hadi->SetRecoilEnergyThreshold(tcut);
00142
00143
00144
00145 G4HadProjectile theProj(track);
00146 if(verboseLevel>1) {
00147 G4cout << "G4HadronElasticProcess::PostStepDoIt for "
00148 << part->GetParticleName()
00149 << " in " << material->GetName()
00150 << " Target Z= " << targNucleus->GetZ_asInt()
00151 << " A= " << targNucleus->GetA_asInt() << G4endl;
00152 }
00153
00154 G4HadFinalState* result = 0;
00155 try
00156 {
00157 result = hadi->ApplyYourself( theProj, *targNucleus);
00158 }
00159 catch(G4HadronicException aR)
00160 {
00161 G4ExceptionDescription ed;
00162 ed << "Call for " << hadi->GetModelName() << G4endl;
00163 ed << "Target element "<< elm->GetName()<<" Z= "
00164 << targNucleus->GetZ_asInt()
00165 << " A= " << targNucleus->GetA_asInt() << G4endl;
00166 DumpState(track,"ApplyYourself",ed);
00167 ed << " ApplyYourself failed" << G4endl;
00168 G4Exception("G4HadronElasticProcess::PostStepDoIt", "had006",
00169 FatalException, ed);
00170 }
00171
00172
00173
00174
00175
00176
00177
00178 G4ThreeVector indir = track.GetMomentumDirection();
00179 G4double phi = CLHEP::twopi*G4UniformRand();
00180 G4ThreeVector it(0., 0., 1.);
00181 G4ThreeVector outdir = result->GetMomentumChange();
00182
00183 if(verboseLevel>1) {
00184 G4cout << "Efin= " << result->GetEnergyChange()
00185 << " de= " << result->GetLocalEnergyDeposit()
00186 << " nsec= " << result->GetNumberOfSecondaries()
00187 << " dir= " << outdir
00188 << G4endl;
00189 }
00190
00191
00192 G4double edep = result->GetLocalEnergyDeposit();
00193 G4double efinal = result->GetEnergyChange();
00194 if(efinal < 0.0) { efinal = 0.0; }
00195 if(edep < 0.0) { edep = 0.0; }
00196
00197
00198
00199 if(efinal <= lowestEnergy) {
00200 edep += efinal;
00201 efinal = 0.0;
00202 }
00203
00204
00205 theTotalResult->ProposeEnergy(efinal);
00206
00207 G4TrackStatus status = track.GetTrackStatus();
00208 if(efinal > 0.0) {
00209 outdir.rotate(phi, it);
00210 outdir.rotateUz(indir);
00211 theTotalResult->ProposeMomentumDirection(outdir);
00212 } else {
00213 if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
00214 { status = fStopButAlive; }
00215 else { status = fStopAndKill; }
00216 theTotalResult->ProposeTrackStatus(status);
00217 }
00218
00219
00220
00221 theTotalResult->SetNumberOfSecondaries(0);
00222
00223
00224 if(result->GetNumberOfSecondaries() > 0) {
00225 G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
00226
00227 if(p->GetKineticEnergy() > tcut) {
00228 theTotalResult->SetNumberOfSecondaries(1);
00229 G4ThreeVector pdir = p->GetMomentumDirection();
00230
00232 pdir.rotate(phi, it);
00233 pdir.rotateUz(indir);
00234
00235 p->SetMomentumDirection(pdir);
00236
00237
00238 G4Track* t = new G4Track(p, track.GetGlobalTime(),
00239 track.GetPosition());
00240 t->SetWeight(weight);
00241 t->SetTouchableHandle(track.GetTouchableHandle());
00242 theTotalResult->AddSecondary(t);
00243
00244 } else {
00245 edep += p->GetKineticEnergy();
00246 delete p;
00247 }
00248 }
00249 theTotalResult->ProposeLocalEnergyDeposit(edep);
00250 theTotalResult->ProposeNonIonizingEnergyDeposit(edep);
00251 result->Clear();
00252
00253 return theTotalResult;
00254 }
00255
00256 void
00257 G4HadronElasticProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
00258 {
00259 if(!isInitialised) {
00260 isInitialised = true;
00261 if(G4Neutron::Neutron() == &part) { lowestEnergy = 1.e-6*eV; }
00262 }
00263 G4HadronicProcess::PreparePhysicsTable(part);
00264 }
00265
00266 void G4HadronElasticProcess::SetLowestEnergy(G4double val)
00267 {
00268 lowestEnergy = val;
00269 }
00270
00271 void
00272 G4HadronElasticProcess::SetLowestEnergyNeutron(G4double val)
00273 {
00274 lowestEnergy = val;
00275 G4HadronicDeprecate("G4HadronElasticProcess::SetLowestEnergyNeutron()");
00276 }
00277