Geant4: G4WHadronElasticProcess.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 21 April 2006 V.Ivanchenko
00031 //  
00032 // Modified:
00033 // 24.04.06 V.Ivanchenko add neutron scattering on hydrogen from CHIPS
00034 // 07.06.06 V.Ivanchenko fix problem of rotation of final state
00035 // 25.07.06 V.Ivanchenko add 19 MeV low energy for CHIPS
00036 // 26.09.06 V.Ivanchenko add lowestEnergy
00037 // 20.10.06 V.Ivanchenko initialise lowestEnergy=0 for neitrals, eV for charged
00038 // 23.01.07 V.Ivanchenko add cross section interfaces with Z and A
00039 // 02.05.07 V.Ivanchenko add He3
00040 // 13.01.10 M.Kosov: Commented not used G4QElasticCrossSection & G4QCHIPSWorld
00041 // 24.02.11 V.Ivanchenko use particle name in IfApplicable, 
00042 //                       added anti particles for light ions
00043 // 07.09.11 M.Kelsey: Follow chanhe to G4HadFinalState interface
00044 // 14.09.12 M.Kelsey: Pass subType code to base ctor
00045 //
00046 
00047 #include <iostream>
00048 #include <typeinfo>
00049 
00050 #include "G4WHadronElasticProcess.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4Nucleus.hh"
00053 #include "G4ProcessManager.hh"
00054 #include "G4CrossSectionDataStore.hh"
00055 #include "G4HadronElasticDataSet.hh"
00056 #include "G4ProductionCutsTable.hh"
00057 #include "G4HadronicException.hh"
00058 #include "G4HadronicDeprecate.hh"
00059 
00060 G4WHadronElasticProcess::G4WHadronElasticProcess(const G4String& pName)
00061   : G4HadronicProcess(pName, fHadronElastic) {
00062   AddDataSet(new G4HadronElasticDataSet);
00063   theNeutron  = G4Neutron::Neutron();
00064   lowestEnergy = 1.*keV;
00065   lowestEnergyNeutron = 1.e-6*eV;
00066   G4HadronicDeprecate("G4WHadronElasticProcess");
00067 }
00068 
00069 G4WHadronElasticProcess::~G4WHadronElasticProcess()
00070 {}
00071 
00072 void G4WHadronElasticProcess::Description() const
00073 {
00074   char* dirName = getenv("G4PhysListDocDir");
00075   if (dirName) {
00076     std::ofstream outFile;
00077     G4String outFileName = GetProcessName() + ".html";
00078     G4String pathName = G4String(dirName) + "/" + outFileName;
00079     outFile.open(pathName);
00080     outFile << "<html>\n";
00081     outFile << "<head>\n";
00082 
00083     outFile << "<title>Description of G4WHadronElasticProcess</title>\n";
00084     outFile << "</head>\n";
00085     outFile << "<body>\n";
00086 
00087     outFile << "G4WHadronElasticProcess handles the elastic scattering of\n"
00088             << "hadrons by invoking one or more hadronic models and one or\n"
00089             << "more hadronic cross sections.\n";
00090 
00091     outFile << "</body>\n";
00092     outFile << "</html>\n";
00093     outFile.close();
00094   }
00095 }
00096 
00097 
00098 G4VParticleChange* 
00099 G4WHadronElasticProcess::PostStepDoIt(const G4Track& track, 
00100                                       const G4Step& /*step*/)
00101 {
00102   theTotalResult->Clear();
00103   theTotalResult->Initialize(track);
00104   G4double weight = track.GetWeight();
00105   theTotalResult->ProposeWeight(weight);
00106 
00107   // For elastic scattering, _any_ result is considered an interaction
00108   ClearNumberOfInteractionLengthLeft();
00109 
00110   G4double kineticEnergy = track.GetKineticEnergy();
00111   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
00112   const G4ParticleDefinition* part = dynParticle->GetDefinition();
00113 
00114   // NOTE:  Very low energy scatters were causing numerical (FPE) errors
00115   //        in earlier releases; these limits have not been changed since.
00116   if (part == theNeutron) {
00117     if(kineticEnergy <= lowestEnergyNeutron) return theTotalResult;
00118   } else if(kineticEnergy <= lowestEnergy)   return theTotalResult;
00119 
00120   G4Material* material = track.GetMaterial();
00121   G4Nucleus* targNucleus = GetTargetNucleusPointer();
00122 
00123   // Select element
00124   G4Element* elm = 0;
00125   try
00126     {
00127       elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material, 
00128                                                     *targNucleus);
00129     }
00130   catch(G4HadronicException & aR)
00131     {
00132       G4ExceptionDescription ed;
00133       DumpState(track,"SampleZandA",ed); 
00134       ed << " PostStepDoIt failed on element selection" << G4endl;
00135       G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had003", 
00136                   FatalException, ed);
00137     }
00138   G4HadronicInteraction* hadi = 0;
00139   try
00140     {
00141       hadi = ChooseHadronicInteraction( kineticEnergy, material, elm );
00142     }
00143   catch(G4HadronicException & aE)
00144     {
00145       G4ExceptionDescription ed;
00146       ed << "Target element "<< elm->GetName()<<"  Z= " 
00147          << targNucleus->GetZ_asInt() << "  A= " 
00148          << targNucleus->GetA_asInt() << G4endl;
00149       DumpState(track,"ChooseHadronicInteraction",ed);
00150       ed << " No HadronicInteraction found out" << G4endl;
00151       G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had005", 
00152                   FatalException, ed);
00153     }
00154 
00155   size_t idx = track.GetMaterialCutsCouple()->GetIndex();
00156   G4double tcut = (*(G4ProductionCutsTable::GetProductionCutsTable()
00157                      ->GetEnergyCutsVector(3)))[idx];
00158   hadi->SetRecoilEnergyThreshold(tcut);
00159 
00160   // Initialize the hadronic projectile from the track
00161   //  G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
00162   G4HadProjectile theProj(track);
00163   if(verboseLevel>1) {
00164     G4cout << "G4WHadronElasticProcess::PostStepDoIt for " 
00165            << part->GetParticleName()
00166            << " in " << material->GetName() 
00167            << " Target Z= " << targNucleus->GetZ_asInt() 
00168            << " A= " << targNucleus->GetA_asInt() << G4endl; 
00169   }
00170 
00171   G4HadFinalState* result = 0;
00172   try
00173     {
00174       result = hadi->ApplyYourself( theProj, *targNucleus);
00175     }
00176   catch(G4HadronicException aR)
00177     {
00178       G4ExceptionDescription ed;
00179       ed << "Call for " << hadi->GetModelName() << G4endl;
00180       ed << "Target element "<< elm->GetName()<<"  Z= " 
00181          << targNucleus->GetZ_asInt() 
00182          << "  A= " << targNucleus->GetA_asInt() << G4endl;
00183       DumpState(track,"ApplyYourself",ed);
00184       ed << " ApplyYourself failed" << G4endl;
00185       G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had006", 
00186                   FatalException, ed);
00187     }
00188 
00189   // Check the result for catastrophic energy non-conservation
00190   // cannot be applied because is not guranteed that recoil 
00191   // nucleus is created
00192   // result = CheckResult(theProj, targNucleus, result);
00193 
00194   // directions
00195   G4ThreeVector indir = track.GetMomentumDirection();
00196   G4double phi = CLHEP::twopi*G4UniformRand();
00197   G4ThreeVector it(0., 0., 1.);
00198   G4ThreeVector outdir = result->GetMomentumChange();
00199 
00200   if(verboseLevel>1) {
00201     G4cout << "Efin= " << result->GetEnergyChange()
00202            << " de= " << result->GetLocalEnergyDeposit()
00203            << " nsec= " << result->GetNumberOfSecondaries()
00204            << " dir= " << outdir
00205            << G4endl;
00206   }
00207 
00208   // energies  
00209   G4double edep = result->GetLocalEnergyDeposit();
00210   G4double efinal = result->GetEnergyChange();
00211   if(efinal < 0.0) { efinal = 0.0; }
00212   if(edep < 0.0)   { edep = 0.0; }
00213 
00214   // NOTE:  Very low energy scatters were causing numerical (FPE) errors
00215   //        in earlier releases; these limits have not been changed since.
00216   if((part == theNeutron && efinal <= lowestEnergyNeutron) || 
00217      (part != theNeutron && efinal <= lowestEnergy)) {
00218     edep += efinal;
00219     efinal = 0.0;
00220   }
00221 
00222   // primary change
00223   theTotalResult->ProposeEnergy(efinal);
00224 
00225   G4TrackStatus status = track.GetTrackStatus();
00226   if(efinal > 0.0) {
00227     outdir.rotate(phi, it);
00228     outdir.rotateUz(indir);
00229     theTotalResult->ProposeMomentumDirection(outdir);
00230   } else {
00231     if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)
00232          { status = fStopButAlive; }
00233     else { status = fStopAndKill; }
00234     theTotalResult->ProposeTrackStatus(status);
00235   }
00236 
00237   //G4cout << "Efinal= " << efinal << "  TrackStatus= " << status << G4endl;
00238 
00239   theTotalResult->SetNumberOfSecondaries(0);
00240 
00241   // recoil
00242   if(result->GetNumberOfSecondaries() > 0) {
00243     G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
00244 
00245     if(p->GetKineticEnergy() > lowestEnergy) {
00246       theTotalResult->SetNumberOfSecondaries(1);
00247       G4ThreeVector pdir = p->GetMomentumDirection();
00248       // G4cout << "recoil " << pdir << G4endl;
00250       pdir.rotate(phi, it);
00251       pdir.rotateUz(indir);
00252       // G4cout << "recoil rotated " << pdir << G4endl;
00253       p->SetMomentumDirection(pdir);
00254 
00255       // in elastic scattering time and weight are not changed
00256       G4Track* t = new G4Track(p, track.GetGlobalTime(), 
00257                                track.GetPosition());
00258       t->SetWeight(weight);
00259       t->SetTouchableHandle(track.GetTouchableHandle());
00260       theTotalResult->AddSecondary(t);
00261 
00262     } else {
00263       edep += p->GetKineticEnergy();
00264       delete p;
00265     }
00266   }
00267   theTotalResult->ProposeLocalEnergyDeposit(edep);
00268   theTotalResult->ProposeNonIonizingEnergyDeposit(edep);
00269   result->Clear();
00270 
00271   return theTotalResult;
00272 }