dd/d33/G4HadronElasticProcess_8cc_source.html

//

// ********************************************************************

// * License and Disclaimer                                           *

// *                                                                  *

// * The  Geant4 software  is  copyright of the Copyright Holders  of *

// * the Geant4 Collaboration.  It is provided  under  the terms  and *

// * conditions of the Geant4 Software License,  included in the file *

// * LICENSE and available at  http://cern.ch/geant4/license .  These *

// * include a list of copyright holders.                             *

// *                                                                  *

// * Neither the authors of this software system, nor their employing *

// * institutes,nor the agencies providing financial support for this *

// * work  make  any representation or  warranty, express or implied, *

// * regarding  this  software system or assume any liability for its *

// * use.  Please see the license in the file  LICENSE  and URL above *

// * for the full disclaimer and the limitation of liability.         *

// *                                                                  *

// * This  code  implementation is the result of  the  scientific and *

// * technical work of the GEANT4 collaboration.                      *

// * By using,  copying,  modifying or  distributing the software (or *

// * any work based  on the software)  you  agree  to acknowledge its *

// * use  in  resulting  scientific  publications,  and indicate your *

// * acceptance of all terms of the Geant4 Software license.          *

// ********************************************************************

//

// $Id$

//

// Geant4 Hadron Elastic Scattering Process

//

// Created 26 July 2012 V.Ivanchenko from G4WHadronElasticProcess

//

// Modified:

// 14-Sep-12 M.Kelsey -- Pass subType code to base ctor


#include <iostream>

#include <typeinfo>


#include "G4HadronElasticProcess.hh"

#include "G4SystemOfUnits.hh"

#include "G4Nucleus.hh"

#include "G4ProcessManager.hh"

#include "G4CrossSectionDataStore.hh"

#include "G4HadronElasticDataSet.hh"

#include "G4ProductionCutsTable.hh"

#include "G4HadronicException.hh"

#include "G4HadronicDeprecate.hh"


G4HadronElasticProcess::G4HadronElasticProcess(const G4String& pName)

  : G4HadronicProcess(pName, fHadronElastic), isInitialised(false)

{

  AddDataSet(new G4HadronElasticDataSet);

  lowestEnergy = 1.*keV;

}


G4HadronElasticProcess::~G4HadronElasticProcess()

{}


void G4HadronElasticProcess::Description() const

{

  char* dirName = getenv("G4PhysListDocDir");

  if (dirName) {

    std::ofstream outFile;

    G4String outFileName = GetProcessName() + ".html";

    G4String pathName = G4String(dirName) + "/" + outFileName;

    outFile.open(pathName);

    outFile << "<html>\n";

    outFile << "<head>\n";


    outFile << "<title>Description of G4HadronElasticProcess</title>\n";

    outFile << "</head>\n";

    outFile << "<body>\n";


    outFile << "G4HadronElasticProcess handles the elastic scattering of\n"

            << "hadrons by invoking one or more hadronic models and one or\n"

            << "more hadronic cross sections.\n";


    outFile << "</body>\n";

    outFile << "</html>\n";

    outFile.close();

  }

}


G4VParticleChange*

G4HadronElasticProcess::PostStepDoIt(const G4Track& track,

                     const G4Step& /*step*/)

