MollerPrimaryGenAction Class Reference

#include <MollerPrimaryGenAction.hh>

Collaboration diagram for MollerPrimaryGenAction:

[legend]

Public Member Functions
	MollerPrimaryGenAction ()
	~MollerPrimaryGenAction ()
void	GeneratePrimaries (G4Event *anEvent)
void	GeneratePrimaries_moller_dustin (G4Event *anEvent)
void	GeneratePrimaries_elasticep_dustin (G4Event *anevent)
void	GeneratePrimaries_inelasticep_dustin (G4Event *anevent)
void	GeneratePrimaries_phasespace_mollers (G4Event *anevent)
void	GeneratePrimaries_phasespace_motts (G4Event *anevent)
void	GeneratePrimaries_pion (G4Event *anevent)
void	SetGenerator (G4int val)
G4int	GetGenerator ()
double	EnergNumInt (double btt, double a0, double b0)
double	RadProfile (double eloss, double btt)
void	setSeedValue (G4bool value)
Double_t	wiser_sigma (Double_t ebeam, Double_t pf, Double_t thf, Double_t rl, Int_t type)
Static Public Member Functions
static Double_t	wiser_all_fit (Double_t *x, Double_t *par)
Private Attributes
G4ParticleGun *	particleGun
MollerPrimaryGeneratorMessenger *	gunMessenger
G4int	GeneratorNumber
G4int	eventnumber
cms_cross_section_t	XStable [nbin]
double_differential_t	DDtable [nbinIneEP][nbinIneEP]
const G4double	ebeam
const G4double	mp
const G4double	bremcut
const G4double	me
G4bool	fSetSeedConst
Static Private Attributes
static const G4int	nbin = 1000
static const G4int	nbinIneEP = 100

---

Detailed Description

Definition at line 28 of file MollerPrimaryGenAction.hh.

---

Constructor & Destructor Documentation

MollerPrimaryGenAction::MollerPrimaryGenAction

(

)

Definition at line 13 of file MollerPrimaryGenAction.cc.

References eventnumber, gunMessenger, and particleGun.

                                               : 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

(

)

Definition at line 36 of file MollerPrimaryGenAction.cc.

References gunMessenger, and particleGun.

{
  delete particleGun;
  delete gunMessenger;
}

---

Member Function Documentation

double MollerPrimaryGenAction::EnergNumInt	(	double	btt,
		double	a0,
		double	b0
	)

Definition at line 101 of file MollerPrimaryGenAction.cc.

References bremcut, nbin, and RadProfile().

Referenced by GeneratePrimaries_elasticep_dustin().

{
    //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;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void MollerPrimaryGenAction::GeneratePrimaries

(

G4Event *

anEvent

)

Definition at line 42 of file MollerPrimaryGenAction.cc.

References fSetSeedConst, GeneratePrimaries_elasticep_dustin(), GeneratePrimaries_inelasticep_dustin(), GeneratePrimaries_moller_dustin(), GeneratePrimaries_phasespace_mollers(), GeneratePrimaries_phasespace_motts(), and GeneratorNumber.

{
  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;
  }
}

Here is the call graph for this function:

void MollerPrimaryGenAction::GeneratePrimaries_elasticep_dustin

(

G4Event *

anevent

)

Definition at line 41 of file MollerPrimaryGenAction_elasticep.cc.

References alpha, bremcut, d2r, ebeam, EnergNumInt(), Euler, gevfm, gRootAnalysis, me, mp, nbin, NINTERVAL, particleGun, pi, pmu, RadProfile(), RootAnalysis::SetMomentum0(), RootAnalysis::SetMomentum1(), RootAnalysis::SetRate(), RootAnalysis::SetTotXS(), cms_cross_section_t::theta, cms_cross_section_t::xs, XStable, and cms_cross_section_t::y.

Referenced by GeneratePrimaries().

