gBeamHNLValidationApp.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/HNLFluxContainer.h"
#include "Physics/BeamHNL/HNLDecayer.h"
#include "Physics/BeamHNL/HNLProductionMode.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 gBeamHNLValidationApp.cxx:

Go to the source code of this file.

Typedefs
typedef enum t_HNLValidation	HNLValidation_t

Enumerations
enum	t_HNLValidation {   kValNone = 0, kValFluxFromDk2nu = 1, kValDecay = 2, kValGeom = 3,   kValFullSim = 4 }

Functions
void	GetCommandLineArgs (int argc, char **argv)
void	ReadInConfig (void)
void	PrintSyntax (void)
const EventRecordVisitorI *	HNLGenerator (void)
int	TestDecay (void)
int	main (int argc, char **argv)

Variables
string	kDefOptGeomLUnits = "mm"
string	kDefOptGeomAUnits = "rad"
string	kDefOptGeomTUnits = "ns"
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
HNLValidation_t	gOptValidationMode = kValNone
double	gCfgMassHNL = -1
double	gCfgECoupling = -1
double	gCfgMCoupling = -1
double	gCfgTCoupling = -1
bool	gCfgIsMajorana = false
int	gCfgHNLKind = 2
double	CoMLifetime = -1
HNLProd_t	gCfgProdMode = kHNLProdNull
HNLDecayMode_t	gCfgDecayMode = kHNLDcyNull
std::vector< HNLDecayMode_t >	gCfgIntChannels
double	gCfgUserOx = -1
double	gCfgUserOy = -1
double	gCfgUserOz = -1
double	gCfgUserAx1 = -1
double	gCfgUserAz = -1
double	gCfgUserAx2 = -1
double	gCfgBoxLx = -1
double	gCfgBoxLy = -1
double	gCfgBoxLz = -1
double	gCfgParentEnergy = -1
double	gCfgParentTheta = -1
double	gCfgParentPhi = -1
double	gCfgParentOx = -1
double	gCfgParentOy = -1
double	gCfgParentOz = -1
double	gCfgHNLCx = -1
double	gCfgHNLCy = -1
double	gCfgHNLCz = -1
double	gOptEnergyHNL = -1
bool	gOptIsMonoEnFlux = true
string	gOptEvFilePrefix = kDefOptEvFilePrefix
bool	gOptUsingRootGeom = false
string	gOptRootGeom
string	geom
string	lunits
string	aunits
string	tunits
double	gOptGeomLUnits = 0
double	gOptGeomAUnits = 0
double	gOptGeomTUnits = 0
string	gOptRootGeomTopVol = ""
double	fdx = 0
double	fdy = 0
double	fdz = 0
double	fox = 0
double	foy = 0
double	foz = 0
long int	gOptRanSeed = -1

Typedef Documentation

typedef enum t_HNLValidation HNLValidation_t

Enumeration Type Documentation

enum t_HNLValidation

Enumerator
kValNone
kValFluxFromDk2nu
kValDecay
kValGeom
kValFullSim

Definition at line 134 of file gBeamHNLValidationApp.cxx.

                             {
  
  kValNone            = 0,
  kValFluxFromDk2nu   = 1,
  kValDecay           = 2,
  kValGeom            = 3,
  kValFullSim         = 4
  
} HNLValidation_t;

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

const EventRecordVisitorI* HNLGenerator

(

void

)

int main	(	int	argc,
		char **	argv
	)

Definition at line 255 of file gBeamHNLValidationApp.cxx.

References GetCommandLineArgs(), gOptRanSeed, gOptValidationMode, genie::RunOpt::Instance(), kValDecay, kValFluxFromDk2nu, kValGeom, LOG, genie::utils::app_init::MesgThresholds(), pFATAL, genie::utils::app_init::RandGen(), ReadInConfig(), and TestDecay().

