Dear Colleagues,
I am trying to create an H2 with a differential cross-section for muon interaction processes in iron (in cm^{2}/g) on the y-axis vs. normalized energy transferred to secondaries on the x-axis. I am using as a reference example TestEm17, since it’s written in the README: “It allows to compute differential cross sections (as function of the energy transfered to secondaries), total cross sections and to compare with analytic calculations.”
I see that the following is printed in the output file:
Simulation: total CrossSection = 0.49418 /cm MeanFreePath = 2.0235 cm massicCrossSection = 0.062794 cm2/g
However, the massicCrossSection in cm^2/g printed here seems to be for the entire run, while I would like to have it for each ν = ε/E.
I have tried to define the H2 in HistoManager and to modify SteppingAction.cc of TestEm17 in the following manner (adding after fHistoManager->FillHisto(id,lgepsE);):
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
G4Material* material = fDetector->GetMaterial();
G4double length = fDetector->GetSize();
G4double density = material->GetDensity();
G4double massCrossSection;
G4double ekin = fPrimary->GetParticleGun()->GetParticleEnergy();
MuCrossSections crossSections;
G4int id2 = 0; G4double cut = 0.;
if (procName == “muIoni”) {id2 = 1; cut = GetEnergyCut(material,1);}
else if (procName == “muPairProd”) {id2 = 2; cut = 2*(GetEnergyCut(material,1)
+ electron_mass_c2); }
else if (procName == “muBrems”) {id2 = 3; cut = GetEnergyCut(material,0);}
//else if (procName == “muonNuclear”){id2 = 4; }
else if ((procName != “muIoni”) && (procName != “muPairProd”) && (procName != “muBrems”) && (procName != “hIoni”) && (procName != “hPairProd”) && (procName != “hBrems”)) id2 = 4;
G4int nbOfBinsX = 100;
G4double binMinX = -10.;
G4double binMaxX = 0.;
G4double binWidthX = (binMaxX-binMinX)/G4double(nbOfBinsX);
G4int nbOfBinsY = 100;
G4double binMinY = -10.;
G4double binMaxY = 0.;
G4double binWidthY = (binMaxY-binMinY)/G4double(nbOfBinsY);
//create histo for theoritical crossSections, with same bining as simulation
//
G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
G4H2* histoTh = 0;
// G4H2* histoMC = 0;
if (fHistoManager->HistoExist2(id2)) {
histoTh = analysisManager->GetH2(fHistoManager->GetHistoID2(id2));
nbOfBinsX = fHistoManager->GetNbinsX2(id2);
binMinX = fHistoManager->GetVminX2 (id2);
binMaxX = fHistoManager->GetVmaxX2 (id2);
binWidthX = fHistoManager->GetBinWidthX2(id2);
nbOfBinsY = fHistoManager->GetNbinsY2(id2);
binMinY = fHistoManager->GetVminY2 (id2);
binMaxY = fHistoManager->GetVmaxY2 (id2);
binWidthY = fHistoManager->GetBinWidthY2(id2);
}
//compute and plot differential crossSection, as function of energy transfert.
//compute and return integrated crossSection for a given process.
//(note: to compare with simulation, the integrated crossSection is function
// of the energy cut.)
//
G4double lgeps, etransf, sigmaE, dsigma;
G4double diffMassCrSec;
G4double sigmaTot = 0.;
const G4double ln10 = std::log(10.);
for (G4int ibin=0; ibin<nbOfBinsX; ibin++) {
lgeps = binMinX + (ibin+0.5)binWidthX;
etransf = ekinstd::pow(10.,lgeps);
sigmaE = crossSections.CR_Macroscopic(procName,material,ekin,etransf);
dsigma = sigmaEetransfbinWidthXln10;
if (etransf > cut) sigmaTot += dsigma;
if (histoTh) {
G4double NbProcess = 1length*dsigma;
diffMassCrSec = dsigma/density;
//histoTh->fill(lgeps, diffMassCrSec, NbProcess);
}
}
//return integrated crossSection
//
massCrossSection = sigmaTot/density;
G4double massCrossSectionE = 0.;
if (sigmaTot > 0.) massCrossSectionE = std::log10(massCrossSection);
fHistoManager->FillHisto2(id2,lgepsE,massCrossSectionE);
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
#include “G4ProductionCutsTable.hh”
G4double SteppingAction::GetEnergyCut(G4Material* material, G4int idParticle)
{
G4ProductionCutsTable* table = G4ProductionCutsTable::GetProductionCutsTable();
size_t index = 0;
while ( (table->GetMaterialCutsCouple(index)->GetMaterial() != material) &&
(index < table->GetTableSize())) index++;
return (*(table->GetEnergyCutsVector(idParticle)))[index];
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
Unfortunately, I fail to see the same plots as the ones obtained by other groups (e.g. Fig. 1 in “Geant4 Simulation of Production and Interaction of Muons” by A. G. Bogdanov et al. Do you have any suggestions on what could be going wrong?
Thank you very much in advance. And to those of you who celebrate, Happy Easter.