{

  theTotalResult->Clear();

  theTotalResult->Initialize(track);

  G4double weight = track.GetWeight();

  theTotalResult->ProposeWeight(weight);


  // For elastic scattering, _any_ result is considered an interaction

  ClearNumberOfInteractionLengthLeft();


  G4double kineticEnergy = track.GetKineticEnergy();

  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();

  const G4ParticleDefinition* part = dynParticle->GetDefinition();


  // NOTE:  Very low energy scatters were causing numerical (FPE) errors

  //        in earlier releases; these limits have not been changed since.

  if (kineticEnergy <= lowestEnergy)   return theTotalResult;


  G4Material* material = track.GetMaterial();

  G4Nucleus* targNucleus = GetTargetNucleusPointer();


  // Select element

  G4Element* elm = 0;

  try

    {

      elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,

                            *targNucleus);

    }

  catch(G4HadronicException & aR)

    {

      G4ExceptionDescription ed;

      DumpState(track,"SampleZandA",ed);

      ed << " PostStepDoIt failed on element selection" << G4endl;

      G4Exception("G4HadronElasticProcess::PostStepDoIt", "had003",

          FatalException, ed);

    }

  G4HadronicInteraction* hadi = 0;

  try

    {

      hadi = ChooseHadronicInteraction( kineticEnergy, material, elm );

    }

  catch(G4HadronicException & aE)

    {

      G4ExceptionDescription ed;

      ed << "Target element "<< elm->GetName()<<"  Z= "

     << targNucleus->GetZ_asInt() << "  A= "

     << targNucleus->GetA_asInt() << G4endl;

      DumpState(track,"ChooseHadronicInteraction",ed);

      ed << " No HadronicInteraction found out" << G4endl;

      G4Exception("G4HadronElasticProcess::PostStepDoIt", "had005",

          FatalException, ed);

    }


  size_t idx = track.GetMaterialCutsCouple()->GetIndex();

  G4double tcut = (*(G4ProductionCutsTable::GetProductionCutsTable()

             ->GetEnergyCutsVector(3)))[idx];

  hadi->SetRecoilEnergyThreshold(tcut);


  // Initialize the hadronic projectile from the track

  //  G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;

  G4HadProjectile theProj(track);

  if(verboseLevel>1) {

    G4cout << "G4HadronElasticProcess::PostStepDoIt for "

       << part->GetParticleName()

       << " in " << material->GetName()

       << " Target Z= " << targNucleus->GetZ_asInt()

       << " A= " << targNucleus->GetA_asInt() << G4endl;

  }


  G4HadFinalState* result = 0;

  try

    {

      result = hadi->ApplyYourself( theProj, *targNucleus);

    }

  catch(G4HadronicException aR)

    {

      G4ExceptionDescription ed;

      ed << "Call for " << hadi->GetModelName() << G4endl;

      ed << "Target element "<< elm->GetName()<<"  Z= "

     << targNucleus->GetZ_asInt()

     << "  A= " << targNucleus->GetA_asInt() << G4endl;

      DumpState(track,"ApplyYourself",ed);

      ed << " ApplyYourself failed" << G4endl;

      G4Exception("G4HadronElasticProcess::PostStepDoIt", "had006",

          FatalException, ed);

    }


  // Check the result for catastrophic energy non-conservation

  // cannot be applied because is not guranteed that recoil

  // nucleus is created

  // result = CheckResult(theProj, targNucleus, result);


  // directions

  G4ThreeVector indir = track.GetMomentumDirection();

  G4double phi = CLHEP::twopi*G4UniformRand();

  G4ThreeVector it(0., 0., 1.);

  G4ThreeVector outdir = result->GetMomentumChange();


  if(verboseLevel>1) {

    G4cout << "Efin= " << result->GetEnergyChange()

       << " de= " << result->GetLocalEnergyDeposit()

       << " nsec= " << result->GetNumberOfSecondaries()

       << " dir= " << outdir

       << G4endl;

  }


  // energies

  G4double edep = result->GetLocalEnergyDeposit();

  G4double efinal = result->GetEnergyChange();

  if(efinal < 0.0) { efinal = 0.0; }

  if(edep < 0.0)   { edep = 0.0; }


  // NOTE:  Very low energy scatters were causing numerical (FPE) errors

  //        in earlier releases; these limits have not been changed since.

  if(efinal <= lowestEnergy) {

    edep += efinal;

    efinal = 0.0;

  }


  // primary change

  theTotalResult->ProposeEnergy(efinal);


  G4TrackStatus status = track.GetTrackStatus();

  if(efinal > 0.0) {

    outdir.rotate(phi, it);

    outdir.rotateUz(indir);

    theTotalResult->ProposeMomentumDirection(outdir);

  } else {

    if(part->GetProcessManager()->GetAtRestProcessVector()->size() > 0)

         { status = fStopButAlive; }

    else { status = fStopAndKill; }

    theTotalResult->ProposeTrackStatus(status);

  }


  //G4cout << "Efinal= " << efinal << "  TrackStatus= " << status << G4endl;


  theTotalResult->SetNumberOfSecondaries(0);


  // recoil

  if(result->GetNumberOfSecondaries() > 0) {

    G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();


    if(p->GetKineticEnergy() > tcut) {

      theTotalResult->SetNumberOfSecondaries(1);

      G4ThreeVector pdir = p->GetMomentumDirection();

      // G4cout << "recoil " << pdir << G4endl;

      //!! is not needed for models inheriting G4HadronElastic

      pdir.rotate(phi, it);

      pdir.rotateUz(indir);

      // G4cout << "recoil rotated " << pdir << G4endl;

      p->SetMomentumDirection(pdir);


      // in elastic scattering time and weight are not changed

      G4Track* t = new G4Track(p, track.GetGlobalTime(),

                   track.GetPosition());

      t->SetWeight(weight);

      t->SetTouchableHandle(track.GetTouchableHandle());

      theTotalResult->AddSecondary(t);


    } else {

      edep += p->GetKineticEnergy();

      delete p;

    }

  }

  theTotalResult->ProposeLocalEnergyDeposit(edep);

  theTotalResult->ProposeNonIonizingEnergyDeposit(edep);

  result->Clear();


  return theTotalResult;

}


