musrPrimaryGeneratorAction.cc

Go to the documentation of this file.

#include "musrPrimaryGeneratorAction.hh"
#include "musrDetectorConstruction.hh"
#include "musrPrimaryGeneratorMessenger.hh"
#include "G4Event.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "Randomize.hh"
#include "G4ios.hh"
#include "musrParticleGun.hh"
#include "G4UnitsTable.hh"
#include "globals.hh"
#include "G4Gamma.hh"
#include "G4ThreeVector.hh"
#include "G4RunManager.hh"
#include "time.h"
#include <iomanip>

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
musrPrimaryGeneratorAction::musrPrimaryGeneratorAction(
                                            musrDetectorConstruction* musrDC)
  :musrDetector(musrDC),rndmFlag("off"),xvertex(0.),yvertex(0.),zvertex(0.), 
 xyvertexdefined(false), zvertexdefined(false),pxanglevertex(0.),pyanglevertex(0.),pvertexdefined(false),
   pMomentum(0.),pMomentumdefined(false),TheSetWidth(0.),Widthdefined(false)
{
  G4int n_particle = 1;
  particleGun  = new musrParticleGun(n_particle);
  
  //create a messenger for this class
  gunMessenger = new musrPrimaryGeneratorMessenger(this);

  // default particle kinematic

  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  G4ParticleDefinition* particle= particleTable->FindParticle("mu+");
  particleGun->SetParticleDefinition(particle);

} 

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

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

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

void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
 


 //this function is called at the begining of event
  // 


    G4double tau=0.000002197;
    G4double rnddecaytime = -tau*log(1-G4UniformRand());
    decaytime = rnddecaytime*s;
    particleGun->SetDecayTime(decaytime);


  //Set the initial particle position
     const double stdFWHMfactor=sqrt(2*log(2));
     G4double TheWidth;

     if(Widthdefined)
       {
    TheWidth=TheSetWidth;
       }
     else
       {
    TheWidth=40.*mm;
       }


    G4double xbeamwidthFWHM = TheWidth;
    G4double ybeamwidthFWHM = TheWidth;

    G4double x0,y0,z0 ;


    G4double xstdposition =  xbeamwidthFWHM/(2*stdFWHMfactor);
    G4double ystdposition =  ybeamwidthFWHM/(2*stdFWHMfactor);
    G4RandGauss* jRndGauss = new G4RandGauss(theEngine, 0, xstdposition);
    G4double xrndposition = jRndGauss->shoot(0, xstdposition);
    G4double yrndposition = jRndGauss->shoot(0, ystdposition);
    G4double xposition = xrndposition*mm;
    G4double yposition = yrndposition*mm;


    // temporary lines
    //          xposition = xbeamwidthFWHM/2.0;
    //  yposition = ybeamwidthFWHM/2.0;



        //      G4cout<<"Initial muon position = " << xposition <<", " 
        //    << yposition << "mm \n" <<G4endl;



  if(xyvertexdefined)
  {
    x0 = xvertex ;
    y0 = yvertex ;
  }
  else
  {
    x0 = xposition ;
    y0 = yposition ;
  }


  if(zvertexdefined)
  {
    z0 = zvertex ;
  }
  else
  {
    z0 = 0.*mm ;
  }


  //Set the initial particle momentum

  G4double pxangle0,pyangle0 ;

  if(pvertexdefined)
  {
    pxangle0 = pxanglevertex ;
    pyangle0 = pyanglevertex ;
  }
  else
  {
    pxangle0 = 0.*deg ;
    pyangle0 = 0.*deg ;
  }

  G4double myMomentum;

  if(pMomentumdefined)
    {
    myMomentum=pMomentum;      
    }
  else
    {
      myMomentum=0.0;
    }


  particleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));



    G4double xpitch = sin(pxangle0);
    G4double ypitch = sin(pyangle0);
    G4double zpitch = sqrt(1-(sin(pxangle0)*sin(pxangle0)
                           +sin(pyangle0)*sin(pyangle0)));

    // Then orient the spin relative to the beam pitch


  //longitudinal
    particleGun->SetParticlePolarization(G4ThreeVector(xpitch,ypitch,-zpitch));

  //transverse
  //     particleGun->SetParticlePolarization(G4ThreeVector(zpitch,0,-xpitch));


  G4double energy = sqrt(1e6*myMomentum*myMomentum+11163.82428)-105.659; 
  particleGun->SetParticleEnergy(energy*MeV);
  particleGun->SetParticleMomentumDirection(G4ThreeVector(xpitch,ypitch,zpitch));  particleGun->GeneratePrimaryVertex(anEvent);



}

//
//


void musrPrimaryGeneratorAction::Setxvertex(G4double x)
{
  xyvertexdefined = true ;
  xvertex = x ;

}

void musrPrimaryGeneratorAction::Setyvertex(G4double y)
{
  xyvertexdefined = true ;
  yvertex = y ;

}
void musrPrimaryGeneratorAction::Setzvertex(G4double z)
{
  zvertexdefined = true ;
  zvertex = z ;

}

void musrPrimaryGeneratorAction::Setpxanglevertex(G4double px)
{
  pvertexdefined = true ;
  pxanglevertex = px ;

}

void musrPrimaryGeneratorAction::Setpyanglevertex(G4double py)
{
  pvertexdefined = true ;
  pyanglevertex = py ;

}
void musrPrimaryGeneratorAction::SetMomentum(G4double p)
{
  pMomentumdefined = true ;
  pMomentum = p ;

}


void musrPrimaryGeneratorAction::SetBeamWidth(G4double w)
{
  Widthdefined = true;
  TheSetWidth=w;
G4cout<< "FWHM beam width set to "<< TheSetWidth <<"mm" <<G4endl;

}

---