In This Package:

GtInverseBeta.cc

Go to the documentation of this file.

#include "GtInverseBeta.h"
#include "KRLInverseBeta.hh"
#include "KRLReactorFlux.hh"
#include "GaudiKernel/IRndmGenSvc.h"

#include "HepMC/GenVertex.h"
#include "HepMC/GenEvent.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include <CLHEP/Units/PhysicalConstants.h>
#include <TROOT.h>
#include <TRandom.h>
#include <sstream>

using namespace std;
using namespace CLHEP;

Double_t inverse_beta_prob(Double_t* x, Double_t* par)
{
  Double_t engInMeV = x[0];
  Double_t cosThEL = x[1];
  Double_t prob;
  //do NOT need to multiply by the jacobian sin(theta) here!

  KRLReactorFlux gReactorFlux;
  KRLInverseBeta gInverseBeta;
  prob = gReactorFlux.Flux(engInMeV)
    *gInverseBeta.DSigDCosTh(engInMeV, cosThEL);
  //prob = gReactorFlux->Flux(engInMeV)*gInverseBeta->SigmaTot(engInMeV);
  return prob;
}


GtInverseBeta::GtInverseBeta(const std::string& type,
                             const std::string& name, 
                             const IInterface* parent)
  : GaudiTool(type,name,parent)
{
  declareInterface<IHepMCEventMutator>(this);
  
  declareProperty("NeutrinoAngle", m_neutrino_angle_in_deg=0, 
                  "neutrino angle in degree");
  declareProperty("PositronOnly",m_positron_only=0,"Only positron");
  declareProperty("NeutronOnly",m_neutron_only=0,"Only neutron");
}

GtInverseBeta::~GtInverseBeta()
{
}

HepMC::GenParticle* GtInverseBeta::make_particle(Hep3Vector& vect, 
                                                 double mass, int pid)
{
  double energy = sqrt(vect.mag2()+mass*mass);
  HepMC::FourVector fourmom(vect.x()/MeV,vect.y()/MeV,vect.z()/MeV,
                            energy/MeV);
  //the new GenParticle will be deleted by GenVertex destructor
  return new HepMC::GenParticle(fourmom,pid,1/*=status*/);
}

StatusCode GtInverseBeta::mutate(HepMC::GenEvent& event)
{
  //  gReactorFlux = new KRLReactorFlux(); // defaults to Petr's fit
  //  gInverseBeta = new KRLInverseBeta();  
    
  Hep3Vector pNeutrino; //incoming neutrino momentum
  Hep3Vector pPositron, pNeutron; // the e+ and n momentum vectors. Unit GeV/c
  
  GenerateInverseBetaPrimaries(funcInvBetaProb, m_neutrino_angle_in_deg,
                               pNeutrino, 
                               pPositron, pNeutron);
  
  //the new GenVertex will be deleted by GenEvent destructor
  HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(0,0,0,0));
  event.add_vertex(vertex);
  event.set_signal_process_vertex(vertex);
  
  int pid = 0;
  HepMC::GenParticle* particle1 = 0;
  HepMC::GenParticle* particle2 = 0;
  if(!m_positron_only && !m_neutron_only){
    pid = -11;
    particle1 = make_particle(pPositron, electron_mass_c2, pid);
    vertex->add_particle_out(particle1);
    
    pid = 2112;
    particle2 = make_particle(pNeutron, neutron_mass_c2, pid);
    vertex->add_particle_out(particle2);
  }  
  else if(m_positron_only){
    pid = -11;
    particle1 = make_particle(pPositron, electron_mass_c2, pid);
    vertex->add_particle_out(particle1);
  }
  else if(m_neutron_only){
    pid = 2112;
    particle1 = make_particle(pPositron, neutron_mass_c2, pid);
    vertex->add_particle_out(particle1);
  }
  
  int mo_pid = -12;
  HepMC::FourVector fourmom(pNeutrino.x()/MeV,pNeutrino.y()/MeV,
                              pNeutrino.z()/MeV,pNeutrino.mag()/MeV);
  HepMC::GenParticle* particle_in = 
    new HepMC::GenParticle(fourmom, mo_pid, 2);
  // Save mother particle as info particle (default status=2)
  vertex->add_particle_in(particle_in);
  
  return StatusCode::SUCCESS;
}

