gBeamHNLEvGen.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/HNLFluxContainer.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/TuneId.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/CmdLnArgParser.h"

Include dependency graph for gBeamHNLEvGen.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	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	kDefOptFluxFilePath = "./input-flux.root"
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
bool	gOptIsMonoEnFlux = true
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	decayMod = 1.0
double	evWeight = 1.0

Function Documentation

TLorentzVector GeneratePosition

(

GHepRecord *

event

)

Definition at line 711 of file gBeamHNLEvGen.cxx.

References genie::hnl::selector::__attribute__(), fdx, fdy, fdz, fox, foy, foz, genie::RandomGen::Instance(), and genie::RandomGen::RndGeom().

Referenced by main().

{ 
  __attribute__((unused)) RandomGen * rnd = RandomGen::Instance();
  TRandom3 & rnd_geo = rnd->RndGeom();
  
  double rndx = 2 * (rnd_geo.Rndm() - 0.5); // [-1,1]
  double rndy = 2 * (rnd_geo.Rndm() - 0.5); // [-1,1]
  double rndz = 2 * (rnd_geo.Rndm() - 0.5); // [-1,1]
  
  double t_gen = 0;
  double x_gen = fox + rndx * fdx;
  double y_gen = foy + rndy * fdy;
  double z_gen = foz + rndz * fdz;
  
  TLorentzVector pos(x_gen, y_gen, z_gen, t_gen);
  return pos;
}

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

const EventRecordVisitorI * HNLGenerator

(

void

)

Definition at line 729 of file gBeamHNLEvGen.cxx.

References genie::gAbortingInErr, genie::AlgFactory::GetAlgorithm(), genie::AlgFactory::Instance(), kDefOptSConfig, kDefOptSName, LOG, pDEBUG, pFATAL, and pINFO.

Referenced by main(), and TestDecay().

