204 lines
6.9 KiB
C++
204 lines
6.9 KiB
C++
#include "GenePrimaryGeneratorAction.hh"
|
|
#include "TGraph.h"
|
|
#include "TH1F.h"
|
|
|
|
#include "GeneAnalysisManager.hh"
|
|
#include "GenePrimaryGeneratorActionMessenger.hh"
|
|
#include "Randomize.hh"
|
|
#include <G4Event.hh>
|
|
#include <G4ParticleDefinition.hh>
|
|
#include <G4ParticleGun.hh>
|
|
#include <G4ParticleTable.hh>
|
|
#include <G4PhysicalConstants.hh>
|
|
#include <G4SystemOfUnits.hh>
|
|
#include <globals.hh>
|
|
|
|
#include "fstream"
|
|
|
|
GenePrimaryGeneratorAction::GenePrimaryGeneratorAction() {
|
|
G4int n_particle = 1;
|
|
particleGun = new G4ParticleGun(n_particle);
|
|
|
|
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
|
|
G4String particleName;
|
|
// 这里模拟的是反应生成的中子
|
|
particleGun->SetParticleDefinition(particleTable->FindParticle(particleName = "neutron"));
|
|
particleGun->SetParticlePosition(G4ThreeVector());
|
|
|
|
// 通过角分布采样获取随机数
|
|
hAng = new TH1F("hAng", "hAng;theta(deg);prob.", 180, 0, 180);
|
|
|
|
fBeamEnergy = 0;
|
|
pMessenger = new GenePrimaryGeneratorActionMessenger(this);
|
|
}
|
|
|
|
GenePrimaryGeneratorAction::~GenePrimaryGeneratorAction() {
|
|
delete particleGun;
|
|
delete hAng;
|
|
delete pMessenger;
|
|
}
|
|
|
|
void GenePrimaryGeneratorAction::SetReaction(G4String& name) {
|
|
if (name == "V51pn")
|
|
iReaction = 0;
|
|
else if (name == "C13an")
|
|
iReaction = 1;
|
|
else {
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cin >> iReaction;
|
|
}
|
|
}
|
|
|
|
// ???
|
|
void GenePrimaryGeneratorAction::SetBeamEnergy(G4double val) {
|
|
fBeamEnergy = val / MeV;
|
|
|
|
G4double umass = 931.494;
|
|
G4double massHe4 = 4.00260325413 * umass;
|
|
G4double massC13 = 13.00335483521 * umass;
|
|
G4double massO16 = 15.9949146196 * umass;
|
|
G4double massn = 939.56563;
|
|
G4double massp = 938.78294;
|
|
G4double massV51 = 50.9439569 * umass;
|
|
G4double massCr51 = 50.9447647 * umass;
|
|
|
|
if (iReaction == 0) {
|
|
mass[0] = massp;
|
|
mass[1] = massV51;
|
|
mass[2] = massCr51;
|
|
mass[3] = massn;
|
|
} else if (iReaction == 1) {
|
|
mass[0] = massHe4;
|
|
mass[1] = massC13;
|
|
mass[2] = massO16;
|
|
mass[3] = massn;
|
|
} else {
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cout << "Unknown reaction type!!!!" << G4endl;
|
|
G4cin >> iReaction;
|
|
}
|
|
|
|
G4double eb = fBeamEnergy + mass[0];
|
|
G4double pb2 = fBeamEnergy * fBeamEnergy + 2.0 * fBeamEnergy * mass[0];
|
|
G4double pb = std::sqrt(pb2);
|
|
G4double esum = eb + mass[1];
|
|
beta = pb / esum;
|
|
gamma = 1.0 / std::sqrt(1.0 - beta * beta);
|
|
G4double e_cm2 = esum * esum - pb2;
|
|
G4double e_cm = std::sqrt(e_cm2);
|
|
G4double t_cm = e_cm - mass[2] - mass[3];
|
|
if (t_cm < 0.0) {
|
|
G4cout << "No solution!" << G4endl;
|
|
G4cin >> iReaction;
|
|
return;
|
|
}
|
|
G4double t_cm2 = t_cm * t_cm;
|
|
t3_cm = (t_cm2 + 2.0 * mass[3] * t_cm) / e_cm / 2.0;
|
|
t4_cm = (t_cm2 + 2.0 * mass[2] * t_cm) / e_cm / 2.0;
|
|
G4double p3_cm2 = t3_cm * t3_cm + 2.0 * t3_cm * mass[2];
|
|
G4double p4_cm2 = t4_cm * t4_cm + 2.0 * t4_cm * mass[3];
|
|
p3_cm = std::sqrt(p3_cm2);
|
|
p4_cm = std::sqrt(p4_cm2);
|
|
}
|
|
|
|
// ???
|
|
void GenePrimaryGeneratorAction::SetAngType(G4String& name) {
|
|
if (name == "ENDF") {
|
|
std::ifstream infile("xt_nov21.a");
|
|
G4int iN = 0;
|
|
if (infile.fail()) {
|
|
G4cout << "ifstream: Failure to open input file!!!" << G4endl;
|
|
G4cout << "ifstream: Failure to open input file!!!" << G4endl;
|
|
G4cout << "ifstream: Failure to open input file!!!" << G4endl;
|
|
G4cout << "ifstream: Failure to open input file!!!" << G4endl;
|
|
G4cout << "ifstream: Failure to open input file!!!" << G4endl;
|
|
G4cin >> iN;
|
|
return;
|
|
}
|
|
G4double Elab[5000], A1[5000], A2[5000], A3[5000];
|
|
G4double felab, fecm, fa0, fa1, fa2, fa3;
|
|
while (infile >> felab >> fecm >> fa0 >> fa1 >> fa2 >> fa3) {
|
|
Elab[iN] = felab;
|
|
A1[iN] = fa1;
|
|
A2[iN] = fa2;
|
|
A3[iN] = fa3;
|
|
iN++;
|
|
}
|
|
infile.close();
|
|
TGraph gr1(iN, Elab, A1);
|
|
TGraph gr2(iN, Elab, A2);
|
|
TGraph gr3(iN, Elab, A3);
|
|
fa1 = gr1.Eval(fBeamEnergy, 0, "S");
|
|
fa2 = gr2.Eval(fBeamEnergy, 0, "S");
|
|
fa3 = gr3.Eval(fBeamEnergy, 0, "S");
|
|
G4double thcm, costh, pthcm;
|
|
for (G4int i = 0; i < 180; i++) {
|
|
thcm = pi * i / 180.0;
|
|
costh = std::cos(thcm);
|
|
pthcm = 1.0 + fa1 * costh + fa2 * 0.5 * (3.0 * costh * costh - 1.0) +
|
|
fa3 * 0.5 * (5.0 * std::pow(costh, 3.0) - 3.0 * costh);
|
|
pthcm = pthcm * 2.0 * pi * std::sin(thcm);
|
|
if (pthcm < 0.0) pthcm = 0.0;
|
|
hAng->SetBinContent(i + 1, pthcm);
|
|
}
|
|
} else if (name == "ISO") {
|
|
G4double thcm, pthcm;
|
|
for (G4int i = 0; i < 180; i++) {
|
|
thcm = pi * i / 180.0;
|
|
pthcm = 2.0 * pi * std::sin(thcm);
|
|
hAng->SetBinContent(i + 1, pthcm);
|
|
}
|
|
} else {
|
|
G4cout << "Wrong Angular type!!!!" << G4endl;
|
|
G4cout << "Wrong Angular type!!!!" << G4endl;
|
|
G4cout << "Wrong Angular type!!!!" << G4endl;
|
|
G4cout << "Wrong Angular type!!!!" << G4endl;
|
|
G4cout << "Wrong Angular type!!!!" << G4endl;
|
|
G4int iN;
|
|
G4cin >> iN;
|
|
return;
|
|
}
|
|
}
|
|
|
|
void GenePrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {
|
|
G4double thetacmr, thetalabr;
|
|
thetacmr = hAng->GetRandom() / 180.0 * pi;
|
|
|
|
G4double tg_thetalabr =
|
|
p4_cm * std::sin(thetacmr) /
|
|
(gamma * (p4_cm * std::cos(thetacmr) + beta * std::sqrt(p4_cm * p4_cm + mass[3] * mass[3])));
|
|
if (tg_thetalabr >= 1.0e6)
|
|
thetalabr = pi / 2.0;
|
|
else
|
|
thetalabr = std::atan(tg_thetalabr);
|
|
if (thetalabr < 0.0) thetalabr = pi + thetalabr;
|
|
|
|
G4double p4_cmx = p4_cm * std::sin(thetacmr);
|
|
G4double p4_cmz = p4_cm * std::cos(thetacmr);
|
|
G4double p4_labx = p4_cmx;
|
|
G4double p4_labz = gamma * (p4_cmz + beta * (t4_cm + mass[3]));
|
|
G4double p4_lab = std::sqrt(p4_labx * p4_labx + p4_labz * p4_labz);
|
|
G4double e4_lab = sqrt(p4_lab * p4_lab + mass[3] * mass[3]) - mass[3];
|
|
|
|
G4double phi = (G4UniformRand() * 2 * pi) - pi;
|
|
G4double x0 = std::sin(thetalabr) * std::cos(phi);
|
|
G4double y0 = std::sin(thetalabr) * std::sin(phi);
|
|
G4double z0 = std::cos(thetalabr);
|
|
|
|
GeneAnalysisManager* analysisManager = GeneAnalysisManager::GetInstance();
|
|
analysisManager->FillAng(thetacmr, thetalabr);
|
|
|
|
particleGun->SetParticleMomentumDirection(G4ThreeVector(x0, y0, z0));
|
|
|
|
particleGun->SetParticleEnergy(e4_lab * MeV);
|
|
|
|
particleGun->GeneratePrimaryVertex(anEvent);
|
|
}
|