GammaRayTelPrimaryGeneratorAction.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// $Id: GammaRayTelPrimaryGeneratorAction.cc 100697 2016-10-31 11:32:35Z gcosmo $
// ------------------------------------------------------------
//      GEANT 4 class implementation file
//      CERN Geneva Switzerland
//
//
//      ------------ GammaRayTelPrimaryGeneratorAction  ------
//           by  G.Santin, F.Longo & R.Giannitrapani (13 nov 2000)
//
// ************************************************************

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

#include "G4RunManager.hh"
#include "GammaRayTelPrimaryGeneratorAction.hh"

#include "GammaRayTelDetectorConstruction.hh"
#include "GammaRayTelPrimaryGeneratorMessenger.hh"

#include "G4PhysicalConstants.hh"
#include "G4SystemOfUnits.hh"
#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4GeneralParticleSource.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "Randomize.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

GammaRayTelPrimaryGeneratorAction::GammaRayTelPrimaryGeneratorAction()
  :rndmFlag("off"),nSourceType(0),nSpectrumType(0),sourceGun(false)
{
  GammaRayTelDetector = static_cast<const GammaRayTelDetectorConstruction*>
    (G4RunManager::GetRunManager()->GetUserDetectorConstruction());

  //create a messenger for this class
  
  gunMessenger = new GammaRayTelPrimaryGeneratorMessenger(this);

  G4int n_particle = 1;

  particleGun  = new G4ParticleGun(n_particle);     
  // default particle kinematic
  
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  G4String particleName;
  G4ParticleDefinition* particle
    = particleTable->FindParticle(particleName="e-");
  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.));
  particleGun->SetParticleEnergy(30.*MeV);
  G4double position = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
  particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
  particleSource = new G4GeneralParticleSource();
   
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

GammaRayTelPrimaryGeneratorAction::~GammaRayTelPrimaryGeneratorAction()
{
 
  delete particleGun;
  delete particleSource;

  delete gunMessenger;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

void GammaRayTelPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
   if (sourceGun)
    {

      //this function is called at the begining of event
      // 
      G4double z0 = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
      G4double x0 = 0.*cm, y0 = 0.*cm;
  
      G4ThreeVector pos0;
      G4ThreeVector dir0;
      G4ThreeVector vertex0 = G4ThreeVector(x0,y0,z0);
      
      dir0 = G4ThreeVector(0.,0.,-1.);
      
      G4double theta, phi, y, f;
      G4double theta0=0.;
      G4double phi0=0.;
      
      switch(nSourceType) {
      case 0:
    particleGun->SetParticlePosition(vertex0);
    particleGun->SetParticleMomentumDirection(dir0);
    break;
      case 1:
    // GS: Generate random position on the 4PIsphere to create a unif. distrib.
    // GS: on the sphere
    phi = G4UniformRand() * twopi;
    do {
      y = G4UniformRand()*1.0;
      theta = G4UniformRand() * pi;
      f = std::sin(theta);
    } while (y > f);
    vertex0 = G4ThreeVector(1.,0.,0.);
    vertex0.setMag(dVertexRadius);
    vertex0.setTheta(theta);
    vertex0.setPhi(phi);
    particleGun->SetParticlePosition(vertex0);
    
    dir0 = G4ThreeVector(1.,0.,0.);
    do {
      phi = G4UniformRand() * twopi;
      do {
        y = G4UniformRand()*1.0;
        theta = G4UniformRand() * pi;
        f = std::sin(theta);
      } while (y > f);
      dir0.setPhi(phi);
      dir0.setTheta(theta);
    } while (vertex0.dot(dir0) >= -0.7 * vertex0.mag());
    particleGun->SetParticleMomentumDirection((G4ParticleMomentum)dir0);
    
    break;
      case 2:
    // GS: Generate random position on the upper semi-sphere z>0 to create a unif. distrib.
    // GS: on a plane
    phi = G4UniformRand() * twopi;
    do {
      y = G4UniformRand()*1.0;
      theta = G4UniformRand() * halfpi;
      f = std::sin(theta) * std::cos(theta);
    } while (y > f);
    vertex0 = G4ThreeVector(1.,0.,0.);
    
    G4double xy = GammaRayTelDetector->GetWorldSizeXY();
    G4double z = GammaRayTelDetector->GetWorldSizeZ();
    
    if (dVertexRadius > xy*0.5)
      { 
        G4cout << "vertexRadius too big " << G4endl;
        G4cout << "vertexRadius setted to " << xy*0.45 << G4endl;
        dVertexRadius = xy*0.45;
      }
    
    if (dVertexRadius > z*0.5)
      { 
        G4cout << "vertexRadius too high " << G4endl;
        G4cout << "vertexRadius setted to " << z*0.45 << G4endl;
        dVertexRadius = z*0.45;
      }
    
    
    vertex0.setMag(dVertexRadius);
    vertex0.setTheta(theta);
    vertex0.setPhi(phi);
    
    // GS: Get the user defined direction for the primaries and
    // GS: Rotate the random position according to the user defined direction for the particle
    
    dir0 = particleGun->GetParticleMomentumDirection();
    if (dir0.mag() > 0.001) 
      {
        theta0 = dir0.theta();
        phi0   = dir0.phi();   
      }
    
    if (theta0!=0.) 
      {
        G4ThreeVector rotationAxis(1.,0.,0.);
        rotationAxis.setPhi(phi0+halfpi);
        vertex0.rotate(theta0+pi,rotationAxis);
      }
    particleGun->SetParticlePosition(vertex0);
    break;
      }
      
      
      G4double pEnergy;
      
      switch(nSpectrumType) {
      case 0:
    break;
      case 1:
    break;
      case 2:
    do {
      y = G4UniformRand()*100000.0;
      pEnergy = G4UniformRand() * 10. * GeV;
      f = std::pow(pEnergy * (1/GeV), -4.);
    } while (y > f);
    
    particleGun->SetParticleEnergy(pEnergy);
    
    break;
      case 3:
    break;
      }
      G4cout << particleGun->GetParticleDefinition()->GetParticleName() << G4endl;
      particleGun->GeneratePrimaryVertex(anEvent);
    }
   else
     {
       particleSource->GeneratePrimaryVertex(anEvent);
     }
   
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....