{
  //string sname   = "genie::EventGenerator";
  //string sconfig = "BeamHNL";
  AlgFactory * algf = AlgFactory::Instance();

  LOG("gevgen_hnl", pINFO)
    << "Instantiating HNL generator.";

  const Algorithm * algmcgen = algf->GetAlgorithm(kDefOptSName, kDefOptSConfig);
  LOG("gevgen_hnl", pDEBUG)
    << "Got algorithm " << kDefOptSName.c_str() << "/" << kDefOptSConfig.c_str();;

  const EventRecordVisitorI * mcgen = 
    dynamic_cast< const EventRecordVisitorI * >( algmcgen );
  if(!mcgen) {
     LOG("gevgen_hnl", pFATAL) << "Couldn't instantiate the HNL generator";
     gAbortingInErr = true;
     exit(1);
  }

  LOG("gevgen_hnl", pINFO)
    << "HNL generator instantiated successfully.";

  return mcgen;
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 219 of file gBeamHNLEvGen.cxx.

References genie::hnl::selector::__attribute__(), genie::NtpWriter::AddEventRecord(), CoMLifetime, genie::NtpWriter::CustomizeFilenamePrefix(), decayMod, e, genie::NtpWriter::EventTree(), evWeight, GeneratePosition(), genie::AlgFactory::GetAlgorithm(), genie::hnl::SimpleHNL::GetCoMLifetime(), GetCommandLineArgs(), genie::hnl::Decayer::GetHNLCouplings(), genie::hnl::Decayer::GetHNLInterestingChannels(), genie::hnl::Decayer::GetHNLLifetime(), genie::hnl::FluxCreator::GetNFluxEntries(), genie::hnl::SimpleHNL::GetValidChannels(), gOptDecayMode, gOptECoupling, gOptEnergyHNL, gOptEvFilePrefix, gOptIntChannels, gOptIsMajorana, gOptIsMonoEnFlux, 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, 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, genie::utils::print::P4AsString(), pDEBUG, pINFO, pNOTICE, genie::hnl::Decayer::ProcessEventRecord(), genie::hnl::VertexGenerator::ProcessEventRecord(), genie::hnl::FluxCreator::ProcessEventRecord(), pWARN, genie::utils::app_init::RandGen(), genie::hnl::FluxCreator::RetrieveFluxInfo(), genie::RandomGen::RndGen(), genie::NtpWriter::Save(), SelectDecayMode(), genie::hnl::FluxCreator::SetFirstFluxEntry(), genie::hnl::VertexGenerator::SetGeomFile(), genie::hnl::FluxCreator::SetGeomFile(), genie::hnl::FluxCreator::SetInputFluxPath(), genie::GHepRecord::SetPrintLevel(), genie::InitialState::SetProbeP4(), genie::InitialState::SetProbePdg(), genie::GMCJMonitor::SetRefreshRate(), genie::GHepRecord::SetWeight(), 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 
  // calls LoadConfig() of each sub-alg
  __attribute__((unused)) const EventRecordVisitorI * mcgen = HNLGenerator();
  const Algorithm * algFluxCreator = AlgFactory::Instance()->GetAlgorithm("genie::hnl::FluxCreator", "Default");
  const Algorithm * algHNLGen = AlgFactory::Instance()->GetAlgorithm("genie::hnl::Decayer", "Default");
  const Algorithm * algVtxGen = AlgFactory::Instance()->GetAlgorithm("genie::hnl::VertexGenerator", "Default");

  const FluxCreator * fluxCreator = dynamic_cast< const FluxCreator * >( algFluxCreator );
  const Decayer * hnlgen = dynamic_cast< const Decayer * >( algHNLGen );
  const VertexGenerator * vtxGen = dynamic_cast< const VertexGenerator * >( algVtxGen );
  
  //string confString = kDefOptSName + "/" + kDefOptSConfig;
  string confString = kDefOptSConfig;

  CoMLifetime = hnlgen->GetHNLLifetime(); // GeV^{-1}

  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() );
  }

  //gOptIntChannels = confIntChan;

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

  // add flux info to the tree
  FluxContainer gnmf;
  if( gOptIsUsingDk2nu ) {
    // fill the flux object with nonsense to start with
    FluxContainer * ptGnmf = new FluxContainer();
    gnmf = *ptGnmf;
    delete ptGnmf;
    FillFluxNonsense( gnmf );
    TBranch * flux = ntpw.EventTree()->Branch( "flux",
                                               "genie::hnl::FluxContainer",
                                               &gnmf, 32000, 1 );
    flux->SetAutoDelete(kFALSE);
  }

  LOG("gevgen_hnl", 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");

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

  LOG("gevgen_hnl", 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__

  if( !gOptIsMonoEnFlux ){
    LOG( "gevgen_hnl", pWARN )
      << "Using input flux files. These are *flat dk2nu-like ROOT trees, so far...*";

  }
  // Event loop
  int iflux = (gOptFirstEvent < 0) ? 0 : gOptFirstEvent; int ievent = iflux;
  int maxFluxEntries = -1;
  fluxCreator->SetInputFluxPath( gOptFluxFilePath );
  fluxCreator->SetGeomFile( gOptRootGeom );
  fluxCreator->SetFirstFluxEntry( iflux );
  vtxGen->SetGeomFile( gOptRootGeom );
  
  bool tooManyEntries = false;
  while (1)
  {
    if( tooManyEntries ){
      if( gOptNev >= 10000 ){
        if( (ievent-gOptFirstEvent) % (gOptNev / 1000) == 0 ){
          int irat = (iflux-gOptFirstEvent) / ( gOptNev / 1000 );
          std::cerr << 0.1 * irat << " % " << " ( " << (iflux-gOptFirstEvent)
                    << " seen ), ( " << (ievent-gOptFirstEvent) << " / " << gOptNev  << " processed ) \r" << std::flush;
        }
      } else if( gOptNev >= 100 ) {
        if( (ievent-gOptFirstEvent) % (gOptNev / 10) == 0 ){
          int irat = (iflux-gOptFirstEvent) / ( gOptNev / 10 );
          std::cerr << 10.0 * irat << " % " << " ( " << (iflux-gOptFirstEvent)
                    << " seen ), ( " << (ievent-gOptFirstEvent) << " / " << gOptNev  << " processed ) \r" << std::flush;
        }
      }
    } else {
      if( gOptNev >= 10000 ){
        if( (ievent-gOptFirstEvent) % (gOptNev / 1000) == 0 ){
          int irat = (ievent-gOptFirstEvent) / ( gOptNev / 1000 );
          std::cerr << 0.1 * irat << " % " << " ( " << (iflux-gOptFirstEvent)
                    << " seen ), ( " << (ievent-gOptFirstEvent) << " / " << gOptNev <<  " processed ) \r" << std::flush;
        }
      } else if( gOptNev >= 100 ) {
        if( (ievent-gOptFirstEvent) % (gOptNev / 10) == 0 ){
          int irat = (ievent-gOptFirstEvent) / ( gOptNev / 10 );
          std::cerr << 10.0 * irat << " % " << " ( " << (iflux-gOptFirstEvent)
                    << " seen ), ( " << (ievent-gOptFirstEvent) << " / " << gOptNev  << " processed ) \r" << std::flush;
        }
      }
    }

    if( tooManyEntries && ((iflux-gOptFirstEvent) == gOptNev) ) break;
    else if( (ievent-gOptFirstEvent) == gOptNev ) break;
    
    if( ievent < gOptFirstEvent ){ ievent++; continue; }
    
    assert( ievent >= gOptFirstEvent && gOptFirstEvent >= 0 );
    
    LOG("gevgen_hnl", pNOTICE)
      << " *** Generating event............ " << (ievent-gOptFirstEvent);
    
    EventRecord * event = new EventRecord;
    event->SetWeight(1.0);
    event->SetProbability( CoMLifetime );
    event->SetXSec( iflux ); // will be overridden, use as handy container
    evWeight = 1.0;
    
    if( !gOptIsMonoEnFlux ){
       fluxCreator->ProcessEventRecord( event );

       // fluxCreator->ProcessEventRecord now tells us how many entries there are
       maxFluxEntries = fluxCreator->GetNFluxEntries();
       if( gOptNev > maxFluxEntries - gOptFirstEvent ){
         LOG( "gevgen_hnl", pWARN )
           << "You have asked for " << gOptNev << " events, but only provided "
           << maxFluxEntries - gOptFirstEvent << " flux entries. Truncating events to " << maxFluxEntries - gOptFirstEvent << ".";
         gOptNev = maxFluxEntries - gOptFirstEvent;
         tooManyEntries = true;
       }
       if( iflux >= maxFluxEntries - 1 ){
         LOG( "gevgen_hnl", pWARN )
           << "Reached end of flux input (iflux = " << iflux << "), about to stop.";
         break;
       }
       
       FluxContainer retGnmf = fluxCreator->RetrieveFluxInfo();
       FillFlux( gnmf, retGnmf );
       
       // check to see if this was nonsense
       if( ! event->Particle(0) ){ iflux++; delete event; continue; }
       
       gOptEnergyHNL = event->Particle(0)->GetP4()->E();
       iflux++;
     } else { // monoenergetic HNL. Add it with energy and momentum pointing on z axis
       
       assert( gOptEnergyHNL > gOptMassHNL );
       double HNLP = std::sqrt( gOptEnergyHNL*gOptEnergyHNL - gOptMassHNL*gOptMassHNL );
       TLorentzVector probeP4( 0.0, 0.0, HNLP, gOptEnergyHNL );
       TLorentzVector v4( 0.0, 0.0, 0.0, 0.0 );
       GHepParticle ptHNL( kPdgHNL, kIStInitialState, -1, -1, -1, -1, probeP4, v4 );
       event->AddParticle( ptHNL );

     }
     assert( gOptEnergyHNL > gOptMassHNL );
     
     int hpdg = genie::kPdgHNL;
     int typeMod = 1;
     if( gOptIsMajorana ) typeMod = ( rnd->RndGen().Uniform( 0.0, 1.0 ) >= 0.5 ) ? -1 : 1;
     else if( event->Particle(0)->Pdg() > 0 ) typeMod = 1;
     else if( event->Particle(0)->Pdg() < 0 ) typeMod = -1;

     // int target = SelectInitState();
     int decay  = (int) gOptDecayMode;
     
     assert( gOptECoupling >= 0.0 && gOptMCoupling >= 0.0 && gOptTCoupling >= 0.0 );
     
     // RETHERE assuming all these HNL have K+- parent. This is wrong 
     // (but not very wrong for interesting masses)
     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_hnl", pNOTICE)
         << "No decay mode specified - sampling from all available modes.";

       std::vector< HNLDecayMode_t > * intChannels = &gOptIntChannels;

       decay = SelectDecayMode( intChannels, sh );
     }

     Interaction * interaction = Interaction::HNL(typeMod * genie::kPdgHNL, gOptEnergyHNL, decay);

     if( event->Particle(0) ){ // we have an HNL with definite momentum, so let's set it now
       interaction->InitStatePtr()->SetProbeP4( *(event->Particle(0)->P4()) );
       interaction->InitStatePtr()->SetProbePdg( event->Particle(0)->Pdg() );
       LOG( "gevgen_hnl", pDEBUG )
         << "\nsetting probe p4 = " << utils::print::P4AsString( event->Particle(0)->P4() );
     }

     event->AttachSummary(interaction);

     // Simulate decay
     hnlgen->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 geometrical weight
     TLorentzVector x4mm;
     if( gOptUsingRootGeom ){
       TLorentzVector * p4HNL = event->Particle(0)->GetP4();
       NTP_IS_E = p4HNL->E(); NTP_IS_PX = p4HNL->Px(); NTP_IS_PY = p4HNL->Py(); NTP_IS_PZ = p4HNL->Pz();
       vtxGen->ProcessEventRecord( event );
       x4mm = *(event->Vertex());
     } else {
       x4mm = GeneratePosition( event );
     }

     // update weight to scale for couplings, inhibited decays
     // acceptance is already handled in FluxCreator
     // geometry handled in VertexGenerator
     evWeight = event->Weight();
     evWeight *= 1.0 / ( gOptECoupling + gOptMCoupling + gOptTCoupling );
     evWeight *= 1.0 / decayMod;
     event->SetWeight( evWeight / 1.0e+20 ); // in units of 1e+20 POT

     // store the decayMod in the second daughter's polarisation
     event->Particle(2)->SetPolarization( 1.0, decayMod );

     // why does InitState show the wrong p4 here?
     interaction->InitStatePtr()->SetProbeP4( *(event->Particle(0)->P4()) );
     
     LOG("gevgen_hnl", pDEBUG) << "Weight = " << evWeight;

     LOG("gevgen_hnl", 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_hnl", pNOTICE) << "Done!";

  return 0;
}

void PrintSyntax

(

void

)

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

Definition at line 756 of file gBeamHNLEvGen.cxx.

References genie::utils::hnl::AsString(), CoMLifetime, decayMod, genie::hnl::SimpleHNL::GetCoMLifetime(), genie::hnl::selector::GetProbabilities(), genie::hnl::SimpleHNL::GetValidChannels(), gOptMassHNL, genie::RandomGen::Instance(), genie::utils::hnl::IsKinematicallyAllowed(), LOG, pDEBUG, pFATAL, genie::RandomGen::RndGen(), genie::hnl::selector::SelectChannelInclusive(), genie::hnl::selector::SetInterestingChannels(), and genie::hnl::SimpleHNL::SetInterestingChannels().

Referenced by main().

{
  const std::map< HNLDecayMode_t, double > gammaMap = sh.GetValidChannels();

  // set CoM lifetime now if unset
  if( CoMLifetime < 0.0 ){
    CoMLifetime = sh.GetCoMLifetime();
  }

  std::vector< HNLDecayMode_t > intAndValidChannels;
  for( std::vector< HNLDecayMode_t >::iterator it = intChannels->begin(); it != intChannels->end(); ++it ){
    HNLDecayMode_t mode = *it;
    
    //check if this is a valid mode
    if( !utils::hnl::IsKinematicallyAllowed( mode, gOptMassHNL ) ) continue;

    intAndValidChannels.emplace_back( mode );
    auto mapG = gammaMap.find( mode );
    double theGamma = mapG->second;
    LOG("gevgen_hnl", pDEBUG)
      << "For mode " << utils::hnl::AsString( mode ) << " gamma = " << theGamma;
  }

  if( intAndValidChannels.size() == 0 ){ // all the modes picked by user are too heavy. Abort.
    LOG( "gevgen_hnl", pFATAL )
      << "None of the channels specified as interesting are kinematically allowed. Please either increase the HNL mass or change interesting channels in config.";
    exit(1);
  }
  
  std::map< HNLDecayMode_t, double > intMap =
    selector::SetInterestingChannels( intAndValidChannels, gammaMap );
     
  sh.SetInterestingChannels( intMap );

  // update fraction of total decay width that is not in inhibited channels
  double gammaAll = 0.0, gammaInt = 0.0;
  for( std::map< HNLDecayMode_t, double >::const_iterator itall = gammaMap.begin() ;
       itall != gammaMap.end() ; ++itall ){
    gammaAll += (*itall).second;
  }
  for( std::map< HNLDecayMode_t, double >::iterator itint = intMap.begin() ;
       itint != intMap.end() ; ++itint ){
    gammaInt += (*itint).second;
  }
  assert( gammaInt > 0.0 && gammaAll >= gammaInt );
  decayMod = gammaInt / gammaAll;

  // get probability that channels in intAndValidChannels will be selected
  std::map< HNLDecayMode_t, double > PMap = 
    selector::GetProbabilities( intMap );
     
  // now do a random throw
  RandomGen * rnd = RandomGen::Instance();
  double ranThrow = rnd->RndGen().Uniform( 0., 1. ); // HNL's fate is sealed.

  HNLDecayMode_t selectedDecayChan =
    selector::SelectChannelInclusive( PMap, ranThrow );

  int decay = ( int ) selectedDecayChan;
  return decay;
}

Variable Documentation

double CoMLifetime = -1.0

Definition at line 212 of file gBeamHNLEvGen.cxx.

Referenced by genie::hnl::Decayer::LoadConfig(), main(), ReadInConfig(), SelectDecayMode(), and TestDecay().

double decayMod = 1.0

Definition at line 214 of file gBeamHNLEvGen.cxx.

Referenced by main(), and SelectDecayMode().

double evWeight = 1.0

Definition at line 216 of file gBeamHNLEvGen.cxx.

Referenced by main().

double fdx = 0

Definition at line 198 of file gBeamHNLEvGen.cxx.

Referenced by GeneratePosition().

double fdy = 0

Definition at line 199 of file gBeamHNLEvGen.cxx.

Referenced by GeneratePosition().

double fdz = 0

Definition at line 200 of file gBeamHNLEvGen.cxx.

Referenced by GeneratePosition().

double fox = 0

Definition at line 201 of file gBeamHNLEvGen.cxx.

Referenced by GeneratePosition().

double foy = 0

Definition at line 202 of file gBeamHNLEvGen.cxx.

Referenced by GeneratePosition().

double foz = 0

Definition at line 203 of file gBeamHNLEvGen.cxx.

Referenced by GeneratePosition().

HNLDecayMode_t gOptDecayMode = kHNLDcyNull

Definition at line 181 of file gBeamHNLEvGen.cxx.

Referenced by main(), and SelectAnnihilationMode().

double gOptECoupling = -1

Definition at line 167 of file gBeamHNLEvGen.cxx.

Referenced by main().

double gOptEnergyHNL = -1

Definition at line 165 of file gBeamHNLEvGen.cxx.

Referenced by GenerateOriginMomentum(), main(), ReadInConfig(), and TestDecay().

string gOptEvFilePrefix = kDefOptEvFilePrefix

Definition at line 184 of file gBeamHNLEvGen.cxx.

double gOptGeomLUnits = 0

Definition at line 194 of file gBeamHNLEvGen.cxx.

std::vector< HNLDecayMode_t > gOptIntChannels

Definition at line 182 of file gBeamHNLEvGen.cxx.

Referenced by main().

bool gOptIsMajorana = false

Definition at line 171 of file gBeamHNLEvGen.cxx.

Referenced by main().

bool gOptIsMonoEnFlux = true

Definition at line 179 of file gBeamHNLEvGen.cxx.

Referenced by main().

double gOptMassHNL = -1

Definition at line 166 of file gBeamHNLEvGen.cxx.

Referenced by GenerateOriginMomentum(), main(), and SelectDecayMode().

double gOptMCoupling = -1

Definition at line 168 of file gBeamHNLEvGen.cxx.

Referenced by main().

int gOptNev = 10

Definition at line 163 of file gBeamHNLEvGen.cxx.

long int gOptRanSeed = -1

Definition at line 195 of file gBeamHNLEvGen.cxx.

string gOptRootGeom

Definition at line 186 of file gBeamHNLEvGen.cxx.

string gOptRootGeomTopVol = ""

Definition at line 193 of file gBeamHNLEvGen.cxx.

Long_t gOptRunNu = 1000

Definition at line 162 of file gBeamHNLEvGen.cxx.

double gOptTCoupling = -1

Definition at line 169 of file gBeamHNLEvGen.cxx.

Referenced by main().

bool gOptUsingRootGeom = false

Definition at line 185 of file gBeamHNLEvGen.cxx.

string kDefOptEvFilePrefix = "gntp"

Definition at line 154 of file gBeamHNLEvGen.cxx.

string kDefOptFluxFilePath = "./input-flux.root"

Definition at line 155 of file gBeamHNLEvGen.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 152 of file gBeamHNLEvGen.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 151 of file gBeamHNLEvGen.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 153 of file gBeamHNLEvGen.cxx.

string kDefOptSConfig = "BeamHNL"

Definition at line 158 of file gBeamHNLEvGen.cxx.

Referenced by HNLGenerator(), and main().

string kDefOptSName = "genie::EventGenerator"

Definition at line 157 of file gBeamHNLEvGen.cxx.

Referenced by HNLGenerator(), and main().

string kDefOptSTune = "GHNL20_00a_00_000"

Definition at line 159 of file gBeamHNLEvGen.cxx.

double NTP_FS0_E = 0.

Definition at line 206 of file gBeamHNLEvGen.cxx.

Referenced by main().

int NTP_FS0_PDG = 0

Definition at line 209 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS0_PX = 0.

Definition at line 206 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS0_PY = 0.

Definition at line 206 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS0_PZ = 0.

Definition at line 206 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS1_E = 0.

Definition at line 207 of file gBeamHNLEvGen.cxx.

Referenced by main().

int NTP_FS1_PDG = 0

Definition at line 209 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS1_PX = 0.

Definition at line 207 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS1_PY = 0.

Definition at line 207 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS1_PZ = 0.

Definition at line 207 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS2_E = 0.

Definition at line 208 of file gBeamHNLEvGen.cxx.

Referenced by main().

int NTP_FS2_PDG = 0

Definition at line 209 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS2_PX = 0.

Definition at line 208 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS2_PY = 0.

Definition at line 208 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_FS2_PZ = 0.

Definition at line 208 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_IS_E = 0.

Definition at line 205 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_IS_PX = 0.

Definition at line 205 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_IS_PY = 0.

Definition at line 205 of file gBeamHNLEvGen.cxx.

Referenced by main().

double NTP_IS_PZ = 0.

Definition at line 205 of file gBeamHNLEvGen.cxx.

Referenced by main().