dnadamage2/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 // This example is provided by the Geant4-DNA collaboration
0027 // Any report or published results obtained using the Geant4-DNA software
0028 // shall cite the following Geant4-DNA collaboration publication:
0029 // Med. Phys. 37 (2010) 4692-4708
0030 // J. Comput. Phys. 274 (2014) 841-882
0031 // The Geant4-DNA web site is available at http://geant4-dna.org
0032 //
0033 // Authors: J. Naoki D. Kondo (UCSF, US)
0034 //          J. Ramos-Mendez and B. Faddegon (UCSF, US)
0035 //
0036 /// \file PrimaryGeneratorAction.cc
0037 /// \brief Implementation of the PrimaryGeneratorAction class
0038
0039 #include "PrimaryGeneratorAction.hh"
0040
0041 #include "G4AffineTransform.hh"
0042 #include "G4LogicalVolumeStore.hh"
0043 #include "G4ParticleDefinition.hh"
0044 #include "G4ParticleGun.hh"
0045 #include "G4ParticleTable.hh"
0046 #include "G4RandomDirection.hh"
0047 #include "G4SystemOfUnits.hh"
0048 #include "G4VSolid.hh"
0049 #include "G4VoxelLimits.hh"
0050 #include "Randomize.hh"
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0053
0054 PrimaryGeneratorAction::PrimaryGeneratorAction() : G4VUserPrimaryGeneratorAction(), fParticleGun(0)
0055 {
0056   G4int n_particle = 1;
0057   fParticleGun = new G4ParticleGun(n_particle);
0058
0059   // default particle kinematic
0060   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
0061   G4ParticleDefinition* particle = particleTable->FindParticle("e-");
0062   fParticleGun->SetParticleDefinition(particle);
0063   fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 0.));
0064   fParticleGun->SetParticleEnergy(100 * keV);
0065   fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
0066
0067   fpSourceFileUI = new G4UIcmdWithAString("/fpGun/SourceFile", this);
0068   fpSourceFileUI->SetGuidance("Select the secondary e- source data file.");
0069
0070   fpVertexUI = new G4UIcmdWithAnInteger("/fpGun/PrimariesPerEvent", this);
0071   fpVertexUI->SetGuidance("Number of primary particles per event.");
0072 }
0073
0074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0075
0076 PrimaryGeneratorAction::~PrimaryGeneratorAction()
0077 {
0078   delete fpVertexUI;
0079   delete fParticleGun;
0080   delete fpSourceFileUI;
0081 }
0082
0083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0084
0085 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
0086 {
0087   for (G4int i = 0; i < fVertex; i++) {
0088     G4ThreeVector Position = SamplePosition();
0089     G4double Energy = SampleSpectrum();
0090     G4ThreeVector Direction = G4RandomDirection();
0091
0092     fParticleGun->SetParticleMomentumDirection(Direction);
0093     fParticleGun->SetParticlePosition(Position);
0094     fParticleGun->SetParticleEnergy(Energy);
0095     fParticleGun->GeneratePrimaryVertex(anEvent);
0096   }
0097 }
0098
0099 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0100
0101 void PrimaryGeneratorAction::SetNewValue(G4UIcommand* command, G4String newValue)
0102 {
0103   if (command == fpSourceFileUI) {
0104     fSourceFile = newValue;
0105     LoadSpectrum();
0106   }
0107
0108   if (command == fpVertexUI) {
0109     fVertex = fpVertexUI->GetNewIntValue(newValue);
0110   }
0111 }
0112
0113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0114
0115 void PrimaryGeneratorAction::LoadSpectrum()
0116 {
0117   fEnergyWeights.clear();
0118   fEnergies.clear();
0119
0120   std::ifstream SpectrumFile;
0121   SpectrumFile.open(fSourceFile);
0122   G4double x, y;
0123
0124   if (!SpectrumFile) {
0125     std::cout << "Source File: " + fSourceFile + " not found!!!" << std::endl;
0126     exit(1);
0127   }
0128
0129   else {
0130     while (SpectrumFile >> x >> y) {
0131       fEnergies.push_back(x);
0132       fEnergyWeights.push_back(y);
0133     }
0134   }
0135 }
0136
0137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0138
0139 G4double PrimaryGeneratorAction::SampleSpectrum()
0140 {
0141   G4double Energy = 1 * eV;
0142   G4int Bins = fEnergyWeights.size();
0143
0144   while (true) {
0145     G4double X = ((fEnergies[Bins - 1] - fEnergies[0]) * G4UniformRand()) + fEnergies[0];
0146     G4double Y = SampleProbability(X);
0147     if (G4UniformRand() < Y) {
0148       Energy = X * MeV;
0149       break;
0150     }
0151   }
0152
0153   if (Energy > 1 * MeV) Energy = 0.9999999 * MeV;
0154
0155   return Energy;
0156 }
0157
0158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0159
0160 G4double PrimaryGeneratorAction::SampleProbability(G4double X)
0161 {
0162   G4double Y = 0;
0163
0164   for (std::size_t i = 0; i < fEnergies.size() - 1; i++) {
0165     if (fEnergies[i] < X && X < fEnergies[i + 1]) {
0166       G4double M = (fEnergyWeights[i + 1] - fEnergyWeights[i]) / (fEnergies[i + 1] - fEnergies[i]);
0167       G4double B = fEnergyWeights[i] - (M * fEnergies[i]);
0168       Y = M * X + B;
0169     }
0170   }
0171
0172   return Y;
0173 }
0174
0175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0176
0177 G4ThreeVector PrimaryGeneratorAction::SamplePosition()
0178 {
0179   G4VSolid* solid = G4LogicalVolumeStore::GetInstance()->GetVolume("PlasmidVolume")->GetSolid();
0180   G4VoxelLimits voxelLimits;
0181   G4AffineTransform dontUse;
0182   G4double testX, testY, testZ;
0183   G4double xmin, xmax, ymin, ymax, zmin, zmax;
0184   solid->CalculateExtent(kXAxis, voxelLimits, dontUse, xmin, xmax);
0185   solid->CalculateExtent(kYAxis, voxelLimits, dontUse, ymin, ymax);
0186   solid->CalculateExtent(kZAxis, voxelLimits, dontUse, zmin, zmax);
0187
0188   while (true) {
0189     testX = G4RandFlat::shoot(xmin, xmax);
0190     testY = G4RandFlat::shoot(ymin, ymax);
0191     testZ = G4RandFlat::shoot(zmin, zmax);
0192     if (solid->Inside(G4ThreeVector(testX, testY, testZ)) == kInside) {
0193       break;
0194     }
0195   }
0196
0197   return G4ThreeVector(testX, testY, testZ);
0198 }
0199
0200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....