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