gBeamHNLParticleGun.cxx File Reference

#include <cassert>
#include <cstdlib>
#include <string>
#include <vector>
#include <sstream>
#include <TSystem.h>
#include "Framework/Algorithm/AlgFactory.h"
#include "Framework/Conventions/Controls.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/EventGen/EventGeneratorI.h"
#include "Framework/EventGen/EventRecordVisitorI.h"
#include "Framework/EventGen/GMCJMonitor.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Ntuple/NtpWriter.h"
#include "Physics/BeamHNL/HNLDecayMode.h"
#include "Physics/BeamHNL/HNLDecayUtils.h"
#include "Physics/BeamHNL/HNLVertexGenerator.h"
#include "Physics/BeamHNL/HNLFluxCreator.h"
#include "Physics/BeamHNL/HNLDecayer.h"
#include "Physics/BeamHNL/SimpleHNL.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/UnitUtils.h"
#include "Framework/Utils/PrintUtils.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/CmdLnArgParser.h"

Include dependency graph for gBeamHNLParticleGun.cxx:

Go to the source code of this file.

Functions
void	GetCommandLineArgs (int argc, char **argv)
void	PrintSyntax (void)
int	SelectDecayMode (std::vector< HNLDecayMode_t > *intChannels, SimpleHNL sh)
const EventRecordVisitorI *	HNLGenerator (void)
TLorentzVector *	GenerateOriginPosition (GHepRecord *event)
TLorentzVector *	GenerateOriginMomentum (GHepRecord *event)
TLorentzVector	GeneratePosition (GHepRecord *event)
int	main (int argc, char **argv)

Variables
string	kDefOptGeomLUnits = "mm"
string	kDefOptGeomDUnits = "g_cm3"
NtpMCFormat_t	kDefOptNtpFormat = kNFGHEP
string	kDefOptEvFilePrefix = "gntp"
string	kDefOptSName = "genie::EventGenerator"
string	kDefOptSConfig = "BeamHNL"
string	kDefOptSTune = "GHNL20_00a_00_000"
Long_t	gOptRunNu = 1000
int	gOptNev = 10
double	gOptEnergyHNL = -1
double	gOptMassHNL = -1
double	gOptECoupling = -1
double	gOptMCoupling = -1
double	gOptTCoupling = -1
bool	gOptIsMajorana = false
int	gOptHNLKind = -1
HNLDecayMode_t	gOptDecayMode = kHNLDcyNull
std::vector< HNLDecayMode_t >	gOptIntChannels
string	gOptEvFilePrefix = kDefOptEvFilePrefix
bool	gOptUsingRootGeom = false
string	gOptRootGeom
string	gOptRootGeomTopVol = ""
double	gOptGeomLUnits = 0
long int	gOptRanSeed = -1
double	fdx = 0
double	fdy = 0
double	fdz = 0
double	fox = 0
double	foy = 0
double	foz = 0
double	NTP_IS_E = 0.
double	NTP_IS_PX = 0.
double	NTP_IS_PY = 0.
double	NTP_IS_PZ = 0.
double	NTP_FS0_E = 0.
double	NTP_FS0_PX = 0.
double	NTP_FS0_PY = 0.
double	NTP_FS0_PZ = 0.
double	NTP_FS1_E = 0.
double	NTP_FS1_PX = 0.
double	NTP_FS1_PY = 0.
double	NTP_FS1_PZ = 0.
double	NTP_FS2_E = 0.
double	NTP_FS2_PX = 0.
double	NTP_FS2_PY = 0.
double	NTP_FS2_PZ = 0.
int	NTP_FS0_PDG = 0
int	NTP_FS1_PDG = 0
int	NTP_FS2_PDG = 0
double	CoMLifetime = -1.0
double	evProdVtx [3] = {0.0, 0.0, 0.0}
double	evWeight = 1.0

Function Documentation

TLorentzVector * GenerateOriginMomentum

