jetcounter/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 PrimaryGeneratorAction class
0028
0029 #include "PrimaryGeneratorAction.hh"
0030
0031 #include "G4Event.hh"
0032 #include "G4ParticleGun.hh"
0033 #include "G4PhysicalConstants.hh"
0034 #include "G4SystemOfUnits.hh"
0035 #include "Randomize.hh"
0036
0037 #include <fstream>
0038
0039 PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* detector)
0040   : G4VUserPrimaryGeneratorAction(), fDetector(detector)
0041 {
0042   fpParticleGun = std::make_unique<G4ParticleGun>();
0043   fpParticleGun->SetParticleEnergy(-1);  // default value - can be overridden in the macro file
0044 }
0045
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047
0048 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
0049 {
0050   // user can set the energy in the macro file
0051   G4double energy = fpParticleGun->GetParticleEnergy();
0052
0053   // only first check is important as only user can set energy to -1,
0054   // hopefully...
0055   if (energy == -1) {  // if energy is larger than zero, then the beam is
0056                        // mono-energetic if the energy is set to zero, then the
0057                        // energy is randomized, based on spectrum file
0058     fMonoEnergetic = false;
0059   }
0060
0061   if (!fMonoEnergetic) {
0062     energy = GenerateParticleEnergy();
0063   }
0064
0065   fpParticleGun->SetParticleEnergy(energy);
0066   fpParticleGun->SetParticlePosition(GenerateParticlePosition());  // point of emission
0067   fpParticleGun->SetParticleMomentumDirection(
0068     GenerateParticleDirection());  // direction of emission
0069
0070   fpParticleGun->GeneratePrimaryVertex(anEvent);  // sending the particle
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 G4double PrimaryGeneratorAction::GenerateParticleEnergy()
0076 {
0077   if (fEnergySpectrum_length == 0) {  // reading the spectrum file if not loaded yet
0078     std::ifstream fin(fEnergySpectrumFilename);
0079     G4String len, gain, offset, counts;
0080     fin >> len >> gain >> offset;
0081     fEnergySpectrum_length = std::stoi(len);
0082     fEnergySpectrum_gain = std::stod(gain);
0083     fEnergySpectrum_offset = std::stod(offset);
0084
0085     fEnergySpectrum_counts.resize(fEnergySpectrum_length);
0086     for (G4int i = 0; i < fEnergySpectrum_length; i++) {
0087       fin >> counts;
0088       fEnergySpectrum_counts[i] = std::stoi(counts);
0089     }
0090   }
0091
0092   auto max_en_counter = fEnergySpectrum_counts[fEnergySpectrum_length - 1] - 1;
0093   auto rand_count = 1 + G4int(G4UniformRand() * max_en_counter);
0094
0095   // primitive search for lowest en_id that gives higher count value than
0096   // rand_count:
0097   G4int en_id = 0;
0098   for (; rand_count > fEnergySpectrum_counts[en_id]; en_id++)
0099     ;
0100
0101   G4int left = fEnergySpectrum_counts[en_id - 1];
0102   G4int right = fEnergySpectrum_counts[en_id];
0103
0104   G4double en_left = (en_id - 1) * fEnergySpectrum_gain + fEnergySpectrum_offset;
0105   G4double en_right = en_id * fEnergySpectrum_gain + fEnergySpectrum_offset;
0106
0107   // linear interpolation:
0108   G4double slope = (en_right - en_left) / (right - left);
0109   return (en_left + slope * (rand_count - left)) * MeV;
0110 }
0111
0112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0113
0114 G4ThreeVector PrimaryGeneratorAction::GenerateParticlePosition()
0115 {
0116   // the source is an infinitely thin disk of radius r
0117   G4double r = std::sqrt(G4UniformRand()) * fDetector->GetCollDiameter() / 2.;
0118   G4double phi = G4UniformRand() * twopi;
0119
0120   G4double x = fDetector->GetCollExitPosistion() - fDetector->GetCollLength();
0121   G4double y = r * std::cos(phi);
0122   G4double z = r * std::sin(phi);
0123
0124   return {x, y, z};
0125 }
0126
0127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0128
0129 G4ThreeVector PrimaryGeneratorAction::GenerateParticleDirection()
0130 {
0131   G4double phi, theta;
0132   G4double px, py, pz;
0133
0134   phi = G4UniformRand() * twopi;
0135
0136   // For convenience theta is measured along the beam axis, which is x-axis.
0137   // To reduce the number of events, when the projectile hits the collimator,
0138   // only forward angles in the range [0, max_theta] are considered.
0139   // If theta is larger, then the projectile will hit the collimator anyway.
0140   // Even with this restriction only 1 in 4 projectiles pass through the
0141   // collimator.
0142   G4double cos_max_theta =
0143     std::cos(std::atan(fDetector->GetCollDiameter() / fDetector->GetCollLength()));
0144   theta = std::acos(cos_max_theta + G4UniformRand() * (1 - cos_max_theta));
0145
0146   px = std::cos(theta);
0147   py = std::sin(theta) * std::sin(phi);
0148   pz = std::sin(theta) * std::cos(phi);
0149
0150   return {px, py, pz};
0151 }
0152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......