#include <GammaRayTelPrimaryGeneratorAction.hh>

Inheritance diagram for GammaRayTelPrimaryGeneratorAction:

Public Member Functions
	GammaRayTelPrimaryGeneratorAction ()

	~GammaRayTelPrimaryGeneratorAction ()

void	GeneratePrimaries (G4Event *)

void	SetRndmFlag (G4String val)

void	SetSourceType (G4int val)

void	SetSpectrumType (G4int val)

void	SetVertexRadius (G4double val)

void	SetSourceGen (G4bool val)

Public Member Functions inherited from G4VUserPrimaryGeneratorAction
	G4VUserPrimaryGeneratorAction ()

virtual	~G4VUserPrimaryGeneratorAction ()

Detailed Description

Definition at line 58 of file GammaRayTelPrimaryGeneratorAction.hh.

Constructor & Destructor Documentation

GammaRayTelPrimaryGeneratorAction::GammaRayTelPrimaryGeneratorAction ( )

Definition at line 58 of file GammaRayTelPrimaryGeneratorAction.cc.

References python.hepunit::cm, G4ParticleTable::FindParticle(), G4ParticleTable::GetParticleTable(), G4RunManager::GetRunManager(), G4RunManager::GetUserDetectorConstruction(), GammaRayTelDetectorConstruction::GetWorldSizeZ(), python.hepunit::MeV, G4ParticleGun::SetParticleDefinition(), G4ParticleGun::SetParticleEnergy(), G4ParticleGun::SetParticleMomentumDirection(), and G4VPrimaryGenerator::SetParticlePosition().

   :rndmFlag("off"),nSourceType(0),nSpectrumType(0)
 {
   G4RunManager* runManager = G4RunManager::GetRunManager();
   GammaRayTelDetector =
     (GammaRayTelDetectorConstruction*)(runManager->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();
    
 }

GammaRayTelPrimaryGeneratorAction::~GammaRayTelPrimaryGeneratorAction ( )

Definition at line 90 of file GammaRayTelPrimaryGeneratorAction.cc.

 {
   if (sourceGun)
     delete particleGun;
   else
     delete particleSource;
 
   delete gunMessenger;
 }

Member Function Documentation

void GammaRayTelPrimaryGeneratorAction::GeneratePrimaries ( G4Event * anEvent )

virtual

Implements G4VUserPrimaryGeneratorAction.

Definition at line 102 of file GammaRayTelPrimaryGeneratorAction.cc.

 {
    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;
       }
       
       particleGun->GeneratePrimaryVertex(anEvent);
     }
    else
      {
        particleSource->GeneratePrimaryVertex(anEvent);
      }
    
 }