void GtInverseBeta::GenerateInverseBetaPrimaries(TF2* invBetaProb, 
                                                 const double nu_angle_wrt_y_in_deg,
                                                 Hep3Vector& pnu,
                                                 Hep3Vector& p1, Hep3Vector& p2)
{
  //now create a TF2 object to calculate the reaction probability
  //given the neutrino energy and the inverse beta probability
  //x: eNu in MeV
  //y: cos(th_e)
//  TF2 *invBetaProb = new TF2("InvBeta",inverse_beta_prob,
//                               0.01,60.0,-1.,1., 0); //"0" parameters
//  //set the bins on x and y to help generate pseudo random numbers
//  //according to the reaction probabilities
//  invBetaProb->SetNpx(1000);//energy
//  invBetaProb->SetNpy(100);//costh

  //generate neutrino energy and outgoing electron cos angle
  // djaffe: 12may09 modified to avoid rare case when e+ energy
  //        is returned as identically zero which indicates that below-threshold
  //        production according to KRLInverseBeta.
  Double_t aEngNu(0), aCosThEl(0);
  Double_t Ee(0);
  while ( Ee == 0 ){
    invBetaProb->GetRandom2(aEngNu, aCosThEl);
    
    //cout<<aEngNu<<" : "<<acos(aCosThEl)*360./M_2PI<<endl;
    
    //at this point, assume the incident nu momentum along z
    
    //calculate positron energy in MeV
    //    Ee = gInverseBeta->Ee1(aEngNu, aCosThEl);
    KRLInverseBeta gInverseBeta;
    Ee = gInverseBeta.Ee1(aEngNu, aCosThEl);

  }
  //generate positron phi angle
  Double_t aPhiEl = gRandom->Uniform(0.,2.0*M_PI);
  Double_t kMe = electron_mass_c2/MeV;
  Double_t Pe = sqrt(Ee*Ee-kMe*kMe);
  p1.setRThetaPhi(Pe*MeV, acos(aCosThEl), aPhiEl);
  
  //now use the 3-momentum conservation to set neutron recoil 
  //momentum
  pnu.setRThetaPhi(aEngNu*MeV,0,0);
  p2 = pnu-p1;
  
  //now make the vector rotation. first rotate around x by 90 degrees
  //so that +z axis is now +y axis
  pnu.rotateX(M_PI/2.);
  p1.rotateX(M_PI/2.);
  p2.rotateX(M_PI/2.);
  
  //if nu_angle_wrt_y is non-zero, rotate again
  if(nu_angle_wrt_y_in_deg!=0) {
    pnu.rotateX(2*M_PI/360.*nu_angle_wrt_y_in_deg);
    p1.rotateX(2*M_PI/360.*nu_angle_wrt_y_in_deg);
    p2.rotateX(2*M_PI/360.*nu_angle_wrt_y_in_deg);  
  }

  return;
}

StatusCode GtInverseBeta::initialize()
{
  static int count = 0;
  count++;
  char itoa[100];
  sprintf(itoa,"%i",count);

  string filename("InvBeta");
  filename.append(itoa);  
  funcInvBetaProb = new TF2(filename.c_str(),inverse_beta_prob,
                            0.01,60.0,-1.,1., 0); //"0" parameters
  //set the bins on x and y to help generate pseudo random numbers
  //according to the reaction probabilities
  funcInvBetaProb->SetNpx(1000);//energy
  funcInvBetaProb->SetNpy(100);//costh

  return StatusCode::SUCCESS;
}

StatusCode GtInverseBeta::finalize()
{
  delete funcInvBetaProb;

  return StatusCode::SUCCESS;
}

---