(

GHepRecord *

event

)

Definition at line 507 of file gBeamHNLParticleGun.cxx.

References genie::AlgFactory::GetAlgorithm(), genie::hnl::Decayer::GetPGunDeviation(), genie::hnl::Decayer::GetPGunDirection(), gOptEnergyHNL, gOptMassHNL, genie::Interaction::InitStatePtr(), genie::RandomGen::Instance(), genie::AlgFactory::Instance(), genie::constants::kPi, LOG, genie::utils::print::P4AsString(), pDEBUG, genie::RandomGen::RndGen(), and genie::InitialState::SetProbeP4().

Referenced by main().

{
  double E = gOptEnergyHNL;
  double M = gOptMassHNL;
  double pMag = std::sqrt( E*E - M*M );

  const Algorithm * algHNLGen = AlgFactory::Instance()->GetAlgorithm("genie::hnl::Decayer", "Default");
  const Decayer * hnlgen = dynamic_cast< const Decayer * >( algHNLGen );

  // first map directional cosines to unit sphere
  /*
  double cx = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-cx" );
  double cy = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-cy" );
  double cz = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-cz" );
  */
  std::vector< double > PGDirection = hnlgen->GetPGunDirection();
  double cx = PGDirection.at(0), cy = PGDirection.at(1), cz = PGDirection.at(2);
  double c2 = std::sqrt( cx*cx + cy*cy + cz*cz );
  assert( c2 > 0.0 );

  cx *= 1.0 / c2; cy *= 1.0 / c2; cz *= 1.0 / c2;

  double theta = TMath::ACos( cz / c2 );
  double phi   = ( TMath::Sin( theta ) != 0.0 ) ? TMath::ACos( cx / ( c2 * TMath::Sin( theta ) ) ) : 0.0;

  // apply uniform random deviation
  std::vector< double > PGunDeviation = hnlgen->GetPGunDeviation();
  /*
  double dthetaDeg = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-DTheta" );
  double dphiDeg   = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-DPhi" );
  */
  double dthetaDeg = PGunDeviation.at(0), dphiDeg = PGunDeviation.at(1);
  double dThetaMax = dthetaDeg * constants::kPi / 180.0;
  double dPhiMax   = dphiDeg * constants::kPi / 180.0;

  RandomGen * rng = RandomGen::Instance();
  double dTheta = rng->RndGen().Uniform( -dThetaMax, dThetaMax );
  double dPhi   = rng->RndGen().Uniform( -dPhiMax, dPhiMax );
  
  std::ostringstream asts;
  asts
    << "Output details for the momentum:"
    << "\nDirectional cosines: ( " << cx << ", " << cy << ", " << cz << " )"
    << "\nMapped to ( " << theta * 180.0 / constants::kPi 
    << ", " << phi * 180.0 / constants::kPi << " ) [deg] on unit sphere"
    << "\nApplied deviation: dTheta = " << dTheta * 180.0 / constants::kPi
    << ", dPhi = " << dPhi * 180.0 / constants::kPi << " [deg]";

  // make momentum
  theta += dTheta; phi += dPhi;
  TLorentzVector * p4HNL = new TLorentzVector();
  p4HNL->SetPxPyPzE( pMag * TMath::Sin( theta ) * TMath::Cos( phi ),
                     pMag * TMath::Sin( theta ) * TMath::Sin( phi ),
                     pMag * TMath::Cos( theta ), E );
  
  asts << "\nMomentum = " << utils::print::P4AsString( p4HNL );

  Interaction * interaction = event->Summary();
  interaction->InitStatePtr()->SetProbeP4( *p4HNL );

  LOG( "gevgen_pghnl", pDEBUG ) << asts.str();

  //delete rng;
  return p4HNL;
}

TLorentzVector * GenerateOriginPosition

(

GHepRecord *

event

)

Definition at line 573 of file gBeamHNLParticleGun.cxx.

