B5/src/PrimaryGeneratorAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file B5/src/PrimaryGeneratorAction.cc
0028 /// \brief Implementation of the B5::PrimaryGeneratorAction class
0029
0030 #include "PrimaryGeneratorAction.hh"
0031
0032 #include "G4GenericMessenger.hh"
0033 #include "G4ParticleDefinition.hh"
0034 #include "G4ParticleGun.hh"
0035 #include "G4ParticleTable.hh"
0036 #include "G4SystemOfUnits.hh"
0037 #include "G4ThreeVector.hh"
0038 #include "Randomize.hh"
0039
0040 namespace B5
0041 {
0042
0043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0044
0045 PrimaryGeneratorAction::PrimaryGeneratorAction()
0046 {
0047   G4int nofParticles = 1;
0048   fParticleGun = new G4ParticleGun(nofParticles);
0049
0050   auto particleTable = G4ParticleTable::GetParticleTable();
0051   fPositron = particleTable->FindParticle("e+");
0052   fMuon = particleTable->FindParticle("mu+");
0053   fPion = particleTable->FindParticle("pi+");
0054   fKaon = particleTable->FindParticle("kaon+");
0055   fProton = particleTable->FindParticle("proton");
0056
0057   // default particle kinematics
0058   fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., -8. * m));
0059   fParticleGun->SetParticleDefinition(fPositron);
0060
0061   // define commands for this class
0062   DefineCommands();
0063 }
0064
0065 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0066
0067 PrimaryGeneratorAction::~PrimaryGeneratorAction()
0068 {
0069   delete fParticleGun;
0070   delete fMessenger;
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
0076 {
0077   G4ParticleDefinition* particle;
0078   if (fRandomizePrimary) {
0079     auto i = (int)(5. * G4UniformRand());
0080     switch (i) {
0081       case 0:
0082         particle = fPositron;
0083         break;
0084       case 1:
0085         particle = fMuon;
0086         break;
0087       case 2:
0088         particle = fPion;
0089         break;
0090       case 3:
0091         particle = fKaon;
0092         break;
0093       default:
0094         particle = fProton;
0095         break;
0096     }
0097     fParticleGun->SetParticleDefinition(particle);
0098   }
0099   else {
0100     particle = fParticleGun->GetParticleDefinition();
0101   }
0102
0103   auto pp = fMomentum + (G4UniformRand() - 0.5) * fSigmaMomentum;
0104   auto mass = particle->GetPDGMass();
0105   auto ekin = std::sqrt(pp * pp + mass * mass) - mass;
0106   fParticleGun->SetParticleEnergy(ekin);
0107
0108   auto angle = (G4UniformRand() - 0.5) * fSigmaAngle;
0109   fParticleGun->SetParticleMomentumDirection(G4ThreeVector(std::sin(angle), 0., std::cos(angle)));
0110
0111   fParticleGun->GeneratePrimaryVertex(event);
0112 }
0113
0114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0115
0116 void PrimaryGeneratorAction::DefineCommands()
0117 {
0118   // Define /B5/generator command directory using generic messenger class
0119   fMessenger = new G4GenericMessenger(this, "/B5/generator/", "Primary generator control");
0120
0121   // momentum command
0122   auto& momentumCmd = fMessenger->DeclarePropertyWithUnit("momentum", "GeV", fMomentum,
0123                                                           "Mean momentum of primaries.");
0124   momentumCmd.SetParameterName("p", true);
0125   momentumCmd.SetRange("p>=0.");
0126   momentumCmd.SetDefaultValue("1.");
0127   // ok
0128   // momentumCmd.SetParameterName("p", true);
0129   // momentumCmd.SetRange("p>=0.");
0130
0131   // sigmaMomentum command
0132   auto& sigmaMomentumCmd = fMessenger->DeclarePropertyWithUnit(
0133     "sigmaMomentum", "MeV", fSigmaMomentum, "Sigma momentum of primaries.");
0134   sigmaMomentumCmd.SetParameterName("sp", true);
0135   sigmaMomentumCmd.SetRange("sp>=0.");
0136   sigmaMomentumCmd.SetDefaultValue("50.");
0137
0138   // sigmaAngle command
0139   auto& sigmaAngleCmd = fMessenger->DeclarePropertyWithUnit("sigmaAngle", "deg", fSigmaAngle,
0140                                                             "Sigma angle divergence of primaries.");
0141   sigmaAngleCmd.SetParameterName("t", true);
0142   sigmaAngleCmd.SetRange("t>=0.");
0143   sigmaAngleCmd.SetDefaultValue("2.");
0144
0145   // randomizePrimary command
0146   auto& randomCmd = fMessenger->DeclareProperty("randomizePrimary", fRandomizePrimary);
0147   G4String guidance = "Boolean flag for randomizing primary particle types.\n";
0148   guidance += "In case this flag is false, you can select the primary particle\n";
0149   guidance += "  with /gun/particle command.";
0150   randomCmd.SetGuidance(guidance);
0151   randomCmd.SetParameterName("flg", true);
0152   randomCmd.SetDefaultValue("true");
0153 }
0154
0155 //..oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0156
0157 }  // namespace B5