Hi all, I want to output the double differential cross section for neutron captures in a target material. I guess that GEANT4 calculates these, but probably does not have a spline table that can just be printed containing this information.
If I wanted to calculate the double differential cross section, what is the simplest method to extract this? I tried to use G4HadronicProcessStore just as in Hadr00 but this only gives the cross section per atom and not say per energy and scattering angle.
I can see an example of what I want to achieve here:
Thanks for your suggestions!
_Geant4 Version:_11.0.0
_Operating System:_MacOS
_Compiler/Version:_Apple clang version 14.0.0
_CMake Version:_3.25.1
So in my RunAction I put the following and it seems to be outputting some cross section information, albeit the single cross section and not double differential cross section:
#include "RunAction.hh"
#include "Run.hh"
#include "DetectorConstruction.hh"
#include "PrimaryGeneratorAction.hh"
#include "HistoManager.hh"
#include "G4HadronicProcessStore.hh"
#include "G4Run.hh"
#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"
#include "G4RunManager.hh"
#include "Randomize.hh"
#include <iomanip>
...
void RunAction::EndOfRunAction(const G4Run*)
{
if (isMaster) fRun->EndOfRun();
// Here we are going to define the cross section tables
//save histograms
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
G4RunManager* runManager = G4RunManager::GetRunManager();
G4cout<< "Run manager created!"<< G4endl;
G4double p1 = std::log10(1*MeV);
G4double p2 = std::log10(1*TeV);
G4double e1 = std::log10(1*MeV);
G4double e2 = std::log10(1*TeV);
G4double de = (e2-e1)/G4double(1000);
G4double dp = (p2-p1)/G4double(1000);
G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
/* G4double mass = analysisManager->; */
G4cout<< "e1= "<< e1 <<" e2= "<< e2 <<" p1= "<< p1 <<" p2= "<< p2 <<G4endl;
const G4Material* mat = fDetector->GetMaterial();
G4cout << "Detector Material = " << mat << G4endl;
//Now we want to grab the properties!
G4double x = e1 - de*0.5;
G4double e,p,xs,tot;
G4int i;
G4int Verbose =1;
PrimaryGeneratorAction* generatorAction = const_cast<PrimaryGeneratorAction*>(dynamic_cast<const PrimaryGeneratorAction*>(G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction()));
/* PrimaryGeneratorAction* generatorAction = static_cast<const PrimaryGeneratorAction*>(G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction()); */
/* const PrimaryGeneratorAction* generatorAction = G4RunManager::GetRunManager()->GetUserPrimaryGeneratorAction(); */
const G4ParticleGun* particleGun = generatorAction->GetParticleGun();
G4cout << "Particle Gun" << particleGun << G4endl;
G4cout<< "particle details!"<< G4endl;
G4ParticleDefinition* particleDefinition = particleGun->GetParticleDefinition();
G4String particleName = particleDefinition->GetParticleName();
G4double kineticEnergy = particleGun->GetParticleEnergy(); // in MeV
G4cout << "particle name:" << particleName << G4endl;
G4cout << "particle energy:" << kineticEnergy/GeV <<"GeV" << G4endl;
std::cin.get();
G4double coeff =1.0;
if(mat) {coeff = mat->GetDensity()*cm2/g;}
// Get the number of elements in the material
G4int numElements = mat->GetNumberOfElements();
// Loop through each element in the material
for (G4int i = 0; i < numElements; ++i) {
// Get the ith element in the material
const G4Element* elm = mat->GetElement(i);
// Extract information about the element
G4String elmName = elm->GetName();
G4double elmFraction = mat->GetFractionVector()[i];
// Print or use the information as needed
G4cout << "Element " << i << ": Name = " << elmName << ", Fraction = " << elmFraction << G4endl;
std::cin.get();
for (int j=0; j<1000; j++){
try{//try bracket
//x=0.;
x += de;
e = std::pow(10.,x)*MeV;
if(Verbose>0) G4cout << std::setw(5) << j << std::setw(12) << e;
xs = store->GetElasticCrossSectionPerAtom(particleDefinition,e,elm,mat);
tot = xs;
if(Verbose>0) G4cout << std::setw(12) << xs/barn;
analysisManager->FillH1(8, x, xs/barn);
xs = store->GetInelasticCrossSectionPerAtom(particleDefinition,e,elm,mat);
tot += xs;
throw(j);
}//try bracket
catch (int j) {//catch bracket
G4cout << "Failed for index" << j << G4endl;
}//catch bracket
}
}
which is outputting something like
Particle Gun0x600001738240
particle details!
particle name:tau-
particle energy:60.856GeV
Element 0: Name = IceHydrogen, Fraction = 0.111846
0 1.00693 0Failed for index0
1 1.02094 0Failed for index1
2 1.03514 0Failed for index2
3 1.04954 0Failed for index3
4 1.06414 0Failed for index4
5 1.07895 0Failed for index5
6 1.09396 0Failed for index6
7 1.10917 0Failed for index7
8 1.1246 0Failed for index8
9 1.14025 0Failed for index9
10 1.15611 0Failed for index10
based off the first part of the Hadr00 example. How should I proceed? Thanks!