References genie::AlgFactory::GetAlgorithm(), genie::hnl::Decayer::GetPGunDOrigin(), genie::hnl::Decayer::GetPGunOrigin(), genie::RandomGen::Instance(), genie::AlgFactory::Instance(), and genie::RandomGen::RndGen().

Referenced by main().

{
  const Algorithm * algHNLGen = AlgFactory::Instance()->GetAlgorithm("genie::hnl::Decayer", "Default");
  const Decayer * hnlgen = dynamic_cast< const Decayer * >( algHNLGen );

  // get centre of box from config
  std::vector< double > PGOrigin = hnlgen->GetPGunOrigin();
  /*
  double ox = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-OriginX" );
  double oy = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-OriginY" );
  double oz = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-OriginZ" );
  */
  double ox = PGOrigin.at(0), oy = PGOrigin.at(1), oz = PGOrigin.at(2);

  // allow uniform deviation
  std::vector< double > PGDOrigin = hnlgen->GetPGunDOrigin();
  /*
  double Dxmax = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-OriginDX" );
  double Dymax = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-OriginDY" );
  double Dzmax = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-OriginDZ" );
  */
  double Dxmax = PGDOrigin.at(0), Dymax = PGDOrigin.at(1), Dzmax = PGDOrigin.at(2);

  RandomGen * rng = RandomGen::Instance();
  double dx = rng->RndGen().Uniform( -Dxmax, Dxmax );
  double dy = rng->RndGen().Uniform( -Dymax, Dymax );
  double dz = rng->RndGen().Uniform( -Dzmax, Dzmax );

  // make position
  ox += dx; oy += dy; oz += dz;
  TLorentzVector * x4HNL = new TLorentzVector();
  x4HNL->SetXYZT( ox, oy, oz, 0.0 );

  event->SetVertex( *x4HNL );
  //delete rng;
  return x4HNL;
}

TLorentzVector GeneratePosition

(

GHepRecord *

event

)

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

const EventRecordVisitorI* HNLGenerator

(

void

)

int main	(	int	argc,
		char **	argv
	)

Definition at line 194 of file gBeamHNLParticleGun.cxx.

References genie::hnl::selector::__attribute__(), genie::NtpWriter::AddEventRecord(), genie::utils::hnl::AsString(), CoMLifetime, genie::NtpWriter::CustomizeFilenamePrefix(), genie::NtpWriter::EventTree(), evWeight, GenerateOriginMomentum(), GenerateOriginPosition(), GeneratePosition(), genie::AlgFactory::GetAlgorithm(), genie::hnl::SimpleHNL::GetCoMLifetime(), GetCommandLineArgs(), genie::hnl::Decayer::GetHNLCouplings(), genie::hnl::Decayer::GetHNLInterestingChannels(), genie::hnl::Decayer::GetHNLLifetime(), genie::hnl::Decayer::GetHNLMass(), genie::hnl::Decayer::GetPGunEnergy(), genie::hnl::SimpleHNL::GetValidChannels(), gOptDecayMode, gOptECoupling, gOptEnergyHNL, gOptEvFilePrefix, gOptIntChannels, gOptIsMajorana, gOptMassHNL, gOptMCoupling, gOptNev, gOptRanSeed, gOptRootGeom, gOptRunNu, gOptTCoupling, gOptUsingRootGeom, genie::Interaction::HNL(), HNLGenerator(), genie::NtpWriter::Initialize(), genie::Interaction::InitStatePtr(), genie::RandomGen::Instance(), genie::RunOpt::Instance(), genie::AlgFactory::Instance(), genie::hnl::Decayer::IsHNLMajorana(), kDefOptNtpFormat, kDefOptSConfig, kDefOptSName, genie::hnl::kHNLDcyNull, genie::kIStInitialState, genie::kPdgHNL, genie::kPdgKP, LOG, genie::utils::app_init::MesgThresholds(), NTP_FS0_E, NTP_FS0_PDG, NTP_FS0_PX, NTP_FS0_PY, NTP_FS0_PZ, NTP_FS1_E, NTP_FS1_PDG, NTP_FS1_PX, NTP_FS1_PY, NTP_FS1_PZ, NTP_FS2_E, NTP_FS2_PDG, NTP_FS2_PX, NTP_FS2_PY, NTP_FS2_PZ, NTP_IS_E, NTP_IS_PX, NTP_IS_PY, NTP_IS_PZ, pDEBUG, pFATAL, pINFO, pNOTICE, genie::hnl::Decayer::ProcessEventRecord(), genie::hnl::VertexGenerator::ProcessEventRecord(), genie::utils::app_init::RandGen(), genie::NtpWriter::Save(), SelectDecayMode(), genie::hnl::VertexGenerator::SetGeomFile(), genie::GHepRecord::SetPrintLevel(), genie::GHepRecord::SetProbability(), genie::InitialState::SetProbeP4(), genie::GMCJMonitor::SetRefreshRate(), and genie::GMCJMonitor::Update().