{
  //----------------------
  // GenerateEvent
  //----------------------
  
  double initial_x = 0.0;//*cm;
  double initial_y = 0.0;//*cm;
  double initial_z = 0.0;//*cm;

  double initial_px = 0.0;//*GeV;
  double initial_py = 0.0;//*GeV;
  double initial_pz = 0.0;//*GeV;

  //Find the scattering point in the target
  double vertex; //scattering point in target (cm) measured from the upstream end of target
  double tt = 150.0; // t = thickness of target in cm
  double rho = 0.0715; // rho = density of target in g/cm**3 (from KK)
  double radiationLength = 63.04; // g/cm**2 (from PDG)
  double t_rad = radiationLength/rho; // t_rad = radiation length in cm
  double z = G4UniformRand();
  vertex = z*tt;//just do this for now
  
  double Evertex = ebeam;
  //Now we know how mush material the beam travels through before moller interaction
  //Use this to calculate energy loss
  double  Ekin = ebeam - me;
  double fracrl = vertex/t_rad; //total radiation lenght fraction travelled through
  double bt = fracrl*4./3.;

  double prob, prob_sample, sample, eloss, value;
  //Calculation of probability to have bremsstrahlung effect above 1 keV
  prob = 1.- pow(bremcut/Ekin,bt) - bt/(bt+1.)*(1.- pow(bremcut/Ekin,bt+1.))
    + 0.75*bt/(2.+bt)*(1.- pow(bremcut/Ekin,bt+2.));
  prob = prob/(1.- bt*Euler + bt*bt/2.*(Euler*Euler+pi*pi/6.)); /* Gamma function */
  prob_sample = G4UniformRand();  /* Random sampling */

  value = 1.0;

  double colEcut = 0.080;  // Expected cutoff for collimator

  double Evlo[NINTERVAL] = {
      bremcut,
      (ebeam-bremcut-colEcut)*1.0/(11.0-colEcut-bremcut),
      (ebeam-bremcut-colEcut)*8.0/(11.0-colEcut-bremcut),
      (ebeam-bremcut-colEcut)*(ebeam-colEcut)/(11.0-colEcut-bremcut),
  };

  double Evhi[NINTERVAL] = {
      (ebeam-bremcut-colEcut)*1.0/(11.0-colEcut-bremcut),
      (ebeam-bremcut-colEcut)*8.0/(11.0-colEcut-bremcut),
      (ebeam-bremcut-colEcut)*(ebeam-colEcut)/(11.0-colEcut-bremcut),
      (ebeam-bremcut-colEcut)*ebeam/(11.0-colEcut-bremcut)
  };

  double Eprob[NINTERVAL]  = { 0.07, 0.05, 0.875, 0.005};

  double Enorm[NINTERVAL];
  // Interval normalization
  for( int idx = 0; idx < NINTERVAL; idx++ ){
      Enorm[idx]  = ((Evhi[idx]-Evlo[idx])/(Evhi[NINTERVAL-1]-Evlo[0]))
    /Eprob[idx];
  }

  int Evidx;
  double evsum = 0.0;

  // Averages over the intervals
  double vweight = 1.0;
  double vavg[4] = {
      log(Evhi[0]/Evlo[0])/(Evhi[0]-Evlo[0]),
      1.0,
      log((ebeam-Evlo[2])/(ebeam-Evhi[2]))/(Evhi[2]-Evlo[2]),
      1.0};

  if (prob_sample <= prob) {//Bremsstrahlung has taken place!
      //  We break this into 4 seperate energy loss intervals
      //  with total integrals roughly the size of
      //  what the ep product looks like with 11 GeV beam
      //   cut  -  1000 MeV, 1/x, 7%
      //   1000 -  8000 MeV, flat 5%
      //   8000 - 10920 MeV, 1/(E-x), 87.5%
      //  10920 - 11000 MeV, flat, 0.5%

      sample = G4UniformRand();

      // Identify our region
      //
      // based on the probability distribution
      Evidx = 0;
      evsum  = Eprob[Evidx];
      while( evsum < sample ){
    Evidx++;
    evsum += Eprob[Evidx];
      }

      sample = G4UniformRand();

      if( Evidx == 1 || Evidx == 3 ){
    // sample energy flat
    eloss = (Evhi[Evidx]-Evlo[Evidx])*sample + Evlo[Evidx];
    vweight = 1.0;
      }

      if( Evidx == 0 ){
    eloss = Evlo[Evidx]*pow(Evhi[Evidx]/Evlo[Evidx],sample);
    vweight = eloss;
      }

      if( Evidx == 2 ){
    eloss = ebeam - (ebeam-Evhi[Evidx])*
        pow((ebeam-Evlo[Evidx])/(ebeam-Evhi[Evidx]),sample);
    vweight = (ebeam-eloss);
      }

      vweight *= vavg[Evidx];
      //  mult by ebeam-bremcut for proper normalization
      value = RadProfile(eloss, bt)*
    ((Evhi[NINTERVAL-1]-Evlo[0])/EnergNumInt(bt, Evlo[0], Evhi[NINTERVAL-1])) // average of RadProfile
    *vweight // sampling weighting (flat or ) / average value for normalization
    *Enorm[Evidx]; //  Weight given the region

      Evertex = ebeam - eloss;
  }

  if( Evertex < me ) return;

  double vertex_x = (-0.25 + (0.5*(G4UniformRand())));
  double vertex_y = (-0.25 + (0.5*(G4UniformRand())));
  initial_x = vertex_x;
  initial_y = vertex_y;
  initial_z = vertex - tt/2.; //cm  
  
  //Make realistic beam angle (Not yet)
  //double theta_div = 0.000125*((double)random()/(double)RAND_MAX);
  //double phi_div = 2.0*M_PI*((double)random()/(double)RAND_MAX);
  ////double theta_div = 0.0;
  ////double phi_div = 0.0;
  ////cout << "theta = "<<theta_div*180/M_PI<<", phi = "<<phi_div*180/M_PI<<endl;
    
  initial_px = 0.0;
  initial_py = 0.0;
  initial_pz = Evertex; // GeV
  //cout <<"-------------elastic ep Event-------------------"<<endl;
  //cout <<" vertex = "<<vertex<<" , Evertex = "<<Evertex<<endl;

  //Lab System
  double thetaMott[2];
  thetaMott[0] = 0.1;//degree in lab
  thetaMott[1] = 2.0;//degree in lab
    
  double p0 = sqrt(Evertex*Evertex - me*me);

  for(int i =0;i<nbin;i++) {
    double theta = thetaMott[0] + i*(thetaMott[1]-thetaMott[0])/(nbin-1);
    double costh = cos(theta*d2r);
    double sinth = sin(theta*d2r);
    double sin2 = (1. - costh)/2.;
    double cos2 = 1. - sin2;
    double p = mp*Evertex/(Evertex*(1. - costh) + mp);
    double e = sqrt(p*p+me*me);
    double q2 = -2.*(me*me - Evertex*e + p0*p*costh);
    double t = q2/4./(mp*mp);
    double ge = 1./pow((1. + q2/0.71),2);
    double gm = ge*pmu;
    double cang = (ge*ge + t*gm*gm)/(1. + t) + t*t*gm*gm*sin2/cos2;
    double crmott = alpha*alpha*cos2/4./(p0*p0)/(sin2*sin2)/(1. + 2.*p0/mp*sin2)*gevfm*gevfm/100.;
    double cross = crmott*cang;
    double csmott = cross*sinth;
    XStable[i].y = theta;
    XStable[i].theta = theta;
    XStable[i].xs = csmott;
  }

  cms_cross_section_t *mott = XStable;
    
  //---- Convert cross section table to integral fractions ----
  // First, numerically integrate
  double Im=0.0;
  double save = mott[0].xs;
  mott[0].xs = 0.0;

  for(int j = 1;j<nbin-1;j++) {
    Im += (save+mott[j].xs)/2.0;
    save = mott[j].xs;
    mott[j].xs = Im;
  }
  // Next, normalize tables
  mott = XStable;

  assert( Im > 0.0 );

  for(int k=0;k<nbin;k++){
    mott[k].xs /= Im;
  }
  double totXS = Im*(thetaMott[1]-thetaMott[0])*d2r/(nbin-1.0)*2*pi;
  double factor = 6.0221415*0.1*1.0E-6/1.602176487E-19;
  double rate = totXS*tt*rho/1.008*1.0*factor; //Event rate/sec/ uA
  rate *= value; //mult by value to renorm rate to uneven energy distribution
  gRootAnalysis->SetTotXS(totXS);
  gRootAnalysis->SetRate(rate);
    
  // Determine mott cms scattering angle using ss
  double ss = (G4UniformRand());
    
  //cms_cross_section_t *c=XStable;
  int index;
  double slope, b;
  for(index=0;index<nbin-1;index++) if(ss <= mott[index+1].xs) break;
  assert( index < nbin-1 );
  slope = (mott[index+1].xs - mott[index].xs)/(mott[index+1].theta - mott[index].theta);
  b = mott[index].xs;
    
  double lab_theta = mott[index].theta + (ss - b)/slope;
  double lab_phi = (2.0*pi*(G4UniformRand())) - pi;
    
  //Now add unknown beam divergence
  //pi0_theta += theta_div;
  //pi0_phi += phi_div;
        
  // Define kinematic conditions in lab
  double costh = cos(lab_theta*d2r);
  double sinth = sin(lab_theta*d2r);
  double p = mp*Evertex/(Evertex*(1. - costh) + mp);    
  double e = sqrt(p*p + me*me);
  double ep = Evertex + mp - e;
  //double ep_alt = sqrt(mp*mp+p0*p0+p*p-2*p0*p*costh);
  double pp = sqrt(ep*ep - mp*mp);
  double sinp = p*sinth/pp;
  //double thetp = asin(sinp);
  double cosp = sqrt(1.- sinp*sinp);
  double cosphie = cos( lab_phi );
  double sinphie = sin( lab_phi );
  double cosphip = cos( lab_phi - pi);
  double sinphip = sin( lab_phi - pi);

  //G4cout <<"Evertex + mp = "<<(Evertex+mp)<<", theta = "<<lab_theta<<", p = "<<p<<", e = "<<e<<", pp = "<<pp<<", ep = "<<ep<<", e+ep = "<<(e+ep)<<"\n";

  double e1, p1[3], e2, p2[3];
  e1 = e;
  e2 = ep;
  for(int i=0;i<3;i++) {
    p1[i] = 0.0;
    p2[i] = 0.0;
  }
    
  p1[2] = p*costh;
  p2[2] = pp*cosp;
  p1[0] = p*sinth*cosphie;
  p1[1] = p*sinth*sinphie;
  p2[0] = pp*sinp*cosphip;
  p2[1] = pp*sinp*sinphip;
    
  double pmag1 = p;//sqrt(p1[0]*p1[0]+p1[1]*p1[1]+p1[2]*p1[2]);
  //double pmag2 = pp;//sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]);
  double pmag0 = sqrt(initial_px*initial_px+initial_py*initial_py+initial_pz*initial_pz);
    
  // ------------ Final state particles ---------------
    
  //G4int n_particle = 1;
  //particleGun = new G4ParticleGun(n_particle);
    
  G4ParticleDefinition* particle = G4Electron::ElectronDefinition();
  //We will only throw the electron for now

  assert( pmag1 > 0.0 );

  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleMomentumDirection(G4ThreeVector(p1[0]/pmag1,p1[1]/pmag1,p1[2]/pmag1));
  particleGun->SetParticleEnergy(e1*1000);
  particleGun->SetParticlePosition(G4ThreeVector(initial_x*10, initial_y*10, initial_z*10));
  particleGun->GeneratePrimaryVertex(anEvent);

  gRootAnalysis->SetMomentum0(Evertex*1000,initial_px/pmag0,initial_py/pmag0,initial_pz/pmag0);
  gRootAnalysis->SetMomentum1(e1*1000,p1[0]/pmag1,p1[1]/pmag1,p1[2]/pmag1);

  //anEvent->SetEventID(0);
  anEvent->GetPrimaryVertex(0)->GetPrimary(0);//->SetTrackID(99999);

}

Here is the call graph for this function:

Here is the caller graph for this function:

void MollerPrimaryGenAction::GeneratePrimaries_inelasticep_dustin

(

G4Event *

anevent

)

Definition at line 41 of file MollerPrimaryGenAction_inelasticep.cc.

References alpha, d2r, DDtable, ebeam, double_differential_t::Eprime, Euler, gRootAnalysis, me, mp, nbinIneEP, particleGun, pi, RootAnalysis::SetMomentum0(), RootAnalysis::SetMomentum1(), RootAnalysis::SetRate(), RootAnalysis::SetTotXS(), SigT(), double_differential_t::theta, and double_differential_t::xs.

Referenced by GeneratePrimaries().