{
  // suppress ROOT Info messages
  gErrorIgnoreLevel = kWarning;

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

  // Get the validation configuration
  ReadInConfig();

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

  switch( gOptValidationMode ){
    
  case kValFluxFromDk2nu: return TestFluxFromDk2nu(); break;
  case kValDecay:         return TestDecay();         break;
  case kValGeom:          return TestGeom();          break;
  default: LOG( "gevald_hnl", pFATAL ) << "I didn't recognise this mode. Goodbye world!"; break;
  }

  return 0;
}

void PrintSyntax

(

void

)

void ReadInConfig

(

void

)

Definition at line 1421 of file gBeamHNLValidationApp.cxx.

References genie::hnl::selector::__attribute__(), genie::utils::hnl::AsString(), aunits, CoMLifetime, gCfgECoupling, gCfgHNLCx, gCfgHNLCy, gCfgHNLCz, gCfgIntChannels, gCfgIsMajorana, gCfgMassHNL, gCfgMCoupling, gCfgTCoupling, gCfgUserAx1, gCfgUserAx2, gCfgUserAz, gCfgUserOx, gCfgUserOy, gCfgUserOz, genie::AlgFactory::GetAlgorithm(), genie::hnl::FluxCreator::GetB2URotation(), genie::hnl::FluxCreator::GetB2UTranslation(), genie::hnl::Decayer::GetHNLCouplings(), genie::hnl::Decayer::GetHNLInterestingChannels(), genie::hnl::Decayer::GetHNLLifetime(), genie::hnl::Decayer::GetHNLMass(), genie::hnl::Decayer::GetPGunDirection(), genie::hnl::Decayer::GetPGunEnergy(), gOptEnergyHNL, gOptGeomAUnits, gOptGeomLUnits, genie::AlgFactory::Instance(), genie::hnl::Decayer::IsHNLMajorana(), LOG, lunits, genie::units::m, pDEBUG, pFATAL, and genie::units::rad.

Referenced by main().