void

G4HadronElasticProcess::PreparePhysicsTable(const G4ParticleDefinition& part)

{

  if(!isInitialised) {

    isInitialised = true;

    if(G4Neutron::Neutron() == &part) { lowestEnergy = 1.e-6*eV; }

  }

  G4HadronicProcess::PreparePhysicsTable(part);

}


void G4HadronElasticProcess::SetLowestEnergy(G4double val)

{

  lowestEnergy = val;

}


void

G4HadronElasticProcess::SetLowestEnergyNeutron(G4double val)

{

  lowestEnergy = val;

  G4HadronicDeprecate("G4HadronElasticProcess::SetLowestEnergyNeutron()");

}


G4CrossSectionDataStore.hh

FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:60

G4HadronElasticDataSet.hh

G4HadronElasticProcess.hh

G4HadronicDeprecate.hh

G4HadronicDeprecate
#define G4HadronicDeprecate(name)
Definition: G4HadronicDeprecate.hh:33

G4HadronicException.hh

fHadronElastic
@ fHadronElastic
Definition: G4HadronicProcessType.hh:44

G4Nucleus.hh

G4ProcessManager.hh

G4ProductionCutsTable.hh

G4SystemOfUnits.hh

G4TrackStatus
G4TrackStatus
Definition: G4TrackStatus.hh:51

fStopAndKill
@ fStopAndKill
Definition: G4TrackStatus.hh:56

fStopButAlive
@ fStopButAlive
Definition: G4TrackStatus.hh:54

G4double
double G4double
Definition: G4Types.hh:64

G4endl
#define G4endl
Definition: G4ios.hh:52

G4cout
G4DLLIMPORT std::ostream G4cout

G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53

CLHEP::Hep3Vector
Definition: ThreeVector.h:41

CLHEP::Hep3Vector::rotateUz
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72

CLHEP::Hep3Vector::rotate
Hep3Vector & rotate(double, const Hep3Vector &)
Definition: ThreeVectorR.cc:28

G4CrossSectionDataStore::SampleZandA
G4Element * SampleZandA(const G4DynamicParticle *, const G4Material *, G4Nucleus &target)
Definition: G4CrossSectionDataStore.cc:212

G4DynamicParticle
Definition: G4DynamicParticle.hh:74

G4DynamicParticle::SetMomentumDirection
void SetMomentumDirection(const G4ThreeVector &aDirection)

G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const

G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const

G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const

G4Element
Definition: G4Element.hh:98

G4Element::GetName
const G4String & GetName() const
Definition: G4Element.hh:127

G4HadFinalState
Definition: G4HadFinalState.hh:46

G4HadFinalState::GetEnergyChange
G4double GetEnergyChange() const
Definition: G4HadFinalState.cc:52

G4HadFinalState::Clear
void Clear()
Definition: G4HadFinalState.cc:89

G4HadFinalState::GetLocalEnergyDeposit
G4double GetLocalEnergyDeposit() const
Definition: G4HadFinalState.cc:131

G4HadFinalState::GetNumberOfSecondaries
G4int GetNumberOfSecondaries() const
Definition: G4HadFinalState.cc:40

G4HadFinalState::GetMomentumChange
const G4ThreeVector & GetMomentumChange() const
Definition: G4HadFinalState.cc:68

G4HadFinalState::GetSecondary
G4HadSecondary * GetSecondary(size_t i)
Definition: G4HadFinalState.cc:109

G4HadProjectile
Definition: G4HadProjectile.hh:40

G4HadSecondary::GetParticle
G4DynamicParticle * GetParticle()
Definition: G4HadSecondary.hh:39

G4HadronElasticDataSet
Definition: G4HadronElasticDataSet.hh:45

G4HadronElasticProcess::~G4HadronElasticProcess
virtual ~G4HadronElasticProcess()
Definition: G4HadronElasticProcess.cc:55

G4HadronElasticProcess::Description
virtual void Description() const
Definition: G4HadronElasticProcess.cc:58

G4HadronElasticProcess::PreparePhysicsTable
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
Definition: G4HadronElasticProcess.cc:257

G4HadronElasticProcess::PostStepDoIt
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
Definition: G4HadronElasticProcess.cc:84

G4HadronElasticProcess::G4HadronElasticProcess
G4HadronElasticProcess(const G4String &procName="hadElastic")
Definition: G4HadronElasticProcess.cc:48

