Geant4: G4HadronElasticProcess.cc Source File

00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // Geant4 Hadron Elastic Scattering Process 
00029 // 
00030 // Created 26 July 2012 V.Ivanchenko from G4WHadronElasticProcess
00031 //  
00032 // Modified:
00033 // 14-Sep-12 M.Kelsey -- Pass subType code to base ctor
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& /*step*/)
00086 {
00087   theTotalResult->Clear();
00088   theTotalResult->Initialize(track);
00089   G4double weight = track.GetWeight();
00090   theTotalResult->ProposeWeight(weight);
00091 
00092   // For elastic scattering, _any_ result is considered an interaction
00093   ClearNumberOfInteractionLengthLeft();
00094 
00095   G4double kineticEnergy = track.GetKineticEnergy();
00096   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
00097   const G4ParticleDefinition* part = dynParticle->GetDefinition();
00098 
00099   // NOTE:  Very low energy scatters were causing numerical (FPE) errors
00100   //        in earlier releases; these limits have not been changed since.
00101   if (kineticEnergy <= lowestEnergy)   return theTotalResult;
00102 
00103   G4Material* material = track.GetMaterial();
00104   G4Nucleus* targNucleus = GetTargetNucleusPointer();
00105 
00106   // Select element
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   // Initialize the hadronic projectile from the track
00144   //  G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
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   // Check the result for catastrophic energy non-conservation
00173   // cannot be applied because is not guranteed that recoil 
00174   // nucleus is created
00175   // result = CheckResult(theProj, targNucleus, result);
00176 
00177   // directions
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   // energies  
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   // NOTE:  Very low energy scatters were causing numerical (FPE) errors
00198   //        in earlier releases; these limits have not been changed since.
00199   if(efinal <= lowestEnergy) {
00200     edep += efinal;
00201     efinal = 0.0;
00202   }
00203 
00204   // primary change
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   //G4cout << "Efinal= " << efinal << "  TrackStatus= " << status << G4endl;
00220 
00221   theTotalResult->SetNumberOfSecondaries(0);
00222 
00223   // recoil
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       // G4cout << "recoil " << pdir << G4endl;
00232       pdir.rotate(phi, it);
00233       pdir.rotateUz(indir);
00234       // G4cout << "recoil rotated " << pdir << G4endl;
00235       p->SetMomentumDirection(pdir);
00236 
00237       // in elastic scattering time and weight are not changed
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