gammaray_telescope/src/GammaRayTelPrimaryGeneratorAction.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 // ------------------------------------------------------------
0028 //      GEANT 4 class implementation file
0029 //      CERN Geneva Switzerland
0030 //
0031 //
0032 //      ------------ GammaRayTelPrimaryGeneratorAction  ------
0033 //           by  G.Santin, F.Longo & R.Giannitrapani (13 nov 2000)
0034 //
0035 // ************************************************************
0036
0037 #include "G4RunManager.hh"
0038 #include "GammaRayTelPrimaryGeneratorAction.hh"
0039
0040 #include "GammaRayTelDetectorConstruction.hh"
0041 #include "GammaRayTelPrimaryGeneratorMessenger.hh"
0042
0043 #include "G4PhysicalConstants.hh"
0044 #include "G4SystemOfUnits.hh"
0045 #include "G4Event.hh"
0046 #include "G4ParticleGun.hh"
0047 #include "G4GeneralParticleSource.hh"
0048 #include "G4ParticleTable.hh"
0049 #include "G4ParticleDefinition.hh"
0050 #include "Randomize.hh"
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0053
0054 GammaRayTelPrimaryGeneratorAction::GammaRayTelPrimaryGeneratorAction() {
0055     detector = static_cast<const GammaRayTelDetectorConstruction*>(G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0056
0057     // create a messenger for this class
0058
0059     gunMessenger = new GammaRayTelPrimaryGeneratorMessenger(this);
0060
0061     constexpr auto NUMBER_OF_PARTICLES{1};
0062     particleGun = new G4ParticleGun(NUMBER_OF_PARTICLES);
0063
0064     // default particle kinematic
0065
0066     auto *particleTable = G4ParticleTable::GetParticleTable();
0067     auto *particle = particleTable->FindParticle("e-");
0068     particleGun->SetParticleDefinition(particle);
0069     particleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., -1.));
0070
0071     constexpr auto PARTICLE_ENERGY{30. * MeV};
0072     particleGun->SetParticleEnergy(PARTICLE_ENERGY);
0073
0074     auto position = 0.5 * (detector->GetWorldSizeZ());
0075     particleGun->SetParticlePosition(G4ThreeVector(0. * cm, 0. * cm, position));
0076     particleSource = new G4GeneralParticleSource();
0077 }
0078
0079 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0080
0081 GammaRayTelPrimaryGeneratorAction::~GammaRayTelPrimaryGeneratorAction() {
0082     delete particleGun;
0083     delete particleSource;
0084     delete gunMessenger;
0085 }
0086
0087 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0088
0089 void GammaRayTelPrimaryGeneratorAction::GeneratePrimaries(G4Event *event) {
0090     if (sourceGun) {
0091         G4cout << "Using G4ParticleGun... " << G4endl;
0092
0093         // this function is called at the begining of event
0094         //
0095         G4double x0 = 0. * cm;
0096         G4double y0 = 0. * cm;
0097         G4double z0 = 0.5 * (detector->GetWorldSizeZ());
0098
0099         G4ThreeVector pos0;
0100         auto vertex0 = G4ThreeVector(x0, y0, z0);
0101         auto momentumDirection0 = G4ThreeVector(0., 0., -1.);
0102
0103         G4double theta;
0104         G4double phi;
0105         G4double y = 0.;
0106         G4double f = 0.;
0107         G4double theta0 = 0.;
0108         G4double phi0 = 0.;
0109
0110         switch (sourceType) {
0111         case 0:
0112             particleGun->SetParticlePosition(vertex0);
0113             particleGun->SetParticleMomentumDirection(momentumDirection0);
0114             break;
0115         case 1:
0116             // GS: Generate random position on the 4PIsphere to create a uniform distribution
0117             // GS: on the sphere
0118             phi = G4UniformRand() * twopi;
0119             do {
0120                 y = G4UniformRand() * 1.0;
0121                 theta = G4UniformRand() * pi;
0122                 f = std::sin(theta);
0123             } while (y > f);
0124             vertex0 = G4ThreeVector(1., 0., 0.);
0125             vertex0.setMag(vertexRadius);
0126             vertex0.setTheta(theta);
0127             vertex0.setPhi(phi);
0128             particleGun->SetParticlePosition(vertex0);
0129
0130             momentumDirection0 = G4ThreeVector(1., 0., 0.);
0131
0132             do {
0133                 phi = G4UniformRand() * twopi;
0134                 do {
0135                     y = G4UniformRand() * 1.0;
0136                     theta = G4UniformRand() * pi;
0137                     f = std::sin(theta);
0138                 } while (y > f);
0139                 momentumDirection0.setPhi(phi);
0140                 momentumDirection0.setTheta(theta);
0141             } while (vertex0.dot(momentumDirection0) >= -0.7 * vertex0.mag());
0142
0143             particleGun->SetParticleMomentumDirection((G4ParticleMomentum) momentumDirection0);
0144
0145             break;
0146         case 2:
0147             // GS: Generate random position on the upper semi-sphere z > 0 to create a uniform distribution
0148             // GS: on a plane
0149             phi = G4UniformRand() * twopi;
0150
0151             do {
0152                 y = G4UniformRand() * 1.0;
0153                 theta = G4UniformRand() * halfpi;
0154                 f = std::sin(theta) * std::cos(theta);
0155             } while (y > f);
0156
0157             vertex0 = G4ThreeVector(1., 0., 0.);
0158
0159             auto xy = detector->GetWorldSizeXY();
0160             auto z = detector->GetWorldSizeZ();
0161
0162             if (vertexRadius > xy * 0.5) {
0163                 G4cout << "vertexRadius too big " << G4endl;
0164                 G4cout << "vertexRadius set to " << xy * 0.45 << G4endl;
0165                 vertexRadius = xy * 0.45;
0166             }
0167
0168             if (vertexRadius > z * 0.5) {
0169                 G4cout << "vertexRadius too high " << G4endl;
0170                 G4cout << "vertexRadius set to " << z * 0.45 << G4endl;
0171                 vertexRadius = z * 0.45;
0172             }
0173
0174             vertex0.setMag(vertexRadius);
0175             vertex0.setTheta(theta);
0176             vertex0.setPhi(phi);
0177
0178             // GS: Get the user defined direction for the primaries and
0179             // GS: Rotate the random position according to the user defined direction for the particle
0180
0181             momentumDirection0 = particleGun->GetParticleMomentumDirection();
0182             if (momentumDirection0.mag() > 0.001) {
0183                 theta0 = momentumDirection0.theta();
0184                 phi0 = momentumDirection0.phi();
0185             }
0186
0187             if (theta0 != 0.) {
0188                 G4ThreeVector rotationAxis(1., 0., 0.);
0189                 rotationAxis.setPhi(phi0 + halfpi);
0190                 vertex0.rotate(theta0 + pi, rotationAxis);
0191             }
0192             particleGun->SetParticlePosition(vertex0);
0193             break;
0194         }
0195
0196         constexpr auto INITIAL_PARTICLE_ENERGY{100 * MeV};
0197         G4double particleEnergy = INITIAL_PARTICLE_ENERGY;
0198
0199         switch (spectrumType) {
0200         case 0: // Uniform energy (1 GeV - 10 GeV)
0201             y = G4UniformRand();
0202             particleEnergy = y * 9.0 * GeV + 1.0 * GeV;
0203             G4cout << "Particle energy: " << particleEnergy / GeV << " LIN" << G4endl;
0204             break;
0205         case 1: // Logarithmic energy
0206             y = G4UniformRand();
0207             particleEnergy = std::pow(10, y) * GeV;
0208             G4cout << "Particle energy: " << particleEnergy / GeV << " LOG" << G4endl;
0209             break;
0210         case 2: // Power law (-4)
0211             do {
0212                 y = G4UniformRand() * 100000.0;
0213                 particleEnergy = G4UniformRand() * 10. * GeV;
0214                 f = std::pow(particleEnergy * (1 / GeV), -4.);
0215             } while (y > f);
0216             // particleGun->SetParticleEnergy(particleEnergy);
0217             break;
0218         case 3: // Monochromatic
0219             particleEnergy = particleGun->GetParticleEnergy();
0220             // 100 MeV;
0221             G4cout << "Particle energy: " << particleEnergy << " MONOCHROMATIC" << G4endl;
0222             break;
0223         }
0224         particleGun->SetParticleEnergy(particleEnergy);
0225         G4cout << "Particle: " << particleGun->GetParticleDefinition()->GetParticleName() << G4endl;
0226         particleGun->GeneratePrimaryVertex(event);
0227     } else {
0228         particleSource->GeneratePrimaryVertex(event);
0229     }
0230 }