GammaTherapy/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 medical/GammaTherapy/src/PrimaryGeneratorAction.cc
0028 /// \brief Implementation of the PrimaryGeneratorAction class
0029 //
0030
0031 //---------------------------------------------------------------------------
0032 //
0033 // ClassName:   PrimaryGeneratorAction
0034 //
0035 // Description: Generate primary beam
0036 //
0037 // Authors: V.Grichine, V.Ivanchenko
0038 //
0039 // Modified:
0040 //
0041 //----------------------------------------------------------------------------
0042 //
0043
0044 #include "PrimaryGeneratorAction.hh"
0045
0046 #include "DetectorConstruction.hh"
0047 #include "PrimaryGeneratorMessenger.hh"
0048
0049 #include "G4ParticleDefinition.hh"
0050 #include "G4ParticleGun.hh"
0051 #include "G4ParticleTable.hh"
0052 #include "G4PhysicalConstants.hh"
0053 #include "G4SystemOfUnits.hh"
0054 #include "Randomize.hh"
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0058
0059 PrimaryGeneratorAction::PrimaryGeneratorAction(DetectorConstruction* pDet) : fDetector(pDet)
0060 {
0061   InitializeMe();
0062 }
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0065
0066 void PrimaryGeneratorAction::InitializeMe()
0067 {
0068   fVerbose = fDetector->GetVerbose();
0069   fMessenger = new PrimaryGeneratorMessenger(this);
0070   fParticleGun = new G4ParticleGun();
0071   fCounter = 0;
0072   fX0 = 0.0;
0073   fY0 = 0.0;
0074   fZ0 = 0.0;
0075   fSigmaX = 1.5 * mm;
0076   fSigmaY = 1.5 * mm;
0077   fSigmaZ = 0.0;
0078   fSigmaE = 0.0;
0079   fRMax2 = 2.5 * 2.5 * mm * mm;
0080   fSigmaTheta = 0.0;
0081   //  fSigmaTheta = 0.17*degree;
0082   fMinCosTheta = 2.0;
0083   SetBeamEnergy(50.0 * MeV);
0084   fPosition = G4ThreeVector(fX0, fY0, fZ0);
0085   fDirection = G4ThreeVector(0.0, 0.0, 1.0);
0086   fGauss = true;
0087 }
0088
0089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0090
0091 PrimaryGeneratorAction::~PrimaryGeneratorAction()
0092 {
0093   delete fParticleGun;
0094   delete fMessenger;
0095 }
0096
0097 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0098
0099 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
0100 {
0101   fCounter++;
0102
0103   // Simulation of beam position
0104   G4double x = fX0;
0105   G4double y = fY0;
0106   G4double z = fDetector->GetGeneratorPosZ();
0107   do {
0108     if (0.0 < fSigmaX) {
0109       x = G4RandGauss::shoot(fX0, fSigmaX);
0110     }
0111     if (0.0 < fSigmaY) {
0112       y = G4RandGauss::shoot(fY0, fSigmaY);
0113     }
0114   } while (x * x + y * y > fRMax2);
0115
0116   fPosition = G4ThreeVector(x, y, z);
0117   fParticleGun->SetParticlePosition(fPosition);
0118
0119   // Simulation of beam direction
0120   G4double ux = fDirection.x();
0121   G4double uy = fDirection.y();
0122   G4double uz = fDirection.z();
0123
0124   // Beam particles are uniformly distributed over phi, cosTheta
0125   if (1.0 > fMinCosTheta) {
0126     uz = fMinCosTheta + (1.0 - fMinCosTheta) * G4UniformRand();
0127     ux = std::sqrt((1.0 - uz) * (1.0 + uz));
0128   }
0129   else if (fSigmaTheta > 0.0) {
0130     ux = G4RandGauss::shoot(0.0, fSigmaTheta);
0131     uz = std::sqrt((1.0 - ux) * (1.0 + ux));
0132   }
0133
0134   G4double phi = twopi * G4UniformRand();
0135   uy = ux;
0136   ux *= std::cos(phi);
0137   uy *= std::sin(phi);
0138   fDirection = G4ThreeVector(ux, uy, uz);
0139
0140   fParticleGun->SetParticleMomentumDirection(fDirection);
0141
0142   // Simulation of beam kinetic energy
0143   G4double kinEnergy = fEnergy;
0144
0145   if (fGauss == "flatE") {
0146     kinEnergy = fEnergy - fSigmaE + 2. * fSigmaE * G4UniformRand();
0147   }
0148   else if (0.0 < fSigmaE) {
0149     kinEnergy = fEnergy + G4RandGauss::shoot(0.0, fSigmaE);
0150   }
0151   fParticleGun->SetParticleEnergy(kinEnergy);
0152
0153   if (fVerbose > 0) {
0154     G4ParticleDefinition* particle = fParticleGun->GetParticleDefinition();
0155     G4String particleName = particle->GetParticleName();
0156     G4cout << "Event#  " << fCounter << "  Beam particle is generated by PrimaryGeneratorAction "
0157            << G4endl;
0158     G4cout << "ParticleName= " << particleName << "  PDGcode= " << particle->GetPDGEncoding()
0159            << std::setprecision(5) << "   KinEnergy(GeV)= " << fEnergy / GeV
0160            << "   x(mm)= " << x / mm << " y(mm)= " << y / mm << " z(mm)= " << z / mm
0161            << "   ux= " << ux << " uy= " << uy << " uz= " << uz << G4endl;
0162   }
0163
0164   fParticleGun->GeneratePrimaryVertex(anEvent);
0165 }
0166
0167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0168
0169 void PrimaryGeneratorAction::SetBeamEnergy(G4double val)
0170 {
0171   fEnergy = val;
0172   if (fEnergy < fDetector->GetMaxEnergy()) fDetector->SetMaxEnergy(fEnergy);
0173 }
0174
0175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....