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 /// \file PrimaryGeneratorAction.cc
0027 /// \brief Implementation of the B5::PrimaryGeneratorAction class
0028
0029 #include "PrimaryGeneratorAction.hh"
0030
0031 #include "G4GenericMessenger.hh"
0032 #include "G4ParticleDefinition.hh"
0033 #include "G4ParticleGun.hh"
0034 #include "G4ParticleTable.hh"
0035 #include "G4SystemOfUnits.hh"
0036 #include "G4ThreeVector.hh"
0037 #include "Randomize.hh"
0038
0039 namespace B5
0040 {
0041
0042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0043
0044 PrimaryGeneratorAction::PrimaryGeneratorAction()
0045 {
0046   G4int nofParticles = 1;
0047   fParticleGun = new G4ParticleGun(nofParticles);
0048
0049   auto particleTable = G4ParticleTable::GetParticleTable();
0050   fPositron = particleTable->FindParticle("e+");
0051   fMuon = particleTable->FindParticle("mu+");
0052   fPion = particleTable->FindParticle("pi+");
0053   fKaon = particleTable->FindParticle("kaon+");
0054   fProton = particleTable->FindParticle("proton");
0055
0056   // default particle kinematics
0057   fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., -8. * m));
0058   fParticleGun->SetParticleDefinition(fPositron);
0059
0060   // define commands for this class
0061   DefineCommands();
0062 }
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065
0066 PrimaryGeneratorAction::~PrimaryGeneratorAction()
0067 {
0068   delete fParticleGun;
0069   delete fMessenger;
0070 }
0071
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073
0074 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
0075 {
0076   G4ParticleDefinition* particle;
0077   if (fRandomizePrimary) {
0078     auto i = (int)(5. * G4UniformRand());
0079     switch (i) {
0080       case 0:
0081         particle = fPositron;
0082         break;
0083       case 1:
0084         particle = fMuon;
0085         break;
0086       case 2:
0087         particle = fPion;
0088         break;
0089       case 3:
0090         particle = fKaon;
0091         break;
0092       default:
0093         particle = fProton;
0094         break;
0095     }
0096     fParticleGun->SetParticleDefinition(particle);
0097   }
0098   else {
0099     particle = fParticleGun->GetParticleDefinition();
0100   }
0101
0102   auto pp = fMomentum + (G4UniformRand() - 0.5) * fSigmaMomentum;
0103   auto mass = particle->GetPDGMass();
0104   auto ekin = std::sqrt(pp * pp + mass * mass) - mass;
0105   fParticleGun->SetParticleEnergy(ekin);
0106
0107   auto angle = (G4UniformRand() - 0.5) * fSigmaAngle;
0108   fParticleGun->SetParticleMomentumDirection(G4ThreeVector(std::sin(angle), 0., std::cos(angle)));
0109
0110   fParticleGun->GeneratePrimaryVertex(event);
0111 }
0112
0113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0114
0115 void PrimaryGeneratorAction::DefineCommands()
0116 {
0117   // Define /B5/generator command directory using generic messenger class
0118   fMessenger = new G4GenericMessenger(this, "/B5/generator/", "Primary generator control");
0119
0120   // momentum command
0121   auto& momentumCmd = fMessenger->DeclarePropertyWithUnit("momentum", "GeV", fMomentum,
0122                                                           "Mean momentum of primaries.");
0123   momentumCmd.SetParameterName("p", true);
0124   momentumCmd.SetRange("p>=0.");
0125   momentumCmd.SetDefaultValue("1.");
0126   // ok
0127   // momentumCmd.SetParameterName("p", true);
0128   // momentumCmd.SetRange("p>=0.");
0129
0130   // sigmaMomentum command
0131   auto& sigmaMomentumCmd = fMessenger->DeclarePropertyWithUnit(
0132     "sigmaMomentum", "MeV", fSigmaMomentum, "Sigma momentum of primaries.");
0133   sigmaMomentumCmd.SetParameterName("sp", true);
0134   sigmaMomentumCmd.SetRange("sp>=0.");
0135   sigmaMomentumCmd.SetDefaultValue("50.");
0136
0137   // sigmaAngle command
0138   auto& sigmaAngleCmd = fMessenger->DeclarePropertyWithUnit("sigmaAngle", "deg", fSigmaAngle,
0139                                                             "Sigma angle divergence of primaries.");
0140   sigmaAngleCmd.SetParameterName("t", true);
0141   sigmaAngleCmd.SetRange("t>=0.");
0142   sigmaAngleCmd.SetDefaultValue("2.");
0143
0144   // randomizePrimary command
0145   auto& randomCmd = fMessenger->DeclareProperty("randomizePrimary", fRandomizePrimary);
0146   G4String guidance = "Boolean flag for randomizing primary particle types.\n";
0147   guidance += "In case this flag is false, you can select the primary particle\n";
0148   guidance += "  with /gun/particle command.";
0149   randomCmd.SetGuidance(guidance);
0150   randomCmd.SetParameterName("flg", true);
0151   randomCmd.SetDefaultValue("true");
0152 }
0153
0154 //..oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0155
0156 }  // namespace B5