Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/dna/dnadamage2/src/PrimaryGeneratorAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // This example is provided by the Geant4-DNA collaboration
 27 // Any report or published results obtained using the Geant4-DNA software
 28 // shall cite the following Geant4-DNA collaboration publication:
 29 // Med. Phys. 37 (2010) 4692-4708
 30 // J. Comput. Phys. 274 (2014) 841-882
 31 // The Geant4-DNA web site is available at http://geant4-dna.org
 32 //
 33 // Authors: J. Naoki D. Kondo (UCSF, US)
 34 //          J. Ramos-Mendez and B. Faddegon (UCSF, US)
 35 //
 36 /// \file PrimaryGeneratorAction.cc
 37 /// \brief Implementation of the PrimaryGeneratorAction class
 38 
 39 #include "PrimaryGeneratorAction.hh"
 40 
 41 #include "G4AffineTransform.hh"
 42 #include "G4LogicalVolumeStore.hh"
 43 #include "G4ParticleDefinition.hh"
 44 #include "G4ParticleGun.hh"
 45 #include "G4ParticleTable.hh"
 46 #include "G4RandomDirection.hh"
 47 #include "G4SystemOfUnits.hh"
 48 #include "G4VSolid.hh"
 49 #include "G4VoxelLimits.hh"
 50 #include "Randomize.hh"
 51 
 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 53 
 54 PrimaryGeneratorAction::PrimaryGeneratorAction() : G4VUserPrimaryGeneratorAction(), fParticleGun(0)
 55 {
 56   G4int n_particle = 1;
 57   fParticleGun = new G4ParticleGun(n_particle);
 58 
 59   // default particle kinematic
 60   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
 61   G4ParticleDefinition* particle = particleTable->FindParticle("e-");
 62   fParticleGun->SetParticleDefinition(particle);
 63   fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 0.));
 64   fParticleGun->SetParticleEnergy(100 * keV);
 65   fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));
 66 
 67   fpSourceFileUI = new G4UIcmdWithAString("/fpGun/SourceFile", this);
 68   fpSourceFileUI->SetGuidance("Select the secondary e- source data file.");
 69 
 70   fpVertexUI = new G4UIcmdWithAnInteger("/fpGun/PrimariesPerEvent", this);
 71   fpVertexUI->SetGuidance("Number of primary particles per event.");
 72 }
 73 
 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 75 
 76 PrimaryGeneratorAction::~PrimaryGeneratorAction()
 77 {
 78   delete fpVertexUI;
 79   delete fParticleGun;
 80   delete fpSourceFileUI;
 81 }
 82 
 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
 84 
 85 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
 86 {
 87   for (G4int i = 0; i < fVertex; i++) {
 88     G4ThreeVector Position = SamplePosition();
 89     G4double Energy = SampleSpectrum();
 90     G4ThreeVector Direction = G4RandomDirection();
 91 
 92     fParticleGun->SetParticleMomentumDirection(Direction);
 93     fParticleGun->SetParticlePosition(Position);
 94     fParticleGun->SetParticleEnergy(Energy);
 95     fParticleGun->GeneratePrimaryVertex(anEvent);
 96   }
 97 }
 98 
 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
100 
101 void PrimaryGeneratorAction::SetNewValue(G4UIcommand* command, G4String newValue)
102 {
103   if (command == fpSourceFileUI) {
104     fSourceFile = newValue;
105     LoadSpectrum();
106   }
107 
108   if (command == fpVertexUI) {
109     fVertex = fpVertexUI->GetNewIntValue(newValue);
110   }
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
114 
115 void PrimaryGeneratorAction::LoadSpectrum()
116 {
117   fEnergyWeights.clear();
118   fEnergies.clear();
119 
120   std::ifstream SpectrumFile;
121   SpectrumFile.open(fSourceFile);
122   G4double x, y;
123 
124   if (!SpectrumFile) {
125     std::cout << "Source File: " + fSourceFile + " not found!!!" << std::endl;
126     exit(1);
127   }
128 
129   else {
130     while (SpectrumFile >> x >> y) {
131       fEnergies.push_back(x);
132       fEnergyWeights.push_back(y);
133     }
134   }
135 }
136 
137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
138 
139 G4double PrimaryGeneratorAction::SampleSpectrum()
140 {
141   G4double Energy = 1 * eV;
142   G4int Bins = fEnergyWeights.size();
143 
144   while (true) {
145     G4double X = ((fEnergies[Bins - 1] - fEnergies[0]) * G4UniformRand()) + fEnergies[0];
146     G4double Y = SampleProbability(X);
147     if (G4UniformRand() < Y) {
148       Energy = X * MeV;
149       break;
150     }
151   }
152 
153   if (Energy > 1 * MeV) Energy = 0.9999999 * MeV;
154 
155   return Energy;
156 }
157 
158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
159 
160 G4double PrimaryGeneratorAction::SampleProbability(G4double X)
161 {
162   G4double Y = 0;
163 
164   for (std::size_t i = 0; i < fEnergies.size() - 1; i++) {
165     if (fEnergies[i] < X && X < fEnergies[i + 1]) {
166       G4double M = (fEnergyWeights[i + 1] - fEnergyWeights[i]) / (fEnergies[i + 1] - fEnergies[i]);
167       G4double B = fEnergyWeights[i] - (M * fEnergies[i]);
168       Y = M * X + B;
169     }
170   }
171 
172   return Y;
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
176 
177 G4ThreeVector PrimaryGeneratorAction::SamplePosition()
178 {
179   G4VSolid* solid = G4LogicalVolumeStore::GetInstance()->GetVolume("PlasmidVolume")->GetSolid();
180   G4VoxelLimits voxelLimits;
181   G4AffineTransform dontUse;
182   G4double testX, testY, testZ;
183   G4double xmin, xmax, ymin, ymax, zmin, zmax;
184   solid->CalculateExtent(kXAxis, voxelLimits, dontUse, xmin, xmax);
185   solid->CalculateExtent(kYAxis, voxelLimits, dontUse, ymin, ymax);
186   solid->CalculateExtent(kZAxis, voxelLimits, dontUse, zmin, zmax);
187 
188   while (true) {
189     testX = G4RandFlat::shoot(xmin, xmax);
190     testY = G4RandFlat::shoot(ymin, ymax);
191     testZ = G4RandFlat::shoot(zmin, zmax);
192     if (solid->Inside(G4ThreeVector(testX, testY, testZ)) == kInside) {
193       break;
194     }
195   }
196 
197   return G4ThreeVector(testX, testY, testZ);
198 }
199 
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....