{
  //----------------------
  // GenerateEvent
  //----------------------
  
  double initial_x = 0.0;//*cm;
  double initial_y = 0.0;//*cm;
  double initial_z = 0.0;//*cm;

  double initial_px = 0.0;//*GeV;
  double initial_py = 0.0;//*GeV;
  double initial_pz = 0.0;//*GeV;

  //Find the scattering point in the target
  double vertex; //scattering point in target (cm) measured from the upstream end of target
  double t = 150.0; // t = thickness of target in cm
  double rho = 0.0715; // rho = density of target in g/cm**3 (from KK)
  double radiationLength = 63.04; // g/cm**2 (from PDG)
  double t_rad = radiationLength/rho; // t_rad = radiation length in cm
  double z = G4UniformRand();
  vertex = z*t;//just do this for now
  
  double Evertex = ebeam;
  //Now we know how mush material the beam travels through before moller interaction
  //Use this to calculate energy loss
  double cut=0.000001; //lower limit of bremsstrahlung (1 keV)
  double  Ekin = ebeam - me;
  double fracrl = vertex/t_rad; //total radiation lenght fraction travelled through
  double bt = fracrl*4./3.;

  double prob, prob_sample, sample, eloss, env, value, ref;
  //Calculation of probability to have bremsstrahlung effect above 1 keV
  prob = 1.- pow(cut/Ekin,bt) - bt/(bt+1.)*(1.- pow(cut/Ekin,bt+1.))
    + 0.75*bt/(2.+bt)*(1.- pow(cut/Ekin,bt+2.));
  prob = prob/(1.- bt*Euler + bt*bt/2.*(Euler*Euler+pi*pi/6.)); /* Gamma function */

  prob_sample = G4UniformRand();  /* Random sampling */

  if (prob_sample <= prob) {//Bremsstrahlung has taken place!
    //Generate photon energy with sample and reject using 1/x as envelope 
    do {
      sample = G4UniformRand();
      eloss = cut*pow(Ekin/cut,sample);
      
      env = 1./eloss;
      value = 1./eloss*(1.-eloss/Ekin+0.75*pow(eloss/Ekin,2))*pow(eloss/Ekin,bt);
      
      sample = G4UniformRand();
      ref = value/env;
    } while (sample > ref);
    
    Evertex = ebeam - eloss;
  }

  if( Evertex < me ) return;

  double vertex_x = (-0.25 + (0.5*(G4UniformRand())));
  double vertex_y = (-0.25 + (0.5*(G4UniformRand())));
  initial_x = vertex_x;
  initial_y = vertex_y;
  initial_z = vertex - t/2.; //cm  
  
  //Make realistic beam angle (Not yet)
  //double theta_div = 0.000125*((double)random()/(double)RAND_MAX);
  //double phi_div = 2.0*M_PI*((double)random()/(double)RAND_MAX);
  ////double theta_div = 0.0;
  ////double phi_div = 0.0;
  ////cout << "theta = "<<theta_div*180/M_PI<<", phi = "<<phi_div*180/M_PI<<endl;
    
  initial_px = 0.0;
  initial_py = 0.0;
  initial_pz = Evertex; // GeV
  //cout <<"------------inelastic ep Event-------------------"<<endl;
  //cout <<" vertex = "<<vertex<<" , Evertex = "<<Evertex<<endl;
    
  /*based in fit to sigma_T and assuming R=sig_L/sig_T=q2 (good to q2=0.2 or so) */

  //Lab System
  double thetaIneEP[2];
  thetaIneEP[0] = 0.1;//0.2578;//degree in lab
  thetaIneEP[1] = 2.0;//0.5157;//degree in lab

  double e = Evertex;
  double eprime_max = e - 0.150;
  double eprime_min = eprime_max/nbinIneEP;
  if(Evertex<eprime_min) G4cout <<"Evertex is less than eprime_min..."<<"\n";
  else if(Evertex<eprime_max) eprime_max = Evertex;

  if( eprime_min < me ){
      eprime_min = me;
  }

  for(int i=0;i<nbinIneEP;i++) {//loop over theta
    for (int j=0;j<nbinIneEP;j++) {//loop over eprime
      double theta = (thetaIneEP[1] - thetaIneEP[0])/(nbinIneEP-1)*i + thetaIneEP[0];
      double costh = cos(theta*d2r/2.);//Note theta/2 here
      double sinth = sin(theta*d2r/2.);//and here
      double sin2 = sinth*sinth;//(1. - costh)/2.;
      double cos2 = costh*costh;//1. - sin2;
      double eprime = (eprime_max - eprime_min)/(nbinIneEP-1)*j + eprime_min;
      double q2 = 4.*e*eprime*sin2;
      double x = q2/2./mp/(e-eprime);
      double nu = q2/2./mp/x;
      //double wsq = mp*mp + 2.*mp*(e-eprime) - q2;
      double k = nu - q2/2./mp;
      double sigt=SigT(k);/* get sigma_T */
      /* Get F_1 using virtual photon flux */
      double f1 = (nu - q2 / 2. / mp) * mp / 4. / (pi*pi) / alpha * sigt /0.389 ;
      //Note that 0.389 = (hbar*c)^2 in GeV^2 mbarn, a conversion constant
      /* calculate f_2 using R=q2 */
      double r = q2;
      double f2 = f1 / ( (1. + (2. * mp * x) * (2. * mp * x) /q2) / 2. / x / (1. + r )) ; 
      /* sig is in barns/GeV/sr */
      double sig=5.2E-9/e/e/sin2/sin2* (cos2 * f2 / nu + sin2 * 2. * f1 / mp); 
      sig *= sin(theta*d2r);//with phase space factor 
      DDtable[i][j].Eprime = eprime;
      DDtable[i][j].theta = theta;
      DDtable[i][j].xs = sig;
    }
  }

  double_differential_t XStable_theta[nbinIneEP];
  double_differential_t XStable_Eprime[nbinIneEP];
  //---- Convert cross section table to integral fractions ----
  // First, numerically integrate
  double Im=0.0;
  double save;

  for(int i = 0;i<nbinIneEP; i++) {
    XStable_theta[i].xs   = 0.0;
    XStable_theta[i].theta = (thetaIneEP[1] - thetaIneEP[0])/(nbinIneEP-1)*i + thetaIneEP[0];
    for(int j = 0;j<nbinIneEP-1;j++){
  XStable_theta[i].xs += (DDtable[i][j].xs + DDtable[i][j+1].xs)/2.0;
    }//Integrate over Eprime
    if( i == 0 ){ 
  save = XStable_theta[i].xs;
  XStable_theta[i].xs = 0.0;
    } else {
  Im += (save+XStable_theta[i].xs)/2.0;
  save = XStable_theta[i].xs;
  XStable_theta[i].xs = Im;
    }
  }
  double totXS = Im*(thetaIneEP[1]-thetaIneEP[0])*d2r*(eprime_max - eprime_min)/(nbinIneEP-1.0)/(nbinIneEP-1.0)*2*pi;
  double factor = 6.0221415*0.1*1.0E-6/1.602176487E-19;
  double rate = totXS*t*rho/1.008*1.0*factor; //Event rate/sec/ uA
  gRootAnalysis->SetTotXS(totXS);
  gRootAnalysis->SetRate(rate);

  if( !(Im > 0.0) ){
      // Cross section is zero
      // just abort

      return;
  }

  assert( Im > 0.0 );

  // Next, normalize tables
  for(int k=0;k<nbinIneEP;k++){
    XStable_theta[k].xs /= Im;
  }

  // Determine scattering angle using ss
  double ss = (G4UniformRand());
  int index;
  double slope, b;
  for(index=0;index<nbinIneEP-1;index++) if(ss <= XStable_theta[index+1].xs) break;
  assert( index < nbinIneEP-1 );

  slope = (XStable_theta[index+1].xs - XStable_theta[index].xs)/(XStable_theta[index+1].theta - XStable_theta[index].theta);
  b = XStable_theta[index].xs;
  double lab_theta = XStable_theta[index].theta + (ss - b)/slope;
  int theIndex  = index;
  // Dimensionless quantity on how far we are to the next bin
  // This is for our bilinear interpolation
  double intpar = (ss - b)/(XStable_theta[index+1].xs - XStable_theta[index].xs);

  assert( 0.0 <= intpar && intpar < 1.0 );

  //Now do the Eprime 
  //For this, simply use the array of the selected theta
  double Im2=0.0;
  // Binlinear interpolation
  save =  DDtable[theIndex][0].xs*(1.0-intpar) + DDtable[theIndex+1][0].xs*intpar;
  XStable_Eprime[0].xs = 0.0;
  XStable_Eprime[0].Eprime = eprime_min;

  for(int j = 1;j<nbinIneEP;j++) {
      // Binlinear interpolation
    Im2 += (save + DDtable[theIndex][j].xs*(1.0-intpar) + DDtable[theIndex+1][j].xs*intpar)/2.0;
    save = DDtable[theIndex][j].xs*(1.0-intpar) + DDtable[theIndex+1][j].xs*intpar;
    XStable_Eprime[j].xs = Im2;
    XStable_Eprime[j].Eprime = (eprime_max - eprime_min)/(nbinIneEP-1)*j + eprime_min;
  }

  assert( Im2 > 0.0 );

  //double_differential_t *ineEP_Eprime = XStable_Eprime;
  // Next, normalize tables
  for(int k=0;k<nbinIneEP;k++){
    XStable_Eprime[k].xs /= Im2;
  }
  // Determine Eprime using ss2
  double ss2 = (G4UniformRand());
  int index2;
  double slope2, b2;
  for(index2=0;index2<nbinIneEP-1;index2++) if(ss2 <= XStable_Eprime[index2+1].xs) break;
  assert( index2 < nbinIneEP-1 );

  slope2 = (XStable_Eprime[index2+1].xs - XStable_Eprime[index2].xs)/(XStable_Eprime[index2+1].Eprime - XStable_Eprime[index2].Eprime);
  b2 = XStable_Eprime[index2].xs;
  double lab_Eprime = XStable_Eprime[index2].Eprime + (ss2 - b2)/slope2;

  //Now do phi
  double lab_phi = (2.0*pi*(G4UniformRand())) - pi;
    
  //Now add unknown beam divergence
  //pi0_theta += theta_div;
  //pi0_phi += phi_div;

  // Define kinematic conditions in lab
  double costh = cos(lab_theta*d2r);
  double sinth = sin(lab_theta*d2r);
  e = lab_Eprime;
  assert( e > 0.0 );
  double p = sqrt(e*e - me*me);
  assert( p > 0.0 );
  //double ep = Evertex + mp - e;
  double cosphie = cos( lab_phi );
  double sinphie = sin( lab_phi );
  //double cosphip = cos( lab_phi - pi);
  //double sinphip = sin( lab_phi - pi);

  //G4cout <<"Evertex + mp = "<<(Evertex+mp)<<", theta = "<<lab_theta<<", p = "<<p<<", e = "<<e<<", pp = "<<pp<<", ep = "<<ep<<", e+ep = "<<(e+ep)<<"\n";
  //G4cout <<"Evertex = "<<(Evertex)<<", theta = "<<lab_theta<<", p = "<<p<<", e = "<<e<<"\n";

  double e1, p1[3];
  e1 = e;
  for(int i=0;i<3;i++) {
    p1[i] = 0.0;
  }
    
  p1[2] = p*costh;
  p1[0] = p*sinth*cosphie;
  p1[1] = p*sinth*sinphie;
    
  double pmag1 = p;//sqrt(p1[0]*p1[0]+p1[1]*p1[1]+p1[2]*p1[2]);
  double pmag0 = sqrt(initial_px*initial_px+initial_py*initial_py+initial_pz*initial_pz);
    
  // ------------ Final state particles ---------------
    
  //G4int n_particle = 1;
  //particleGun = new G4ParticleGun(n_particle);
    
  G4ParticleDefinition* particle_e = G4Electron::ElectronDefinition();
  //We will only throw the electron for now    
  particleGun->SetParticleDefinition(particle_e);
  particleGun->SetParticleMomentumDirection(G4ThreeVector(p1[0]/pmag1,p1[1]/pmag1,p1[2]/pmag1));
  particleGun->SetParticleEnergy(e1*1000);
  particleGun->SetParticlePosition(G4ThreeVector(initial_x*10, initial_y*10, initial_z*10));
  particleGun->GeneratePrimaryVertex(anEvent);
    
  gRootAnalysis->SetMomentum0(Evertex*1000,initial_px/pmag0,initial_py/pmag0,initial_pz/pmag0);
  gRootAnalysis->SetMomentum1(e1*1000,p1[0]/pmag1,p1[1]/pmag1,p1[2]/pmag1);

  //anEvent->SetEventID(0);
  anEvent->GetPrimaryVertex(0)->GetPrimary(0);//->SetTrackID(99999);
    
}