{
  // suppress ROOT Info messages
  gErrorIgnoreLevel = kWarning;

  // Parse command line arguments
  GetCommandLineArgs(argc,argv);

  // Init messenger and random number seed
  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
  utils::app_init::RandGen(gOptRanSeed);

  __attribute__((unused)) RandomGen * rnd = RandomGen::Instance();

  // Get the HNL generator first to load config
  // config loaded upon instantiation of HNLGenerator algorithm 
  // ==> Decayer::LoadConfig()
  __attribute__((unused)) const EventRecordVisitorI * mcgen = HNLGenerator();
  const Algorithm * algHNLGen = AlgFactory::Instance()->GetAlgorithm("genie::hnl::Decayer", "Default");
  const Algorithm * algVtxGen = AlgFactory::Instance()->GetAlgorithm("genie::hnl::VertexGenerator", "Default");
  
  const Decayer * hnlgen = dynamic_cast< const Decayer * >( algHNLGen );
  const VertexGenerator * vtxGen = dynamic_cast< const VertexGenerator * >( algVtxGen );

  bool geom_is_accessible = ! (gSystem->AccessPathName(gOptRootGeom.c_str()));
  if (!geom_is_accessible) {
    LOG("gevgen_pghnl", pFATAL)
      << "The specified ROOT geometry doesn't exist! Initialization failed!";
    exit(1);
  }
  vtxGen->SetGeomFile( gOptRootGeom );

  if( !gOptRootGeoManager ) gOptRootGeoManager = TGeoManager::Import(gOptRootGeom.c_str()); 

  TGeoVolume * top_volume = gOptRootGeoManager->GetTopVolume();
  assert( top_volume );
  __attribute__((unused)) TGeoShape * ts  = top_volume->GetShape();

  //TGeoBBox *  box = (TGeoBBox *)ts;

  string confString = kDefOptSName + kDefOptSConfig;

  CoMLifetime = hnlgen->GetHNLLifetime();

  gOptMassHNL    = hnlgen->GetHNLMass();
  std::vector< double > U4l2s = hnlgen->GetHNLCouplings();
  gOptECoupling  = U4l2s.at(0);
  gOptMCoupling  = U4l2s.at(1);
  gOptTCoupling  = U4l2s.at(2);
  gOptIsMajorana = hnlgen->IsHNLMajorana();

  std::string stIntChannels = hnlgen->GetHNLInterestingChannels(); int iChan = -1;
  if( gOptIntChannels.size() > 0 ) gOptIntChannels.clear();
  while( stIntChannels.size() > 0 ){ // read channels from right (lowest mass) to left (highest mass)
    iChan++;
    HNLDecayMode_t md = static_cast< HNLDecayMode_t >( iChan );
    std::string tmpSt = stIntChannels.substr( stIntChannels.size()-1, stIntChannels.size() );
    if( std::strcmp( tmpSt.c_str(), "1" ) == 0 )
      gOptIntChannels.emplace_back( md );

    stIntChannels.erase( stIntChannels.end()-1, stIntChannels.end() );
  }

  assert( gOptECoupling >= 0.0 && gOptMCoupling >= 0.0 && gOptTCoupling >= 0.0 );

  // Initialize an Ntuple Writer to save GHEP records into a TTree
  NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu, gOptRanSeed);
  ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix);
  ntpw.Initialize();

  LOG("gevgen_pghnl", pNOTICE)
    << "Initialised Ntuple Writer";

  // add another few branches to tree.
  ntpw.EventTree()->Branch("hnl_mass", &gOptMassHNL, "gOptMassHNL/D");
  ntpw.EventTree()->Branch("hnl_coup_e", &gOptECoupling, "gOptECoupling/D");
  ntpw.EventTree()->Branch("hnl_coup_m", &gOptMCoupling, "gOptMCoupling/D");
  ntpw.EventTree()->Branch("hnl_coup_t", &gOptTCoupling, "gOptTCoupling/D");
  ntpw.EventTree()->Branch("hnl_ismaj", &gOptIsMajorana, "gOptIsMajorana/I");

  // let's make HNL-specific FS branches until we get gntpc sorted out
  ntpw.EventTree()->Branch("hnl_IS_E", &NTP_IS_E, "NTP_IS_E/D");
  ntpw.EventTree()->Branch("hnl_IS_PX", &NTP_IS_PX, "NTP_IS_PX/D");
  ntpw.EventTree()->Branch("hnl_IS_PY", &NTP_IS_PY, "NTP_IS_PY/D");
  ntpw.EventTree()->Branch("hnl_IS_PZ", &NTP_IS_PZ, "NTP_IS_PZ/D");
  ntpw.EventTree()->Branch("hnl_FS0_PDG", &NTP_FS0_PDG, "NTP_FS0_PDG/I");
  ntpw.EventTree()->Branch("hnl_FS0_E", &NTP_FS0_E, "NTP_FS0_E/D");
  ntpw.EventTree()->Branch("hnl_FS0_PX", &NTP_FS0_PX, "NTP_FS0_PX/D");
  ntpw.EventTree()->Branch("hnl_FS0_PY", &NTP_FS0_PY, "NTP_FS0_PY/D");
  ntpw.EventTree()->Branch("hnl_FS0_PZ", &NTP_FS0_PZ, "NTP_FS0_PZ/D");
  ntpw.EventTree()->Branch("hnl_FS1_PDG", &NTP_FS1_PDG, "NTP_FS1_PDG/I");
  ntpw.EventTree()->Branch("hnl_FS1_E", &NTP_FS1_E, "NTP_FS1_E/D");
  ntpw.EventTree()->Branch("hnl_FS1_PX", &NTP_FS1_PX, "NTP_FS1_PX/D");
  ntpw.EventTree()->Branch("hnl_FS1_PY", &NTP_FS1_PY, "NTP_FS1_PY/D");
  ntpw.EventTree()->Branch("hnl_FS1_PZ", &NTP_FS1_PZ, "NTP_FS1_PZ/D");
  ntpw.EventTree()->Branch("hnl_FS2_PDG", &NTP_FS2_PDG, "NTP_FS2_PDG/I");
  ntpw.EventTree()->Branch("hnl_FS2_E", &NTP_FS2_E, "NTP_FS2_E/D");
  ntpw.EventTree()->Branch("hnl_FS2_PX", &NTP_FS2_PX, "NTP_FS2_PX/D");
  ntpw.EventTree()->Branch("hnl_FS2_PY", &NTP_FS2_PY, "NTP_FS2_PY/D");
  ntpw.EventTree()->Branch("hnl_FS2_PZ", &NTP_FS2_PZ, "NTP_FS2_PZ/D");

  // Create a MC job monitor for a periodically updated status file
  GMCJMonitor mcjmonitor(gOptRunNu);
  mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());

  LOG("gevgen_pghnl", pNOTICE)
    << "Initialised MC job monitor";

  // Set GHEP print level
  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());