{
  LOG("gevald_hnl", pFATAL)
    << "Reading in validation configuration. . .";

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

  CoMLifetime = hnlgen->GetHNLLifetime();

  gCfgMassHNL    = hnlgen->GetHNLMass();
  std::vector<double> U4l2s = hnlgen->GetHNLCouplings();
  gCfgECoupling  = U4l2s.at(0);
  gCfgMCoupling  = U4l2s.at(1);
  gCfgTCoupling  = U4l2s.at(2);
  gCfgIsMajorana = hnlgen->IsHNLMajorana();

  //gOptEnergyHNL = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-Energy" );
  gOptEnergyHNL = hnlgen->GetPGunEnergy();

  std::vector< double > PGDirection = hnlgen->GetPGunDirection();
  gCfgHNLCx = PGDirection.at(0), gCfgHNLCy = PGDirection.at(1), gCfgHNLCz = PGDirection.at(2);
  /*
  gCfgHNLCx     = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-cx" );
  gCfgHNLCy     = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-cy" );
  gCfgHNLCz     = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-cz" );
  */

  double dircosMag2 = std::pow( gCfgHNLCx, 2.0 ) + 
    std::pow( gCfgHNLCy, 2.0 ) + 
    std::pow( gCfgHNLCz, 2.0 );
  double invdircosmag = 1.0 / std::sqrt( dircosMag2 );
  gCfgHNLCx *= invdircosmag;
  gCfgHNLCy *= invdircosmag;
  gCfgHNLCz *= invdircosmag;
  
  std::string stIntChannels = hnlgen->GetHNLInterestingChannels(); int iChan = -1;
  if( gCfgIntChannels.size() > 0 ) gCfgIntChannels.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 )
      gCfgIntChannels.emplace_back( md );

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

  const Algorithm * algFluxCreator = AlgFactory::Instance()->GetAlgorithm("genie::hnl::FluxCreator", "Default");
  __attribute__((unused)) const FluxCreator * fluxCreator = dynamic_cast< const FluxCreator * >( algFluxCreator );

  std::vector< double > b2uTranslation = fluxCreator->GetB2UTranslation();
  gCfgUserOx = b2uTranslation.at(0); gCfgUserOy = b2uTranslation.at(1); gCfgUserOz = b2uTranslation.at(2);
  std::vector< double > b2uRotation = fluxCreator->GetB2URotation();
  gCfgUserAx1 = b2uRotation.at(0); gCfgUserAz = b2uRotation.at(1); gCfgUserAx2 = b2uRotation.at(2);

  // now transform the lengths and angles to the correct units
  gCfgUserOx   *= units::m / gOptGeomLUnits;
  gCfgUserOy   *= units::m / gOptGeomLUnits;
  gCfgUserOz   *= units::m / gOptGeomLUnits;

  gCfgUserAx1  *= units::rad / gOptGeomAUnits;
  gCfgUserAz   *= units::rad / gOptGeomAUnits;
  gCfgUserAx2  *= units::rad / gOptGeomAUnits;

  ostringstream csts; 
  csts << "Read out the following config:"
       << "\n"
       << "\nHNL mass = " << gCfgMassHNL << " [GeV]"
       << "\n|U_e4|^2 = " << gCfgECoupling
       << "\n|U_m4|^2 = " << gCfgMCoupling
       << "\n|U_t4|^2 = " << gCfgTCoupling
       << "\n"
       << "\nInteresting decay channels:";
  for( std::vector< HNLDecayMode_t >::iterator chit = gCfgIntChannels.begin();
       chit != gCfgIntChannels.end(); ++chit ){ csts << "\n\t" << utils::hnl::AsString(*chit); }
  csts << "\n"
       << "\nUser origin in beam coordinates = ( " << gCfgUserOx
       << ", " << gCfgUserOy << ", " << gCfgUserOz << " ) [" << lunits.c_str() << "]"
       << "\nEuler extrinsic x-z-x rotation = ( " << gCfgUserAx1
       << ", " << gCfgUserAz << ", " << gCfgUserAx2 << " ) [" << aunits.c_str() << "]"
       << "\nHNL particle-gun directional cosines: ( " << gCfgHNLCx << ", " << gCfgHNLCy 
       << ", " << gCfgHNLCz << ") [ GeV / GeV ]";

  LOG("gevald_hnl", pDEBUG) << csts.str();
  
}

int TestDecay

(

void

)

Definition at line 545 of file gBeamHNLValidationApp.cxx.

References genie::hnl::selector::__attribute__(), genie::utils::hnl::AsString(), CoMLifetime, genie::XclsTag::DecayMode(), genie::Interaction::ExclTag(), gCfgECoupling, gCfgHNLCx, gCfgHNLCy, gCfgHNLCz, gCfgMassHNL, gCfgMCoupling, gCfgTCoupling, genie::AlgFactory::GetAlgorithm(), genie::hnl::Decayer::GetPGunEnergy(), genie::hnl::Decayer::GetPGunOrigin(), genie::hnl::SimpleHNL::GetValidChannels(), gOptEnergyHNL, gOptNev, gOptRootGeom, genie::Interaction::HNL(), HNLGenerator(), genie::Interaction::InitStatePtr(), genie::RunOpt::Instance(), genie::AlgFactory::Instance(), genie::hnl::kHNLDcyNuEE, genie::hnl::kHNLDcyNull, genie::hnl::kHNLDcyNuMuE, genie::hnl::kHNLDcyNuMuMu, genie::hnl::kHNLDcyNuNuNu, genie::hnl::kHNLDcyPi0Nu, genie::hnl::kHNLDcyPi0Pi0Nu, genie::hnl::kHNLDcyPiE, genie::hnl::kHNLDcyPiMu, genie::hnl::kHNLDcyPiPi0E, genie::hnl::kHNLDcyPiPi0Mu, genie::kIStInitialState, genie::kPdgHNL, genie::kPdgKP, LOG, genie::units::m, genie::units::mm, genie::utils::print::P4AsString(), pDEBUG, pFATAL, pINFO, genie::hnl::Decayer::ProcessEventRecord(), genie::hnl::VertexGenerator::ProcessEventRecord(), genie::hnl::SimpleHNL::SetEnergy(), genie::hnl::VertexGenerator::SetGeomFile(), genie::hnl::SimpleHNL::SetMomentumDirection(), genie::GHepRecord::SetPrintLevel(), and genie::InitialState::SetProbeP4().