Here is the call graph for this function:

Here is the caller graph for this function:

void MollerPrimaryGenAction::GeneratePrimaries_moller_dustin

(

G4Event *

anEvent

)

Definition at line 38 of file MollerPrimaryGenAction_moller.cc.

References d2r, ebeam, Euler, gRootAnalysis, me, nbin, particleGun, pi, r0, RootAnalysis::SetMomentum0(), RootAnalysis::SetMomentum1(), RootAnalysis::SetMomentum2(), RootAnalysis::SetRate(), RootAnalysis::SetTotXS(), cms_cross_section_t::theta, cms_cross_section_t::xs, XStable, and cms_cross_section_t::y.

Referenced by GeneratePrimaries().

{
  //----------------------
  // GenerateEvent
  //----------------------
  
  double initial_x = 0.0;//*cm;
  double initial_y = 0.0;//*cm;
  double initial_z = 0.0;//*cm;

  double initial_px = 0.0;//*GeV;
  double initial_py = 0.0;//*GeV;
  double initial_pz = 0.0;//*GeV;

  //Find the scattering point in the target
  double vertex; //scattering point in target (cm) measured from the upstream end of target
  double t = 150.0; // t = thickness of target in cm
  double rho = 0.0715; // rho = density of target in g/cm**3 (from KK)
  double radiationLength = 63.04; // g/cm**2 (from PDG)
  double t_rad = radiationLength/rho; // t_rad = radiation length in cm
  double z = G4UniformRand();
  vertex = z*t;//just do this for now
  
  double Evertex = ebeam;
  //Now we know how mush material the beam travels through before moller interaction
  //Use this to calculate energy loss
  double cut=0.000001; //lower limit of bremsstrahlung (1 keV)
  double  Ekin = ebeam - me;
  double fracrl = vertex/t_rad; //total radiation lenght fraction travelled through
  double bt = fracrl*4./3.;

  double prob, prob_sample, sample, eloss, env, value, ref;
  //Calculation of probability to have bremsstrahlung effect above 1 keV
  prob = 1.- pow(cut/Ekin,bt) - bt/(bt+1.)*(1.- pow(cut/Ekin,bt+1.))
    + 0.75*bt/(2.+bt)*(1.- pow(cut/Ekin,bt+2.));
  prob = prob/(1.- bt*Euler + bt*bt/2.*(Euler*Euler+pi*pi/6.)); /* Gamma function */

  prob_sample = G4UniformRand();  /* Random sampling */

  if (prob_sample <= prob) {//Bremsstrahlung has taken place!
    //Generate photon energy with sample and reject using 1/x as envelope 
    do {
      sample = G4UniformRand();
      eloss = cut*pow(Ekin/cut,sample);
      
      env = 1./eloss;
      value = 1./eloss*(1.-eloss/Ekin+0.75*pow(eloss/Ekin,2))*pow(eloss/Ekin,bt);
      
      sample = G4UniformRand();
      ref = value/env;
    } while (sample > ref);
    
    Evertex = ebeam - eloss;
  }

  if( Evertex < me ) return;

  double vertex_x = (-0.25 + (0.5*(G4UniformRand())));
  double vertex_y = (-0.25 + (0.5*(G4UniformRand())));
  initial_x = vertex_x;
  initial_y = vertex_y;
  initial_z = vertex - t/2.; //cm  
  
  //Make realistic beam angle (Not yet)
  //double theta_div = 0.000125*((double)random()/(double)RAND_MAX);
  //double phi_div = 2.0*M_PI*((double)random()/(double)RAND_MAX);
  ////double theta_div = 0.0;
  ////double phi_div = 0.0;
  ////cout << "theta = "<<theta_div*180/M_PI<<", phi = "<<phi_div*180/M_PI<<endl;
    
  initial_px = 0.0;
  initial_py = 0.0;
  initial_pz = Evertex; // GeV
  
  //cout <<"-------------Moller Event-------------------"<<endl;
  //cout <<" vertex = "<<vertex<<" , Evertex = "<<Evertex<<endl;
        
  //Center of Mass System
  double thetaMoller[2];
  thetaMoller[0] = 53.13010;//costheta_cm = +0.6
  thetaMoller[1] = 126.8699;//costheta_cm = -0.6
  double x1 = cos(thetaMoller[1]*d2r);
  double x2 = cos(thetaMoller[0]*d2r);
  double y2 = (1 - x1)/2.;
  double y1 = (1 - x2)/2.;
  double ds = me*me+me*me+2.*me*Evertex;
  double ecm = sqrt(ds);
  double decm[2];
  decm[0] = ds/2./ecm;
  double dps = pow(decm[0],2)-me*me;
  double dp = sqrt(dps);
  decm[1] = sqrt(dps+me*me);

  for(int i =0;i<nbin;i++) {
    double y = y1 + ((y2 - y1)/nbin*i)+(y2 - y1)/nbin/2.;
    double dt = -1.*y*ds;
    double ct = (dt-me*me-me*me+2.*pow(decm[0],2))/(2.*pow(dp,2));
    if( !( -1.0 < ct && ct < 1.0 )){ 
  // SPR 9/20/12
  // Something's wrong here for very low beam energies
  // Not sure how to handle it, but at least this point
  // is kinematically forbidden
  return;
    }
    double ct2 = pow(ct,2);
    double st2 = 1. - ct2;
    double st = sqrt(st2);
    double th = acos(ct);
      
    double ecmp = ecm/2.;
    double vcmp = sqrt(pow(ecmp,2)-me*me)/ecmp;
    double dsig = 0.;
      
    double eom = ecmp/me;
    double eom2 = pow(eom,2);
    double d1 = pow(r0,2)*pow((2.*eom2-1.),2)/4./pow(vcmp,4)/pow(eom2,3);
    double d2 = 4./pow(st2,2)-3./st2+pow((eom2-1.),2)/pow((2.*eom2-1),2)*(1+4./pow(st,2));
      
    dsig = d1*d2;

    assert( dsig > 0.0 );

    dsig = dsig*2.; //Convert from dsigma/dOmega to dsigma/dy/dphi
    dsig = dsig/100; //Convert from fm**2 to Barns
    XStable[i].y = y;
    XStable[i].theta = th;
    XStable[i].xs = dsig;
  }

  cms_cross_section_t *moller = XStable;
    
  //---- Convert cross section table to integral fractions ----
  // First, numerically integrate
  double Im=0.0;
  double save = moller[0].xs;
  moller[0].xs = Im;

  for(int j = 1;j<nbin;j++) {
      Im += (save + moller[j].xs)/2.0;
      save = moller[j].xs;
      moller[j].xs = Im;
  }

  double totXS = Im*(y2-y1)/(nbin-1.0)*2*pi;
  double factor = 6.0221415*0.1*1.0E-6/1.602176487E-19;
  double rate = totXS*t*rho/1.008*1.0*factor; //Event rate/sec/ uA
  gRootAnalysis->SetTotXS(totXS);
  gRootAnalysis->SetRate(rate);

  if( !(Im > 0.0) ){ 
      printf("Im = %f, Evertex = %f\n", Im, Evertex);
  }
  assert( Im > 0.0 ); // idiot check

  // Next, normalize tables
  moller = XStable;
  for(int k=0;k<nbin;k++){
    moller[k].xs /= Im;
  }
    
  // Determine moller cms scattering angle using ss
  double ss = (G4UniformRand());
    
  //cms_cross_section_t *c=XStable;
  int index;
  double slope, b;
  for(index=0;index<nbin-1;index++) if(ss <= moller[index+1].xs) break;
  assert( index < nbin-1 ); // idiot check
  slope = (moller[index+1].xs - moller[index].xs)/(moller[index+1].theta - moller[index].theta);
  b = moller[index].xs;
    
  double cms_theta = moller[index].theta + (ss - b)/slope;

  double cms_phi = (2.0*pi*(G4UniformRand())) - pi;
    
  //Now add unknown beam divergence
  //pi0_theta += theta_div;
  //pi0_phi += phi_div;
        
  // Define kinematic conditions in cms
  double gammap = Evertex/ecm;
  double betap = 1/gammap*sqrt(pow(gammap,2)-1);
    
  double e1cm = gammap*me;
  double e2cm = e1cm;
    
  //Primary (beam) electron
  double costheta1 = cos( cms_theta );
  double sintheta1 = sin( cms_theta );
    
  double cosphi1 = cos( cms_phi );
  double sinphi1 = sin( cms_phi );
    
  double p1cm[3];


  p1cm[0] = e1cm * sintheta1 * cosphi1;
  p1cm[1] = e1cm * sintheta1 * sinphi1;
  p1cm[2] = e1cm * costheta1;
    
  //Secondary (target) electron
  double costheta2 = cos( pi - cms_theta );
  double sintheta2 = sin( pi - cms_theta );
    
  double cosphi2 = cos( cms_phi - pi);
  double sinphi2 = sin( cms_phi - pi);
  if (cms_phi<0.0){
    cosphi2 = cos( cms_phi + pi);
    sinphi2 = sin( cms_phi + pi);
  } else {
      
  }
    
  double p2cm[3];  
  p2cm[0] = e2cm * sintheta2 * cosphi2;
  p2cm[1] = e2cm * sintheta2 * sinphi2;
  p2cm[2] = e2cm * costheta2;
    
    
  // Transform to the lab frame.
  double e1, p1[3], e2, p2[3];
  e1 = 0.0;
  e2 = 0.0;
  for(int i=0;i<3;i++) {
    p1[i] = 0.0;
    p2[i] = 0.0;
  }
    
  p1[2] = gammap * (p1cm[2] + betap*e1cm);
  p2[2] = gammap * (p2cm[2] + betap*e2cm);
  p1[0] = p1cm[0];
  p1[1] = p1cm[1];
  p2[0] = p2cm[0];
  p2[1] = p2cm[1];
    
  double pmag1 = sqrt(p1[0]*p1[0]+p1[1]*p1[1]+p1[2]*p1[2]);
  double pmag2 = sqrt(p2[0]*p2[0]+p2[1]*p2[1]+p2[2]*p2[2]);
  double pmag0 = sqrt(initial_px*initial_px+initial_py*initial_py+initial_pz*initial_pz);
    
  e1 = gammap * ( e1cm + betap*p1cm[2] );
  e2 = gammap * ( e2cm + betap*p2cm[2] );
    
  // Calculate the laboratory polar angles for moller electrons.  
  //double th1 = acos( p1[2] / e1 );
  //double ph1 = atan2( p1[1], p1[0] );
    
  //double th2 = acos( p2[2] / e2 );
  //double ph2 = atan2( p2[1], p2[0] );
    
    
  // ------------ Final state electrons ---------------
    
  //G4int n_particle = 1;
  //particleGun = new G4ParticleGun(n_particle);
    
  G4ParticleDefinition* particle = G4Electron::ElectronDefinition();

  assert( pmag1 > 0.0 && pmag2 > 0.0 && pmag0 > 0.0 );

  assert( e1 > 0.0 && e2 > 0.0 );
    
  if(e1>=e2) {
    particleGun->SetParticleDefinition(particle);
    particleGun->SetParticleMomentumDirection(G4ThreeVector(p1[0]/pmag1,p1[1]/pmag1,p1[2]/pmag1));
    particleGun->SetParticleEnergy(e1*1000);
    particleGun->SetParticlePosition(G4ThreeVector(initial_x*10, initial_y*10, initial_z*10));
    particleGun->GeneratePrimaryVertex(anEvent);
      
    particleGun->SetParticleDefinition(particle);
    particleGun->SetParticleMomentumDirection(G4ThreeVector(p2[0]/pmag2,p2[1]/pmag2,p2[2]/pmag2));
    particleGun->SetParticleEnergy(e2*1000);
    particleGun->SetParticlePosition(G4ThreeVector(initial_x*10, initial_y*10, initial_z*10));
    particleGun->GeneratePrimaryVertex(anEvent);
      
    gRootAnalysis->SetMomentum0(Evertex*1000,initial_px/pmag0,initial_py/pmag0,initial_pz/pmag0);
    gRootAnalysis->SetMomentum1(e1*1000,p1[0]/pmag1,p1[1]/pmag1,p1[2]/pmag1);
    gRootAnalysis->SetMomentum2(e2*1000,p2[0]/pmag2,p2[1]/pmag2,p2[2]/pmag2);
  } else {
    particleGun->SetParticleDefinition(particle);
    particleGun->SetParticleMomentumDirection(G4ThreeVector(p2[0]/pmag2,p2[1]/pmag2,p2[2]/pmag2));
    particleGun->SetParticleEnergy(e2*1000);
    particleGun->SetParticlePosition(G4ThreeVector(initial_x*10, initial_y*10, initial_z*10));
    particleGun->GeneratePrimaryVertex(anEvent);
      
    particleGun->SetParticleDefinition(particle);
    particleGun->SetParticleMomentumDirection(G4ThreeVector(p1[0]/pmag1,p1[1]/pmag1,p1[2]/pmag1));
    particleGun->SetParticleEnergy(e1*1000);
    particleGun->SetParticlePosition(G4ThreeVector(initial_x*10, initial_y*10, initial_z*10));
    particleGun->GeneratePrimaryVertex(anEvent);    
      
    gRootAnalysis->SetMomentum0(Evertex*1000,initial_px/pmag0,initial_py/pmag0,initial_pz/pmag0);
    gRootAnalysis->SetMomentum2(e1*1000,p1[0]/pmag1,p1[1]/pmag1,p1[2]/pmag1);
    gRootAnalysis->SetMomentum1(e2*1000,p2[0]/pmag2,p2[1]/pmag2,p2[2]/pmag2);
  }
  //anEvent->SetEventID(0);
  anEvent->GetPrimaryVertex(0)->GetPrimary(0);//->SetTrackID(99999);

}