#ifdef __CAN_USE_ROOT_GEOM__
  // Read geometry bounding box - for vertex position generation
  if( gOptUsingRootGeom ){
    InitBoundingBox();
  }
#endif // #ifdef __CAN_USE_ROOT_GEOM__

  // read in energy. This is always constant
  //gOptEnergyHNL = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-Energy" );
  gOptEnergyHNL = hnlgen->GetPGunEnergy();
  if( gOptEnergyHNL <= gOptMassHNL ){
    LOG( "gevgen_pghnl", pFATAL )
      << "\nInsufficient energy " << gOptEnergyHNL << " GeV for HNL of mass " << gOptMassHNL << " GeV. Exiting."
      << "\nPlease check ${GENIE}/config/CommonHNL.xml sections \"ParamSpace\" and \"ParticleGun\"";
    exit(1);
  }
  assert( gOptEnergyHNL > gOptMassHNL );

  // Event loop
  int ievent = 0;

  while (1)
  {
    if( gOptNev >= 10000 ){
      if( ievent % (gOptNev / 1000 ) == 0 ){
        int irat = ievent / ( gOptNev / 1000 );
        std::cerr << 0.1 * irat << " % " << " ( " << ievent
                  << " / " << gOptNev << " ) \r" << std::flush;
      }
    } else if( gOptNev >= 100 ){
      if( ievent % (gOptNev / 10 ) == 0 ){
        int irat = ievent / ( gOptNev / 10 );
        std::cerr << 10.0 * irat << " % " << " ( " << ievent
                  << " / " << gOptNev << " ) \r" << std::flush;
      }
    }

    if((ievent) == gOptNev) break;
      
     LOG("gevgen_pghnl", pNOTICE)
          << " *** Generating event............ " << ievent;

     int hpdg = genie::kPdgHNL;
     EventRecord * event = new EventRecord;

     event->SetProbability( CoMLifetime );

     int decay  = (int) gOptDecayMode;

     SimpleHNL sh( "HNL", ievent, hpdg, genie::kPdgKP, 
                   gOptMassHNL, gOptECoupling, gOptMCoupling, gOptTCoupling, false );
     
     const std::map< HNLDecayMode_t, double > gammaMap = sh.GetValidChannels();
     CoMLifetime = sh.GetCoMLifetime();
     
     if( gOptDecayMode == kHNLDcyNull ){ // select from available modes
       LOG("gevgen_pghnl", pNOTICE)
         << "No decay mode specified - sampling from all available modes.";

       std::vector< HNLDecayMode_t > * intChannels = &gOptIntChannels;
       
       decay = SelectDecayMode( intChannels, sh );
     }

     Interaction * interaction = Interaction::HNL(genie::kPdgHNL, gOptEnergyHNL, decay);
     event->AttachSummary(interaction);

     // generate origin position and momentum direction for this HNL
     TLorentzVector * x4orig = GenerateOriginPosition( event );
     TLorentzVector * p4HNL  = GenerateOriginMomentum( event );

     // and add as Particle(0)
     GHepParticle ptHNL( kPdgHNL, kIStInitialState, -1, -1, -1, -1, *p4HNL, *x4orig );
     event->AddParticle( ptHNL );

     LOG("gevgen_pghnl", pDEBUG)
       << "Note decay mode is " << utils::hnl::AsString(gOptDecayMode);

     // Simulate decay
     hnlgen->ProcessEventRecord(event);
     vtxGen->ProcessEventRecord(event);

     // add the FS 4-momenta to special branches
     // Quite inelegant. Gets the job done, though
     NTP_FS0_PDG = (event->Particle(1))->Pdg();
     NTP_FS0_E  = ((event->Particle(1))->P4())->E();
     NTP_FS0_PX = ((event->Particle(1))->P4())->Px();
     NTP_FS0_PY = ((event->Particle(1))->P4())->Py();
     NTP_FS0_PZ = ((event->Particle(1))->P4())->Pz();
     NTP_FS1_PDG = (event->Particle(2))->Pdg();
     NTP_FS1_E  = ((event->Particle(2))->P4())->E();
     NTP_FS1_PX = ((event->Particle(2))->P4())->Px();
     NTP_FS1_PY = ((event->Particle(2))->P4())->Py();
     NTP_FS1_PZ = ((event->Particle(2))->P4())->Pz();
     if( event->Particle(3) ){
       NTP_FS2_PDG = (event->Particle(3))->Pdg();
       NTP_FS2_E  = ((event->Particle(3))->P4())->E();
       NTP_FS2_PX = ((event->Particle(3))->P4())->Px();
       NTP_FS2_PY = ((event->Particle(3))->P4())->Py();
       NTP_FS2_PZ = ((event->Particle(3))->P4())->Pz();
     }
     else{
       NTP_FS2_PDG = 0;
       NTP_FS2_E = 0.0;
       NTP_FS2_PX = 0.0;
       NTP_FS2_PY = 0.0;
       NTP_FS2_PZ = 0.0;
     }

     // Generate (or read) a position for the decay vertex
     // also currently handles the event weight
     TLorentzVector x4mm = GeneratePosition( event );
     event->SetWeight( evWeight );

     // why does InitState show the wrong p4 here?
     interaction->InitStatePtr()->SetProbeP4( *(event->Particle(0)->P4()) );

     LOG("gevgen_pghnl", pINFO)
         << "Generated event: " << *event;

     // Add event at the output ntuple, refresh the mc job monitor & clean-up
     ntpw.AddEventRecord(ievent, event);
     mcjmonitor.Update(ievent,event);
     delete event;

     ievent++;
  } // event loop

  // Save the generated event tree & close the output file
  ntpw.Save();

  LOG("gevgen_pghnl", pNOTICE) << "Done!";

  return 0;
}

