MollerPrimaryGenAction.cc

Go to the documentation of this file.

#include "MollerPrimaryGenAction.hh"
#include "MollerPrimaryGeneratorMessenger.hh"

#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4Electron.hh"
#include "globals.hh"
#include "TRandom2.h"


MollerPrimaryGenAction::MollerPrimaryGenAction() : ebeam(11.000), mp(0.9382796), bremcut(1.0e-6), me(0.511E-3)
{

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

  eventnumber=0; // used only in optics generator?

  G4int n_particle = 1;
  particleGun = new G4ParticleGun(n_particle);

  //set particle
  //G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  G4ParticleDefinition* particle = G4Electron::ElectronDefinition();
  
  G4cout << "Electron gun created!!\n" << G4endl;

  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
  particleGun->SetParticleEnergy(11.0*GeV);
  particleGun->SetParticlePosition(G4ThreeVector(0.*cm, 0.*cm, -125.*cm));
}

MollerPrimaryGenAction::~MollerPrimaryGenAction()
{
  delete particleGun;
  delete gunMessenger;
}

void MollerPrimaryGenAction::GeneratePrimaries(G4Event* anEvent)
{
  if(anEvent->GetEventID()==0){
  TRandom2 *random_num = new TRandom2(0);
  long particle_seed = (long) random_num->GetSeed();
  particle_seed += time(0);

  if (!fSetSeedConst) {
    G4cout<<"Setting the seed to a random value..."<<anEvent->GetEventID()<<G4endl;
    CLHEP::HepRandom::setTheSeed(particle_seed);
  }
  if (fSetSeedConst) {
    G4cout<<"Setting the seed to a constant value..."<<anEvent->GetEventID()<<G4endl;
    //CLHEP::HepRandom::setTheSeed(194738798457);  //First event too low angle
    //CLHEP::HepRandom::setTheSeed(734687698678); //First event too high angle
    CLHEP::HepRandom::setTheSeed(94528347598798); 
  }

  }
  // To add an additional generator:
  // 1. add an extra case statement here
  // 2. create another generator method in a separate file
  // 3. edit the GenCmd by altering the SetGuidance and SetRange method calls
  switch(GeneratorNumber) {
  case 0:
    GeneratePrimaries_moller_dustin(anEvent);
    break;
  case 1:
    GeneratePrimaries_elasticep_dustin(anEvent);
    break;
  case 2:
    GeneratePrimaries_inelasticep_dustin(anEvent);
    break;
  case 3:
    GeneratePrimaries_phasespace_mollers(anEvent);
    break;
  case 4:
    GeneratePrimaries_phasespace_motts(anEvent);
    break;
  default:
    G4cerr << "Generator number " << GeneratorNumber << " does not exist!\n";
    G4cerr << "Choose which generator to use.\n";
    G4cerr << "  Choice : 0. Moller electrons (Dustin) [default]\n";
    G4cerr << "           1. Elastic ep electrons (Dustin)\n";
    G4cerr << "           2. Inelastic ep electrons (Dustin)\n";
    G4cerr << "           3. Optics rays electrons (Moller scattering)\n";
    G4cerr << "           4. Optics rays electrons (Mott scattering)\n";
    G4cerr << "Running with default.\n";
    GeneratePrimaries_moller_dustin(anEvent);
    break;
  }
}

double MollerPrimaryGenAction::RadProfile(double eloss, double btt){
    double Ekin = ebeam - me;
    return 1./eloss*(1.-eloss/Ekin+0.75*pow(eloss/Ekin,2))*pow(eloss/Ekin,btt);
}

//add numerical integration of energy loss
double MollerPrimaryGenAction::EnergNumInt(double btt, double a0, double b0)
{
    //double de = (ebeam-bremcut)/((double) nbin);
    double sum = 0.0;

    int j;
    double boolc[5] = {7.0, 32.0, 12.0, 32.0, 7.0};

    double a, b, thissum;

    for(int i =0;i<nbin;i++) {
  // Integrate over sample spacings that are logarithmic
  a = bremcut*pow(b0/a0, (((double) i)/nbin));
  b = bremcut*pow(b0/a0, (((double) i+1.0)/nbin));

  // Boole's rule
  thissum = 0.0;
  for( j = 0; j < 5; j++ ){
      thissum +=  boolc[j]*RadProfile( (b-a)*j*0.25 + a, btt);
  }
  sum += thissum*(b-a)/90.0;
    }

    return sum;
}