Here is the call graph for this function:

Here is the caller graph for this function:

void MollerPrimaryGenAction::GeneratePrimaries_phasespace_mollers

(

G4Event *

anevent

)

Definition at line 13 of file MollerPrimaryGenAction_phasespace.cc.

References ebeam, eventnumber, gRootAnalysis, particleGun, pi, RootAnalysis::SetMomentum0(), RootAnalysis::SetMomentum1(), RootAnalysis::SetMomentum2(), and RootAnalysis::SetProcess().

Referenced by GeneratePrimaries().

{
  //This generator throws events in phase space

  G4int numsectors = 7; 
  G4double sectorhalfopeningangle = pi/numsectors/2.;
  G4double etheta=0., ephi=0., zpos=0.;

  G4int numthetas=17, numphis=13, numzpos=2, numphicentre=7;
  G4double thetamin=0.004, thetamax=0.020;
  G4double zposmin=-0.75, zposmax=0.75;
  G4double phicentre=0;
  G4double phimin=0, phimax=2*pi;
  G4double phicentremin=0, phicentremax=12./7*pi;

  int thetanum, ephinum, zposnum, phicentrenum;
  double thetafrac, ephifrac, zposfrac, phicentrefrac;
  // step through theta then phi then target position z then sector   
  thetanum = eventnumber%(numthetas);  // % gives remainder after division
  ephinum = (eventnumber/(numthetas))%(numphis);
  zposnum = (eventnumber/(numthetas*numphis))%(numzpos);
  phicentrenum = (eventnumber/(numthetas*numphis*numzpos))%(numphicentre);

  thetafrac = ((double)thetanum)/(numthetas-1);
  ephifrac = ((double)ephinum)/(numphis-1);
  zposfrac = ((double)zposnum)/(numzpos-1);
  phicentrefrac = ((double)phicentrenum)/(numphicentre-1);

  etheta = thetamin + thetafrac*(thetamax-thetamin);
  phicentre = phicentremin + phicentrefrac*(phicentremax-phicentremin);
  phimin = phicentre - sectorhalfopeningangle*0.999;
  phimax = phicentre + sectorhalfopeningangle*0.999;
  ephi = phimin + ephifrac*(phimax-phimin);
  zpos = zposmin + zposfrac*(zposmax-zposmin);

  G4double eenergy;
  eenergy = (2*ebeam*0.000511)/(ebeam*(etheta*etheta) + 2*0.000511); // Moller electrons

  G4double xdir = sin(etheta)*cos(ephi);
  G4double ydir = sin(etheta)*sin(ephi);
  G4double zdir = cos(etheta);
  printf("Primaries:%5i phi:%2i %.2f zpos:%2i %.2f; phicentre:%2i %.2f; theta:%2i %.2f =%.3f deg =%4.1f mrad  eenergy=%.3f GeV  phi=%6.1f deg, xdir=%.3f, ydir=%.3f, zdir=%.3f\n", 
       eventnumber, ephinum, ephifrac, zposnum, zposfrac, phicentrenum, phicentrefrac, 
       thetanum, thetafrac, etheta/degree, etheta/rad*1000, eenergy, ephi/degree, xdir, ydir, zdir);
  
  particleGun->SetParticleEnergy(eenergy*GeV);
  particleGun->SetParticleMomentumDirection(G4ThreeVector(xdir,ydir,zdir));
  particleGun->SetParticlePosition(G4ThreeVector(0.*cm, 0.*cm, zpos*m));
  particleGun->GeneratePrimaryVertex(anEvent);
  eventnumber++;

  gRootAnalysis->SetMomentum0(11000,0,0,11000);
  gRootAnalysis->SetMomentum2(eenergy*MeV,eenergy*MeV*xdir,eenergy*MeV*ydir,eenergy*MeV*zdir);
  gRootAnalysis->SetMomentum1(eenergy*MeV,eenergy*MeV*xdir,eenergy*MeV*ydir,eenergy*MeV*zdir);
  gRootAnalysis->SetProcess(1);

  if (eventnumber==(numphis*numthetas*numzpos*7)) {
    printf("\nFully completed one set of optics rays.\n\n");
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void MollerPrimaryGenAction::GeneratePrimaries_phasespace_motts

(

G4Event *

anevent

)

Definition at line 13 of file MollerPrimaryGenAction_phasespace_motts.cc.

References ebeam, eventnumber, gRootAnalysis, mp, particleGun, pi, RootAnalysis::SetMomentum0(), RootAnalysis::SetMomentum1(), RootAnalysis::SetMomentum2(), and RootAnalysis::SetProcess().

Referenced by GeneratePrimaries().

{
  //This generator throws events in phase space

  G4int numsectors = 7; 
  G4double sectorhalfopeningangle = pi/numsectors/2.;
  G4double etheta=0., ephi=0., zpos=0.;

  G4int numthetas=17, numphis=13, numzpos=2, numphicentre=7;
  G4double thetamin=0.004, thetamax=0.020;
  G4double zposmin=-0.75, zposmax=0.75;
  G4double phicentre=0;
  G4double phimin=0, phimax=2*pi;
  G4double phicentremin=0, phicentremax=12./7*pi;

  int thetanum, ephinum, zposnum, phicentrenum;
  double thetafrac, ephifrac, zposfrac, phicentrefrac;
  // step through theta then phi then target position z then sector   
  thetanum = eventnumber%(numthetas);  // % gives remainder after division
  ephinum = (eventnumber/(numthetas))%(numphis);
  zposnum = (eventnumber/(numthetas*numphis))%(numzpos);
  phicentrenum = (eventnumber/(numthetas*numphis*numzpos))%(numphicentre);

  thetafrac = ((double)thetanum)/(numthetas-1);
  ephifrac = ((double)ephinum)/(numphis-1);
  zposfrac = ((double)zposnum)/(numzpos-1);
  phicentrefrac = ((double)phicentrenum)/(numphicentre-1);

  etheta = thetamin + thetafrac*(thetamax-thetamin);
  phicentre = phicentremin + phicentrefrac*(phicentremax-phicentremin);
  phimin = phicentre - sectorhalfopeningangle*0.999;
  phimax = phicentre + sectorhalfopeningangle*0.999;
  ephi = phimin + ephifrac*(phimax-phimin);
  zpos = zposmin + zposfrac*(zposmax-zposmin);

  G4double eenergy;
    eenergy = ebeam - ebeam/(1 + mp/(2*ebeam*sin(etheta/2.)*sin(etheta/2.))); // ep elastic electrons

  G4double xdir = sin(etheta)*cos(ephi);
  G4double ydir = sin(etheta)*sin(ephi);
  G4double zdir = cos(etheta);
  printf("Primaries:%5i phi:%2i %.2f zpos:%2i %.2f; phicentre:%2i %.2f; theta:%2i %.2f =%.3f deg =%4.1f mrad  eenergy=%.3f GeV  phi=%6.1f deg, xdir=%.3f, ydir=%.3f, zdir=%.3f\n", 
       eventnumber, ephinum, ephifrac, zposnum, zposfrac, phicentrenum, phicentrefrac, 
       thetanum, thetafrac, etheta/degree, etheta/rad*1000, eenergy, ephi/degree, xdir, ydir, zdir);
  
  particleGun->SetParticleEnergy(eenergy*GeV);
  particleGun->SetParticleMomentumDirection(G4ThreeVector(xdir,ydir,zdir));
  particleGun->SetParticlePosition(G4ThreeVector(0.*cm, 0.*cm, zpos*m));
  particleGun->GeneratePrimaryVertex(anEvent);
  eventnumber++;

  gRootAnalysis->SetMomentum0(11000,0,0,11000);
  gRootAnalysis->SetMomentum2(eenergy*MeV,eenergy*MeV*xdir,eenergy*MeV*ydir,eenergy*MeV*zdir);
  gRootAnalysis->SetMomentum1(eenergy*MeV,eenergy*MeV*xdir,eenergy*MeV*ydir,eenergy*MeV*zdir);
  gRootAnalysis->SetProcess(1);

  if (eventnumber==(numphis*numthetas*numzpos*7)) {
    printf("\nFully completed one set of optics rays.\n\n");
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void MollerPrimaryGenAction::GeneratePrimaries_pion

(

G4Event *

anevent

)

Definition at line 16 of file MollerPrimaryGenAction_pion.cc.

References alpha, ebeam, gRootAnalysis, me, particleGun, pi, RootAnalysis::SetMomentum0(), RootAnalysis::SetMomentum1(), RootAnalysis::SetRate(), RootAnalysis::SetTotXS(), and wiser_sigma().

{
    const double max_pi_ang_deg = 5.0; // [deg] basic limit on what the spectrometer
                                       //       will accept in lab frame
    const double max_pi_ang = max_pi_ang_deg*pi/180.0;
    const double costhrange = 1.0 - cos(max_pi_ang);

    //Find the scattering point in the target
    double vertex_z; //scattering point in target (cm) measured from the upstream end of target
    double t = 150.0; // t = thickness of target in cm
    double rho = 0.0715; // rho = density of target in g/cm**3 (from KK)
    double radiationLength = 63.04; // g/cm**2 (from PDG)
    double t_rad = radiationLength/rho; // t_rad = radiation length in cm
    double z = G4UniformRand();
    vertex_z = z*t;//just do this for now

    double fracrl = vertex_z/t_rad; //total radiation lenght fraction travelled through
    double bt = fracrl*4./3.;

    double intrad = 2.0*log(ebeam/me)*alpha/pi;

    double vertex_x = (-0.25 + (0.5*(G4UniformRand())));
    double vertex_y = (-0.25 + (0.5*(G4UniformRand())));

    // Generate pi- angles and energy

    double pimass = 0.139;

    double pip  = G4UniformRand()*(ebeam-pimass);
    double pien = sqrt(pip*pip + pimass*pimass);
    double pith = acos( 1.0 - G4UniformRand()*costhrange );
    double piph = 2.0*pi*G4UniformRand();

    // 1e-7 is to get ub -> cm2 (10^-30) and 10^23 from Avagadro
    double factor = 6.0221415*1e-7*1.0E-6/1.602176487E-19;


    double xs = wiser_sigma(ebeam, pip, pith, bt + intrad, 1 )*1e-3; // ubarn

    // 2pi*cosrange*(Ebeam-pimass) is the phase space for such and event
    // and needs to be included
    double rate = xs*t*rho/1.008*1.0*factor*2.0*pi*costhrange*(ebeam-pimass); //Event rate/sec/ uA


    gRootAnalysis->SetTotXS(xs);
    gRootAnalysis->SetRate(rate);

    G4ParticleDefinition* particle = G4PionMinus::PionMinusDefinition();

    G4ThreeVector pipdir(sin(pith)*cos(piph), sin(pith)*sin(piph), cos(pith));

    particleGun->SetParticleDefinition(particle);
    particleGun->SetParticleMomentumDirection(pipdir);
    particleGun->SetParticleEnergy(pien*GeV);
    particleGun->SetParticlePosition(G4ThreeVector(vertex_x*cm, vertex_y*cm, vertex_z*cm));
    particleGun->GeneratePrimaryVertex(anEvent);

    gRootAnalysis->SetMomentum0(ebeam*GeV, 0.0, 0.0, 1.0);
    gRootAnalysis->SetMomentum1(pip*GeV, pipdir.x(), pipdir.y(), pipdir.z());

    anEvent->GetPrimaryVertex(0)->GetPrimary(0);

}

Here is the call graph for this function:

G4int MollerPrimaryGenAction::GetGenerator

(

)

[inline]

Definition at line 44 of file MollerPrimaryGenAction.hh.

References GeneratorNumber.

Referenced by MollerPrimaryGeneratorMessenger::GetCurrentValue().

{return GeneratorNumber;}

Here is the caller graph for this function:

double MollerPrimaryGenAction::RadProfile	(	double	eloss,
		double	btt
	)

Definition at line 95 of file MollerPrimaryGenAction.cc.

References ebeam, and me.

Referenced by EnergNumInt(), and GeneratePrimaries_elasticep_dustin().

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

Here is the caller graph for this function:

void MollerPrimaryGenAction::SetGenerator

(

G4int

val

)

[inline]

Definition at line 43 of file MollerPrimaryGenAction.hh.

References GeneratorNumber.

Referenced by MollerPrimaryGeneratorMessenger::SetNewValue().

{ GeneratorNumber = val;}

Here is the caller graph for this function:

void MollerPrimaryGenAction::setSeedValue

(

G4bool

value

)

[inline]

Definition at line 48 of file MollerPrimaryGenAction.hh.

References fSetSeedConst.

Referenced by MollerPrimaryGeneratorMessenger::SetNewValue().

                                   { fSetSeedConst = value; 
    G4cout<<"Setting fSetSeedConst"<<G4endl;};

Here is the caller graph for this function:

Double_t MollerPrimaryGenAction::wiser_all_fit	(	Double_t *	x,
		Double_t *	par
	)			[static]

Definition at line 216 of file MollerPrimaryGenAction_pion.cc.

Referenced by wiser_sigma().

                                                                        {
    // Primary variable x[0] is photon energy in [GeV]
    
    // Parameters are:
    // par[0]    Beam energy [GeV]
    Double_t Ebeam = par[0];
    // par[1]    Final particle momentum [GeV/c]
    Double_t pf    = par[1];
    // par[2]    Final particle angle [rad]
    Double_t thf   = par[2];
    // par[3]    Type (as defined in wiser_all_sig)
    Int_t type  = (Int_t) par[3];
    // par[4]    Minimum invariant mass of the residual system [GeV]
    Double_t M_X   = par[4];


    const Double_t mass_p  = 0.9383;
    const Double_t mass_p2 = mass_p*mass_p;
    const Double_t mass_pi = 0.1396;
    const Double_t mass_K  = 0.4973;

    Double_t mass[] = {mass_pi, mass_pi, mass_K, mass_K, mass_p, mass_p};

    Double_t E_gamma = x[0];

    double s = mass_p2 + 2.0*E_gamma*mass_p;

    /*  Wiser's fit    pi+     pi-    k+     k-     p+       p-  */
    Double_t A1[] =  {566.,  486.,   368., 18.2,  1.33E5,  1.63E3 };
    Double_t A2[] =  {829.,  115.,   1.91, 307.,  5.69E4, -4.30E3};
    Double_t A3[] =  {1.79,  1.77,   1.91, 0.98,  1.41,    1.79 };
    Double_t A4[] =  {2.10,  2.18,   1.15, 1.83,   .72,    2.24 };
    Double_t A6 =  1.90;
    Double_t A7 = -.0117;


    // Boost to CoM
    double beta_cm = E_gamma/(E_gamma+mass_p);
    double gamma_cm = 1.0/sqrt(1.0 - beta_cm*beta_cm);

    double p_cm_z = -gamma_cm*beta_cm*sqrt(pf*pf+mass[type]*mass[type])
              +gamma_cm*pf*cos(thf);

    double pT   = pf*sin(thf);
    double p_cm = sqrt( pT*pT + p_cm_z*p_cm_z );
    double Ef = sqrt( p_cm*p_cm + mass[type]*mass[type] );

    double p_cm_max = sqrt(s +pow(M_X*M_X - mass[type]*mass[type],2.0)/s -
     2.0*(M_X*M_X + mass[type]*mass[type]) )/2.0;

    double X_R = p_cm/p_cm_max;

    if( X_R > 1.0 ){ return 0.0; } // Kinematically forbidden


    if( type != 4 ){ // Everything but proton
  return ( A1[type] + A2[type]/sqrt(s) )*
      pow(1.0 - X_R + A3[type]*A3[type]/s, A4[type])/E_gamma;
    } else {
  Double_t U_MAN = abs(2.0*mass_p2 - 2.0*mass_p*Ef);

  return ( A1[type] + A2[type]/sqrt(s) )*
      pow(1.0 - X_R + A3[type]*A3[type]/s, A4[type])/pow(1.0 + U_MAN,A6+A7*s)/E_gamma;
    }

    return 0.0;
}

Here is the caller graph for this function:

Double_t MollerPrimaryGenAction::wiser_sigma	(	Double_t	ebeam,
		Double_t	pf,
		Double_t	thf,
		Double_t	rl,
		Int_t	type
	)

Definition at line 82 of file MollerPrimaryGenAction_pion.cc.

References wiser_all_fit().

Referenced by GeneratePrimaries_pion().

                                                                                                                   {
    /*
        Adapted from the infamous Wiser code from Peter Bosted

  Seamus Riordan
  riordan@jlab.org
  September 4, 2012
       
      type = 0 pi+
         1 pi-
         2 K+
         3 K-
         4 p
         5 p-bar
      
  type defines the type of particle produced
        Ebeam and pf are the beam energy in GeV and final particle momentum 
  in GeV/c respectively thf is the lab angle relative to the beam of 
  the produced particle in radians

  rad_len is the *effective* *total* radiator. [not in percent, typically 
  called bt in the literature]  This means put in the effective radiator for 
  internal radiation BEFORE and AFTER the vertex (total 0.05 for JLab kinematics)
        on top of real radiation lengths from ext brem times 4/3

  Returns the cross section in units of nanobarns/GeV/str
     */

    const int ntype = 6;

    Double_t A5[] = {-5.49,  -5.23, -5.91, -4.45, -6.77,  -6.53 };
    Double_t A6[] = {-1.73,  -1.82, -1.74, -3.23,  1.90,  -2.45 };

    const Double_t mass_p  = 0.9383;
    const Double_t mass_p2 = mass_p*mass_p;
    const Double_t mass_pi = 0.1396;
    const Double_t mass_K  = 0.4973;
    const Double_t mass_Lambda  = 1.116;

    Double_t mass[] = {mass_pi, mass_pi, mass_K, mass_K, mass_p, mass_p};
    Double_t *mass2 = new Double_t[ntype];

    int i;
    for( i = 0; i < 6; i++ ){
  mass2[i] = mass[i]*mass[i];
    }

    // Calculate minimum photon energy required to produce such a system
    // with the particle of transverse momentum pf*sin(thf)
    // This is the case where in the CoM frame, the particle and the
    // residual system are going back to back, transverse to the incoming particles
    // and all components of the residual system are at rest relative to one another

    double M_X;// invariant mass of the minimum residual system

    switch( type ){
  // pi+ can just be N
  case 0:
      M_X = mass_p;
      break;
  // pi- needs a pi+ N
  case 1:
      M_X = mass_p + mass_pi;
      break;
  // K+ can be created with just a Lambda, final state
  // NOTE:  This is different from the original Wiser
  // code!  They had just proton mass, but you still have
  // to have an additional strange quark lying around
  // The lowest energy configuration you can have like this
  // is the Lambda
  case 2:
      M_X = mass_Lambda;
      break;
  // p/p-bar production is twice the mass of the proton
  case 4:
  case 5:
      M_X = 2.0*mass_p;
      break;
  default:
      fprintf(stderr, "%s: %s line %d - Forbidden type passed to Wiser parameterization\n",
        __PRETTY_FUNCTION__, __FILE__, __LINE__ );
      exit(1);
      break;
    }

    double MX2  = M_X*M_X;
    double Ef = sqrt(pf*pf + mass2[type]);

    Double_t E_gamma_min = ( MX2 - mass2[type] - mass_p2 + 2.0*mass_p*Ef )/
             ( 2.0*( mass_p - Ef + pf*cos(thf) ) );

    // Parameterization parameters
    double pT = pf*sin(thf);
    double ML = sqrt(pT*pT + mass2[type]);

    double sigma, sig_e;
    double fitres;


    if( E_gamma_min > 0.0 && E_gamma_min < Ebeam ){

  TF1 *wiserfit = new TF1("wiserfit", wiser_all_fit, E_gamma_min, Ebeam, 5);
  wiserfit->SetParameter(0, Ebeam);
  wiserfit->SetParameter(1, pf);
  wiserfit->SetParameter(2, thf);
  wiserfit->SetParameter(3, (Double_t) type);
  wiserfit->SetParameter(4, M_X);

  fitres = wiserfit->Integral(E_gamma_min, Ebeam);

  delete wiserfit;

  if( type != 4 ){
      sig_e = fitres*exp(A5[type]*ML)*exp(A6[type]*pT*pT/Ef);
  } else {
      sig_e = fitres*exp(A5[type]*ML);
  }

  // Factor of 1000 is required for units
  sigma = pf*pf*sig_e*rad_len*1000.0/Ef;

  delete mass2;
  return sigma;
    } else {
  // Kinematically forbidden
  delete mass2;
  return 0.0;
    }

    delete mass2;
    return 0.0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

---

Field Documentation

const G4double MollerPrimaryGenAction::bremcut [private]

Definition at line 66 of file MollerPrimaryGenAction.hh.

Referenced by EnergNumInt(), and GeneratePrimaries_elasticep_dustin().

double_differential_t MollerPrimaryGenAction::DDtable[nbinIneEP][nbinIneEP] [private]

Definition at line 63 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries_inelasticep_dustin().

const G4double MollerPrimaryGenAction::ebeam [private]

Definition at line 64 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries_elasticep_dustin(), GeneratePrimaries_inelasticep_dustin(), GeneratePrimaries_moller_dustin(), GeneratePrimaries_phasespace_mollers(), GeneratePrimaries_phasespace_motts(), GeneratePrimaries_pion(), and RadProfile().

G4int MollerPrimaryGenAction::eventnumber [private]

Definition at line 59 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries_phasespace_mollers(), GeneratePrimaries_phasespace_motts(), and MollerPrimaryGenAction().

G4bool MollerPrimaryGenAction::fSetSeedConst [private]

Definition at line 70 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries(), and setSeedValue().

G4int MollerPrimaryGenAction::GeneratorNumber [private]

Definition at line 58 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries(), GetGenerator(), and SetGenerator().

MollerPrimaryGeneratorMessenger* MollerPrimaryGenAction::gunMessenger [private]

Definition at line 57 of file MollerPrimaryGenAction.hh.

Referenced by MollerPrimaryGenAction(), and ~MollerPrimaryGenAction().

const G4double MollerPrimaryGenAction::me [private]

Definition at line 67 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries_elasticep_dustin(), GeneratePrimaries_inelasticep_dustin(), GeneratePrimaries_moller_dustin(), GeneratePrimaries_pion(), and RadProfile().

const G4double MollerPrimaryGenAction::mp [private]

Definition at line 65 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries_elasticep_dustin(), GeneratePrimaries_inelasticep_dustin(), and GeneratePrimaries_phasespace_motts().

const G4int MollerPrimaryGenAction::nbin = 1000 [static, private]

Definition at line 60 of file MollerPrimaryGenAction.hh.

Referenced by EnergNumInt(), GeneratePrimaries_elasticep_dustin(), and GeneratePrimaries_moller_dustin().

const G4int MollerPrimaryGenAction::nbinIneEP = 100 [static, private]

Definition at line 61 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries_inelasticep_dustin().

G4ParticleGun* MollerPrimaryGenAction::particleGun [private]

Definition at line 56 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries_elasticep_dustin(), GeneratePrimaries_inelasticep_dustin(), GeneratePrimaries_moller_dustin(), GeneratePrimaries_phasespace_mollers(), GeneratePrimaries_phasespace_motts(), GeneratePrimaries_pion(), MollerPrimaryGenAction(), and ~MollerPrimaryGenAction().

cms_cross_section_t MollerPrimaryGenAction::XStable[nbin] [private]

Definition at line 62 of file MollerPrimaryGenAction.hh.

Referenced by GeneratePrimaries_elasticep_dustin(), and GeneratePrimaries_moller_dustin().

---

The documentation for this class was generated from the following files:

• MollerPrimaryGenAction.hh
• MollerPrimaryGenAction.cc
• MollerPrimaryGenAction_elasticep.cc
• MollerPrimaryGenAction_inelasticep.cc
• MollerPrimaryGenAction_moller.cc
• MollerPrimaryGenAction_phasespace.cc
• MollerPrimaryGenAction_phasespace_motts.cc
• MollerPrimaryGenAction_pion.cc