void PrintSyntax

(

void

)

int SelectDecayMode	(	std::vector< HNLDecayMode_t > *	intChannels,
		SimpleHNL	sh
	)

Variable Documentation

double CoMLifetime = -1.0

Definition at line 187 of file gBeamHNLParticleGun.cxx.

double evProdVtx[3] = {0.0, 0.0, 0.0}

Definition at line 189 of file gBeamHNLParticleGun.cxx.

double evWeight = 1.0

Definition at line 191 of file gBeamHNLParticleGun.cxx.

double fdx = 0

Definition at line 172 of file gBeamHNLParticleGun.cxx.

double fdy = 0

Definition at line 173 of file gBeamHNLParticleGun.cxx.

double fdz = 0

Definition at line 174 of file gBeamHNLParticleGun.cxx.

double fox = 0

Definition at line 175 of file gBeamHNLParticleGun.cxx.

double foy = 0

Definition at line 176 of file gBeamHNLParticleGun.cxx.

double foz = 0

Definition at line 177 of file gBeamHNLParticleGun.cxx.

HNLDecayMode_t gOptDecayMode = kHNLDcyNull

Definition at line 155 of file gBeamHNLParticleGun.cxx.

double gOptECoupling = -1

Definition at line 148 of file gBeamHNLParticleGun.cxx.

