#include "GenePrimaryGeneratorAction.hh" #include "TGraph.h" #include "TH1F.h" #include "G4Event.hh" #include "G4ParticleDefinition.hh" #include "G4ParticleGun.hh" #include "G4ParticleTable.hh" #include "G4PhysicalConstants.hh" #include "G4SystemOfUnits.hh" #include "GeneAnalysisManager.hh" #include "GenePrimaryGeneratorActionMessenger.hh" #include "Randomize.hh" #include "globals.hh" #include "fstream" #define ANALYSIS_USE 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); #ifdef ANALYSIS_USE GeneAnalysisManager* analysisManager = GeneAnalysisManager::GetInstance(); analysisManager->FillAng(thetacmr, thetalabr); #endif particleGun->SetParticleMomentumDirection(G4ThreeVector(x0, y0, z0)); particleGun->SetParticleEnergy(e4_lab * MeV); particleGun->GeneratePrimaryVertex(anEvent); }