G4HadronElasticProcess::SetLowestEnergy
virtual void SetLowestEnergy(G4double)
Definition: G4HadronElasticProcess.cc:266

G4HadronElasticProcess::SetLowestEnergyNeutron
virtual void SetLowestEnergyNeutron(G4double)
Definition: G4HadronElasticProcess.cc:272

G4HadronicException
Definition: G4HadronicException.hh:34

G4HadronicInteraction
Definition: G4HadronicInteraction.hh:64

G4HadronicInteraction::ApplyYourself
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0

G4HadronicInteraction::GetModelName
const G4String & GetModelName() const
Definition: G4HadronicInteraction.hh:119

G4HadronicInteraction::SetRecoilEnergyThreshold
void SetRecoilEnergyThreshold(G4double val)
Definition: G4HadronicInteraction.hh:142

G4HadronicProcess
Definition: G4HadronicProcess.hh:70

G4HadronicProcess::GetTargetNucleusPointer
G4Nucleus * GetTargetNucleusPointer()
Definition: G4HadronicProcess.hh:149

G4HadronicProcess::theTotalResult
G4ParticleChange * theTotalResult
Definition: G4HadronicProcess.hh:222

G4HadronicProcess::AddDataSet
void AddDataSet(G4VCrossSectionDataSet *aDataSet)
Definition: G4HadronicProcess.hh:117

G4HadronicProcess::GetCrossSectionDataStore
G4CrossSectionDataStore * GetCrossSectionDataStore()
Definition: G4HadronicProcess.hh:170

G4HadronicProcess::ChooseHadronicInteraction
G4HadronicInteraction * ChooseHadronicInteraction(G4double kineticEnergy, G4Material *aMaterial, G4Element *anElement)
Definition: G4HadronicProcess.hh:142

G4HadronicProcess::DumpState
void DumpState(const G4Track &, const G4String &, G4ExceptionDescription &)
Definition: G4HadronicProcess.cc:823

G4HadronicProcess::PreparePhysicsTable
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
Definition: G4HadronicProcess.cc:158

G4MaterialCutsCouple::GetIndex
G4int GetIndex() const
Definition: G4MaterialCutsCouple.hh:111

G4Material
Definition: G4Material.hh:119

G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:177

G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104

G4Nucleus
Definition: G4Nucleus.hh:51

G4Nucleus::GetA_asInt
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109

G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115

G4ParticleChange::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChange.cc:218

G4ParticleChange::ProposeEnergy
void ProposeEnergy(G4double finalEnergy)

G4ParticleChange::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)

G4ParticleChange::Initialize
virtual void Initialize(const G4Track &)
Definition: G4ParticleChange.cc:228

G4ParticleDefinition
Definition: G4ParticleDefinition.hh:68

G4ParticleDefinition::GetProcessManager
G4ProcessManager * GetProcessManager() const

G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:115

G4ProcessManager::GetAtRestProcessVector
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const

G4ProcessVector::size
G4int size() const

G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
Definition: G4ProductionCutsTable.hh:250

G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63

G4Step
Definition: G4Step.hh:78

G4String
Definition: G4String.hh:105

G4Track
Definition: G4Track.hh:72

G4Track::GetTrackStatus
G4TrackStatus GetTrackStatus() const

G4Track::GetWeight
G4double GetWeight() const

G4Track::SetWeight
void SetWeight(G4double aValue)

G4Track::GetPosition
const G4ThreeVector & GetPosition() const

G4Track::SetTouchableHandle
void SetTouchableHandle(const G4TouchableHandle &apValue)

G4Track::GetGlobalTime
G4double GetGlobalTime() const

G4Track::GetMaterial
G4Material * GetMaterial() const

G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const

G4Track::GetTouchableHandle
const G4TouchableHandle & GetTouchableHandle() const

G4Track::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const

G4Track::GetKineticEnergy
G4double GetKineticEnergy() const

G4Track::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple() const

G4VParticleChange
Definition: G4VParticleChange.hh:95

G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)

G4VParticleChange::ProposeNonIonizingEnergyDeposit
void ProposeNonIonizingEnergyDeposit(G4double anEnergyPart)

G4VParticleChange::ProposeWeight
void ProposeWeight(G4double finalWeight)

G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)

G4VParticleChange::SetNumberOfSecondaries
void SetNumberOfSecondaries(G4int totSecondaries)

G4VParticleChange::Clear
void Clear()

G4VProcess::ClearNumberOfInteractionLengthLeft
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:418

G4VProcess::verboseLevel
G4int verboseLevel
Definition: G4VProcess.hh:368

G4VProcess::GetProcessName
const G4String & GetProcessName() const
Definition: G4VProcess.hh:379

G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41

G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76