double gOptEnergyHNL = -1

Definition at line 146 of file gBeamHNLParticleGun.cxx.

string gOptEvFilePrefix = kDefOptEvFilePrefix

Definition at line 158 of file gBeamHNLParticleGun.cxx.

double gOptGeomLUnits = 0

Definition at line 168 of file gBeamHNLParticleGun.cxx.

int gOptHNLKind = -1

Definition at line 153 of file gBeamHNLParticleGun.cxx.

std::vector< HNLDecayMode_t > gOptIntChannels

Definition at line 156 of file gBeamHNLParticleGun.cxx.

bool gOptIsMajorana = false

Definition at line 152 of file gBeamHNLParticleGun.cxx.

double gOptMassHNL = -1

Definition at line 147 of file gBeamHNLParticleGun.cxx.

double gOptMCoupling = -1

Definition at line 149 of file gBeamHNLParticleGun.cxx.

int gOptNev = 10

Definition at line 144 of file gBeamHNLParticleGun.cxx.

long int gOptRanSeed = -1

Definition at line 169 of file gBeamHNLParticleGun.cxx.

string gOptRootGeom

Definition at line 160 of file gBeamHNLParticleGun.cxx.

string gOptRootGeomTopVol = ""

Definition at line 167 of file gBeamHNLParticleGun.cxx.