Referenced by main().

{
  string foutName("test_decay.root");

  TFile * fout = TFile::Open( foutName.c_str(), "RECREATE" );

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

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

  if( !gOptRootGeoManager ) gOptRootGeoManager = TGeoManager::Import(gOptRootGeom.c_str()); 
  
  TGeoVolume * top_volume = gOptRootGeoManager->GetTopVolume();
  assert( top_volume );
  TGeoShape * ts  = top_volume->GetShape();
  __attribute__((unused)) TGeoBBox *  box = (TGeoBBox *)ts;
  
  LOG( "gevald_hnl", pDEBUG ) << "Imported box.";
  
  const Algorithm * algVtxGen = AlgFactory::Instance()->GetAlgorithm("genie::hnl::VertexGenerator", "Default");
  
  const VertexGenerator * vtxGen = dynamic_cast< const VertexGenerator * >( algVtxGen );
  vtxGen->SetGeomFile( gOptRootGeom );
  
  SimpleHNL sh = SimpleHNL( "HNLInstance", 0, kPdgHNL, kPdgKP,
                            gCfgMassHNL, gCfgECoupling, gCfgMCoupling, gCfgTCoupling, false );
  std::map< HNLDecayMode_t, double > valMap = sh.GetValidChannels();
  //const double CoMLifetime = sh.GetCoMLifetime();

  assert( valMap.size() > 0 ); // must be able to decay to something!
  assert( (*valMap.begin()).first == kHNLDcyNuNuNu );

  LOG( "gevald_hnl", pINFO )
    << "\n\nTesting decay modes for the HNL."
    << "\nWill process gOptNev = " << gOptNev << " x ( N_channels = " << valMap.size()
    << " ) = " << valMap.size() * gOptNev << " events, 1 for each valid channel."
    << "\nWill produce 1 ROOT file ( " << foutName << " ) with:"
    << "\n--> Energy spectrum for the decay products for each channel"
    << "\n--> Rates of each decay channel";

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

  // first set the 4-momentum of the HNL
  //gOptEnergyHNL = utils::hnl::GetCfgDouble( "HNL", "ParticleGun", "PG-Energy" );
  gOptEnergyHNL = hnlgen->GetPGunEnergy();
  double p3HNL = std::sqrt( gOptEnergyHNL * gOptEnergyHNL - gCfgMassHNL * gCfgMassHNL );
  assert( p3HNL >= 0.0 );
  TLorentzVector * p4HNL = new TLorentzVector( p3HNL * gCfgHNLCx, 
                                               p3HNL * gCfgHNLCy, 
                                               p3HNL * gCfgHNLCz, gOptEnergyHNL );

  LOG( "gevald_hnl", pDEBUG )
    << "\nUsing HNL with 4-momentum " << utils::print::P4AsString( p4HNL );
  sh.SetEnergy( gOptEnergyHNL );
  sh.SetMomentumDirection( gCfgHNLCx, gCfgHNLCy, gCfgHNLCz );

  TLorentzVector * x4HNL = new TLorentzVector( 1.0, 2.0, 3.0, 0.0 ); // dummy

  // now build array with indices of valid decay modes for speedy access
  HNLDecayMode_t validModes[10] = { kHNLDcyNull, kHNLDcyNull, kHNLDcyNull, kHNLDcyNull, kHNLDcyNull, kHNLDcyNull, kHNLDcyNull, kHNLDcyNull, kHNLDcyNull, kHNLDcyNull };
  //double validRates[10] = { -1, -1, -1, -1, -1, -1, -1, -1, -1, -1 };
  std::map< HNLDecayMode_t, double >::iterator vmit = valMap.begin(); int modeIdx = 0;
  std::ostringstream msts;
  for( ; vmit != valMap.end(); ++vmit ){
    validModes[ modeIdx ] = (*vmit).first;
    //validRates[ modeIdx ] = (*vmit).second;
    msts << "\n" << utils::hnl::AsString( (*vmit).first );
    modeIdx++;
  }

  LOG( "gevald_hnl", pDEBUG ) << "Here are the modes in order : " << msts.str();

  // declare histos
  // hSpectrum[i][j]: i iterates over HNLDecayMode_t, j over FS particle in same order as event record
  TH1D hSpectrum[10][3], hCMSpectrum[10][3], hRates;
  TH1D hParamSpace = TH1D( "hParamSpace", "Parameter space", 5, 0., 5. );

  hParamSpace.SetBinContent( 1, 1000.0 * gCfgMassHNL ); // MeV
  hParamSpace.SetBinContent( 2, gCfgECoupling );
  hParamSpace.SetBinContent( 3, gCfgMCoupling );
  hParamSpace.SetBinContent( 4, gCfgTCoupling );

  hRates = TH1D( "hRates", "Rates of HNL decay channels", 10, 0, 10 );

  std::string shortModes[10] = { "vvv", "vee", "vmue", "pi0v", "pie", "vmumu", "pimu", "pi0pi0v", 
                                 "pipi0e", "pipi0mu" };
  std::string part0names[10] = { "v1", "v", "v", "pi0", "pi", "v", "pi", "pi01", "pi", "pi" };
  std::string part1names[10] = { "v2", "e1", "mu", "v", "e", "mu1", "mu", "pi02", "pi0", "pi0" };
  std::string part2names[10] = { "v3", "e2", "e", "None", "None", "mu2", "None", "v", "e", "mu" };
  std::string partNames[3][10] = { part0names, part1names, part2names };
  for( unsigned int iChan = 0; iChan < valMap.size(); iChan++ ){

    std::string shortMode = shortModes[iChan];

    for( Int_t iPart = 0; iPart < 3; iPart++ ){
      std::string ParticleName = partNames[iPart][iChan];
      if( strcmp( ParticleName.c_str(), "None" ) != 0 ){
        hSpectrum[iChan][iPart]   = TH1D( Form( "hSpectrum_%s_%s", shortMode.c_str(), ParticleName.c_str() ),
                                          Form( "Fractional energy of particle: %s  in decay: %s", 
                                                ParticleName.c_str(), 
                                                (utils::hnl::AsString( validModes[iChan] )).c_str() ), 
                                          100, 0., 1.0 );
        hCMSpectrum[iChan][iPart] = TH1D( Form( "hCMSpectrum_%s_%s", shortMode.c_str(), ParticleName.c_str() ),
                                          Form( "Rest frame energy of particle: %s  in decay: %s", 
                                                ParticleName.c_str(), 
                                                (utils::hnl::AsString( validModes[iChan] )).c_str() ), 
                                          100, 0., gCfgMassHNL );
      } // only declare histos of particles that exist in decay
        // and are of allowed decays
    }
  }

  // fill hRates now
  double rnununu = ( (*valMap.find( kHNLDcyNuNuNu )) ).second;
  double rnuee = ( valMap.find( kHNLDcyNuEE ) != valMap.end() ) ? (*(valMap.find( kHNLDcyNuEE ))).second : 0.0;
  double rnumue = ( valMap.find( kHNLDcyNuMuE ) != valMap.end() ) ? (*(valMap.find( kHNLDcyNuMuE ))).second : 0.0;
  double rpi0nu = ( valMap.find( kHNLDcyPi0Nu ) != valMap.end() ) ? (*(valMap.find( kHNLDcyPi0Nu ))).second : 0.0;
  double rpie = ( valMap.find( kHNLDcyPiE ) != valMap.end() ) ? (*(valMap.find( kHNLDcyPiE ))).second : 0.0;
  double rnumumu = ( valMap.find( kHNLDcyNuMuMu ) != valMap.end() ) ? (*(valMap.find( kHNLDcyNuMuMu ))).second : 0.0;
  double rpimu = ( valMap.find( kHNLDcyPiMu ) != valMap.end() ) ? (*(valMap.find( kHNLDcyPiMu ))).second : 0.0;
  double rpi0pi0nu = ( valMap.find( kHNLDcyPi0Pi0Nu ) != valMap.end() ) ? (*(valMap.find( kHNLDcyPi0Pi0Nu ))).second : 0.0;
  double rpipi0e = ( valMap.find( kHNLDcyPiPi0E ) != valMap.end() ) ? (*(valMap.find( kHNLDcyPiPi0E ))).second : 0.0;
  double rpipi0mu = ( valMap.find( kHNLDcyPiPi0Mu ) != valMap.end() ) ? (*(valMap.find( kHNLDcyPiPi0Mu ))).second : 0.0;

  hRates.SetBinContent(  1, rpimu );
  hRates.SetBinContent(  2, rpie );
  hRates.SetBinContent(  3, rpi0nu );
  hRates.SetBinContent(  4, rnununu );
  hRates.SetBinContent(  5, rnumumu );
  hRates.SetBinContent(  6, rnuee );
  hRates.SetBinContent(  7, rnumue );
  hRates.SetBinContent(  8, rpipi0e );
  hRates.SetBinContent(  9, rpipi0mu );
  hRates.SetBinContent( 10, rpi0pi0nu );

  int ievent = 0;
  while( true ){
    
    if( gOptNev >= 10000 ){
      if( ievent % (gOptNev / 1000) == 0 ){
        int irat = ievent / (gOptNev / 1000);
        std::cerr << Form("%2.2f", 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 ){ std::cerr << " \n"; break; }
    
    ostringstream asts;
    for( unsigned int iMode = 0; iMode < valMap.size(); iMode++ ){
      if( ievent == 0 ){
        asts
          << "\nDecay mode " << iMode << " is " << utils::hnl::AsString( validModes[ iMode ] );
      }

      // build an event
      EventRecord * event = new EventRecord;
      Interaction * interaction = Interaction::HNL( genie::kPdgHNL, gOptEnergyHNL, validModes[ iMode ] );
      // set Particle(0) and a dummy vertex
      GHepParticle ptHNL( kPdgHNL, kIStInitialState, -1, -1, -1, -1, *p4HNL, *x4HNL );
      event->AddParticle( ptHNL );
      interaction->InitStatePtr()->SetProbeP4( *p4HNL );
      event->SetVertex( *x4HNL );

      event->SetProbability( CoMLifetime );

      event->AttachSummary( interaction );
      LOG( "gevald_hnl", pDEBUG )
        << "Simulating decay with mode " 
        << utils::hnl::AsString( (HNLDecayMode_t) interaction->ExclTag().DecayMode() );

      // simulate the decay
      hnlgen->ProcessEventRecord( event );

      // now get a weight.
      // = exp( - T_{box} / \tau_{HNL} ) = exp( - L_{box} / ( \beta_{HNL} \gamma_{HNL} c ) * h / \Gamma_{tot} )
      // placing the HNL at a point configured by the user
      std::vector< double > PGOrigin = hnlgen->GetPGunOrigin();
      double ox = PGOrigin.at(0), oy = PGOrigin.at(1), oz = PGOrigin.at(2);
      /*
      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" );
      */
      ox *= units::m / units::mm; oy *= units::m / units::mm; oz *= units::m / units::mm;
      TVector3 startPoint( ox, oy, oz ); TVector3 entryPoint, exitPoint;
      TLorentzVector tmpVtx( ox, oy, oz, 0.0 );
      event->SetVertex( tmpVtx );

      vtxGen->ProcessEventRecord(event);

      LOG( "gevald_hnl", pDEBUG ) << *event;

      // now fill the histos!
      double wgt = event->Weight();
      hSpectrum[iMode][0].Fill( (event->Particle(1))->E() / gOptEnergyHNL, wgt );
      hSpectrum[iMode][1].Fill( (event->Particle(2))->E() / gOptEnergyHNL, wgt );
      if( event->Particle(3) ) hSpectrum[iMode][2].Fill( (event->Particle(3))->E() / gOptEnergyHNL, wgt );

      // let's also fill the CM spectra
      // get particle 4-momenta and boost back to rest frame!
      TLorentzVector * p4p1 = (event->Particle(1))->GetP4();
      TLorentzVector * p4p2 = (event->Particle(2))->GetP4();
      TLorentzVector * p4p3 = 0;
      if( event->Particle(3) ) p4p3 = (event->Particle(3))->GetP4();

      TVector3 boostVec = p4HNL->BoostVector();

      p4p1->Boost( -boostVec );
      p4p2->Boost( -boostVec );
      if( p4p3 ) p4p3->Boost( -boostVec );

      hCMSpectrum[iMode][0].Fill( p4p1->E(), wgt );
      hCMSpectrum[iMode][1].Fill( p4p2->E(), wgt );
      if( p4p3 ) hCMSpectrum[iMode][2].Fill( p4p3->E(), wgt );

      // clean-up
      delete event;

    } // loop over valid decay channels
    
    if( ievent == 0 ) LOG( "gevald_hnl", pDEBUG ) << asts.str();
    LOG( "gevald_hnl", pDEBUG ) << "Finished event: " << ievent;

    ievent++;
  
  } // event loop

  fout->cd();
  hParamSpace.Write();
  hRates.Write();
  for( unsigned int i = 0; i < valMap.size(); i++ ){
    for( Int_t j = 0; j < 3; j++ ){
      std::string ParticleName = partNames[j][i];
      if( strcmp( ParticleName.c_str(), "None" ) != 0 ){
        hSpectrum[i][j].Write();
        hCMSpectrum[i][j].Write();
      }
    }
  }
  fout->Write();
  fout->Close();
  
  delete p4HNL;
  delete x4HNL;
  return 0;
}

Variable Documentation

string aunits

Definition at line 234 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

double CoMLifetime = -1

Definition at line 190 of file gBeamHNLValidationApp.cxx.

double fdx = 0

Definition at line 245 of file gBeamHNLValidationApp.cxx.

double fdy = 0

Definition at line 246 of file gBeamHNLValidationApp.cxx.

double fdz = 0

Definition at line 247 of file gBeamHNLValidationApp.cxx.

double fox = 0

Definition at line 248 of file gBeamHNLValidationApp.cxx.

double foy = 0

Definition at line 249 of file gBeamHNLValidationApp.cxx.

double foz = 0

Definition at line 250 of file gBeamHNLValidationApp.cxx.

double gCfgBoxLx = -1

Definition at line 204 of file gBeamHNLValidationApp.cxx.

double gCfgBoxLy = -1

Definition at line 205 of file gBeamHNLValidationApp.cxx.

double gCfgBoxLz = -1

Definition at line 206 of file gBeamHNLValidationApp.cxx.

HNLDecayMode_t gCfgDecayMode = kHNLDcyNull

Definition at line 194 of file gBeamHNLValidationApp.cxx.

double gCfgECoupling = -1

Definition at line 184 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig(), and TestDecay().

double gCfgHNLCx = -1

Definition at line 215 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig(), and TestDecay().

double gCfgHNLCy = -1

Definition at line 216 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig(), and TestDecay().

double gCfgHNLCz = -1

Definition at line 217 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig(), and TestDecay().

int gCfgHNLKind = 2

Definition at line 188 of file gBeamHNLValidationApp.cxx.

std::vector< HNLDecayMode_t > gCfgIntChannels

Definition at line 195 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

bool gCfgIsMajorana = false

Definition at line 187 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

double gCfgMassHNL = -1

Definition at line 183 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig(), and TestDecay().

double gCfgMCoupling = -1

Definition at line 185 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig(), and TestDecay().

double gCfgParentEnergy = -1

Definition at line 208 of file gBeamHNLValidationApp.cxx.

double gCfgParentOx = -1

Definition at line 211 of file gBeamHNLValidationApp.cxx.

double gCfgParentOy = -1

Definition at line 212 of file gBeamHNLValidationApp.cxx.

double gCfgParentOz = -1

Definition at line 213 of file gBeamHNLValidationApp.cxx.

double gCfgParentPhi = -1

Definition at line 210 of file gBeamHNLValidationApp.cxx.

double gCfgParentTheta = -1

Definition at line 209 of file gBeamHNLValidationApp.cxx.

HNLProd_t gCfgProdMode = kHNLProdNull

Definition at line 193 of file gBeamHNLValidationApp.cxx.

double gCfgTCoupling = -1

Definition at line 186 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig(), and TestDecay().

double gCfgUserAx1 = -1

Definition at line 200 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

double gCfgUserAx2 = -1

Definition at line 202 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

double gCfgUserAz = -1

Definition at line 201 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

double gCfgUserOx = -1

Definition at line 197 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

double gCfgUserOy = -1

Definition at line 198 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

double gCfgUserOz = -1

Definition at line 199 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

string geom

Definition at line 233 of file gBeamHNLValidationApp.cxx.

Referenced by GetCommandLineArgs(), GetTargetCodes(), and main().

double gOptEnergyHNL = -1

Definition at line 221 of file gBeamHNLValidationApp.cxx.

string gOptEvFilePrefix = kDefOptEvFilePrefix

Definition at line 229 of file gBeamHNLValidationApp.cxx.

double gOptGeomAUnits = 0

Definition at line 236 of file gBeamHNLValidationApp.cxx.

Referenced by ReadInConfig().

double gOptGeomLUnits = 0

Definition at line 235 of file gBeamHNLValidationApp.cxx.

double gOptGeomTUnits = 0

Definition at line 237 of file gBeamHNLValidationApp.cxx.

bool gOptIsMonoEnFlux = true

Definition at line 227 of file gBeamHNLValidationApp.cxx.

int gOptNev = 10

Definition at line 177 of file gBeamHNLValidationApp.cxx.

long int gOptRanSeed = -1

Definition at line 252 of file gBeamHNLValidationApp.cxx.

string gOptRootGeom

Definition at line 231 of file gBeamHNLValidationApp.cxx.

string gOptRootGeomTopVol = ""

Definition at line 242 of file gBeamHNLValidationApp.cxx.

Long_t gOptRunNu = 1000

Definition at line 176 of file gBeamHNLValidationApp.cxx.

bool gOptUsingRootGeom = false

Definition at line 230 of file gBeamHNLValidationApp.cxx.

HNLValidation_t gOptValidationMode = kValNone

Definition at line 179 of file gBeamHNLValidationApp.cxx.

Referenced by main().

string kDefOptEvFilePrefix = "gntp"

Definition at line 168 of file gBeamHNLValidationApp.cxx.

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

Definition at line 169 of file gBeamHNLValidationApp.cxx.

string kDefOptGeomAUnits = "rad"

Definition at line 163 of file gBeamHNLValidationApp.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 165 of file gBeamHNLValidationApp.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 162 of file gBeamHNLValidationApp.cxx.

string kDefOptGeomTUnits = "ns"

Definition at line 164 of file gBeamHNLValidationApp.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 167 of file gBeamHNLValidationApp.cxx.

string kDefOptSConfig = "BeamHNL"

Definition at line 172 of file gBeamHNLValidationApp.cxx.

string kDefOptSName = "genie::EventGenerator"

Definition at line 171 of file gBeamHNLValidationApp.cxx.

string kDefOptSTune = "GHNL20_00a_00_000"

Definition at line 173 of file gBeamHNLValidationApp.cxx.

string lunits

Definition at line 234 of file gBeamHNLValidationApp.cxx.

Referenced by GetCommandLineArgs(), and ReadInConfig().

string tunits

Definition at line 234 of file gBeamHNLValidationApp.cxx.