#include "CloverPrimaryGeneratorAction.hh" #include "CloverEventAction.hh" #include "G4RunManager.hh" #include "G4LogicalVolumeStore.hh" #include "G4LogicalVolume.hh" #include "G4Box.hh" #include "G4Event.hh" #include "G4ParticleGun.hh" #include "G4ParticleTable.hh" #include "G4ParticleDefinition.hh" #include "G4SystemOfUnits.hh" #include "G4GenericMessenger.hh" #include "Randomize.hh" #include //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... CloverPrimaryGeneratorAction::CloverPrimaryGeneratorAction() : G4VUserPrimaryGeneratorAction(), fMessenger(nullptr), fParticleGun(nullptr), fAngle(180 * degree) //,fEventAction(eventAction) { G4int nofParticles = 1; fParticleGun = new G4ParticleGun(nofParticles); // default particle kinematic // auto particleDefinition = G4ParticleTable::GetParticleTable()->FindParticle("gamma"); fParticleGun->SetParticleDefinition(particleDefinition); fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.)); fParticleGun->SetParticleEnergy(6.13*MeV); fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 0)); //Cs137 //fEnergyList = {{.663,100}}; //Co60 fEnergyList = {{1.173,99.973}, {1.332,99.986}}; //N16 iso /* fEnergyList = {{0.120, 100}, //energy, branch {1.775, 0.114}, {1.955, 0.036}, {2.742, 0.777}, {8.872, 0.072}, {7.117, 5}, {6.917, 0.036}, {6.130, 68.0+0.777}}; */ //Eu152 /* fEnergyList = {{.12178,28.58}, {.14801,0.037}, {.21257,0.019}, {.24469,7.583}, {.29594,0.477}, {.32913,0.128}, {.34428,26.50}, {.36778,0.861}, {.41112,2.234}, {.44396,0.327}, {.44397,2.821}, {.48868,0.419}, {.50347,0.149}, {.56399,0.489}, {.56643,0.129}, {.58627,0.459}, {.65648,0.145}, {.67468,0.172}, {.67862,0.471}, {.68867,0.857}, {.71935,0.321}, {.76490,0.215}, {.77891,12.942}, {.81045,0.320}, {.86738,4.245}, {.91933,0.427}, {.92632,0.278}, {.96339,0.135}, {.96408,14.605}, {1.00527,0.646}, {1.0840,0.246}, {1.08586,10.207}, {1.089737,1.727}, {1.109174,0.186}, {1.11207,13.644}, {1.21294,1.422}, {1.249938,0.188}, {1.292778,0.105}, {1.29914,1.623}, {1.408006,21.005}, {1.457643,0.502}, {1.528103,0.281}}; */ //accumulate the branch G4int numEnergy = (G4int) fEnergyList.size(); for( int i = 1 ; i < numEnergy ; i++){ fEnergyList[i][1] += fEnergyList[i-1][1]; } //renormalized the branching ratio for( int i = 0 ; i < numEnergy ; i++){ fEnergyList[i][1] = fEnergyList[i][1]/fEnergyList[numEnergy-1][1]; } //Define command for this class DefineCommands(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... CloverPrimaryGeneratorAction::~CloverPrimaryGeneratorAction() { delete fParticleGun; delete fMessenger; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void CloverPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) { G4double cosTheta = 1.- (1.-cos(fAngle ))*G4UniformRand(); G4double phi = 360*degree*G4UniformRand(); G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta); G4double ux = sinTheta*std::cos(phi), uy = sinTheta*std::sin(phi), uz = cosTheta; fParticleGun->SetParticleMomentumDirection(G4ThreeVector(ux,uy,uz)); fParticleGun->GeneratePrimaryVertex(anEvent); //generate a random number form 0 to 1 G4double pop = G4UniformRand(); G4double beamEnergy = 6.13; G4int numEnergy = (G4int) fEnergyList.size(); for( int i = 0; i < numEnergy ; i++){ if( i == 0 && pop <= fEnergyList[i][1] ) beamEnergy = fEnergyList[0][0] ; if( i > 0 && fEnergyList[i-1][1] < pop && pop <= fEnergyList[i][1] ) beamEnergy = fEnergyList[i][0] ; } fParticleGun->SetParticleEnergy(beamEnergy); //G4double energy = fParticleGun->GetParticleEnergy(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void CloverPrimaryGeneratorAction::DefineCommands() { // Define /Clover/generator command directory using generic messenger class fMessenger = new G4GenericMessenger(this, "/Clover/generator/", "Primary generator control"); // OpeningAngle command auto& openingAngleCmd = fMessenger->DeclarePropertyWithUnit("OpeningAngle", "deg", fAngle, "Opening Angle of beam in deg"); openingAngleCmd.SetParameterName("t", true); openingAngleCmd.SetRange("t>=0."); openingAngleCmd.SetDefaultValue("180."); }