Long_t gOptRunNu = 1000

Definition at line 143 of file gBeamHNLParticleGun.cxx.

double gOptTCoupling = -1

Definition at line 150 of file gBeamHNLParticleGun.cxx.

bool gOptUsingRootGeom = false

Definition at line 159 of file gBeamHNLParticleGun.cxx.

string kDefOptEvFilePrefix = "gntp"

Definition at line 136 of file gBeamHNLParticleGun.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 134 of file gBeamHNLParticleGun.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 133 of file gBeamHNLParticleGun.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 135 of file gBeamHNLParticleGun.cxx.

string kDefOptSConfig = "BeamHNL"

Definition at line 139 of file gBeamHNLParticleGun.cxx.

string kDefOptSName = "genie::EventGenerator"

Definition at line 138 of file gBeamHNLParticleGun.cxx.

string kDefOptSTune = "GHNL20_00a_00_000"

Definition at line 140 of file gBeamHNLParticleGun.cxx.

double NTP_FS0_E = 0.

Definition at line 180 of file gBeamHNLParticleGun.cxx.

int NTP_FS0_PDG = 0

Definition at line 183 of file gBeamHNLParticleGun.cxx.

double NTP_FS0_PX = 0.

Definition at line 180 of file gBeamHNLParticleGun.cxx.

double NTP_FS0_PY = 0.

Definition at line 180 of file gBeamHNLParticleGun.cxx.

double NTP_FS0_PZ = 0.

Definition at line 180 of file gBeamHNLParticleGun.cxx.

double NTP_FS1_E = 0.

Definition at line 181 of file gBeamHNLParticleGun.cxx.

int NTP_FS1_PDG = 0

Definition at line 183 of file gBeamHNLParticleGun.cxx.

double NTP_FS1_PX = 0.

Definition at line 181 of file gBeamHNLParticleGun.cxx.

double NTP_FS1_PY = 0.

Definition at line 181 of file gBeamHNLParticleGun.cxx.

double NTP_FS1_PZ = 0.

Definition at line 181 of file gBeamHNLParticleGun.cxx.

double NTP_FS2_E = 0.

Definition at line 182 of file gBeamHNLParticleGun.cxx.

int NTP_FS2_PDG = 0

Definition at line 183 of file gBeamHNLParticleGun.cxx.

double NTP_FS2_PX = 0.

Definition at line 182 of file gBeamHNLParticleGun.cxx.

double NTP_FS2_PY = 0.

Definition at line 182 of file gBeamHNLParticleGun.cxx.

double NTP_FS2_PZ = 0.

Definition at line 182 of file gBeamHNLParticleGun.cxx.

double NTP_IS_E = 0.

Definition at line 179 of file gBeamHNLParticleGun.cxx.

double NTP_IS_PX = 0.

Definition at line 179 of file gBeamHNLParticleGun.cxx.

double NTP_IS_PY = 0.

Definition at line 179 of file gBeamHNLParticleGun.cxx.

double NTP_IS_PZ = 0.

Definition at line 179 of file gBeamHNLParticleGun.cxx.