gNNBarOscEvGen.cxx File Reference

#include <cassert>
#include <cstdlib>
#include <string>
#include <vector>
#include <sstream>
#include <TSystem.h>
#include "Framework/Algorithm/AlgFactory.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/EventGen/EventGeneratorI.h"
#include "Framework/EventGen/EventRecordVisitorI.h"
#include "Framework/EventGen/GMCJMonitor.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Ntuple/NtpWriter.h"
#include "Physics/NNBarOscillation/NNBarOscMode.h"
#include "Physics/NNBarOscillation/NNBarOscUtils.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 gNNBarOscEvGen.cxx:

Go to the source code of this file.

Functions
void	GetCommandLineArgs (int argc, char **argv)
void	PrintSyntax (void)
int	SelectAnnihilationMode (int pdg_code)
int	SelectInitState (void)
const EventRecordVisitorI *	NeutronOscGenerator (void)
int	main (int argc, char **argv)

Variables
string	kDefOptGeomLUnits = "mm"
string	kDefOptGeomDUnits = "g_cm3"
NtpMCFormat_t	kDefOptNtpFormat = kNFGHEP
string	kDefOptEvFilePrefix = "gntp"
Long_t	gOptRunNu = 1000
int	gOptNev = 10
NNBarOscMode_t	gOptDecayMode = kNONull
string	gOptEvFilePrefix = kDefOptEvFilePrefix
bool	gOptUsingRootGeom = false
map< int, double >	gOptTgtMix
string	gOptRootGeom
string	gOptRootGeomTopVol = ""
double	gOptGeomLUnits = 0
double	gOptGeomDUnits = 0
long int	gOptRanSeed = -1

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

int main	(	int	argc,
		char **	argv
	)

Definition at line 183 of file gNNBarOscEvGen.cxx.

References genie::NtpWriter::AddEventRecord(), genie::NtpWriter::CustomizeFilenamePrefix(), GetCommandLineArgs(), gOptEvFilePrefix, gOptNev, gOptRanSeed, gOptRunNu, genie::NtpWriter::Initialize(), genie::RunOpt::Instance(), kDefOptNtpFormat, LOG, genie::utils::app_init::MesgThresholds(), NeutronOscGenerator(), genie::Interaction::NOsc(), pINFO, pNOTICE, genie::EventRecordVisitorI::ProcessEventRecord(), genie::utils::app_init::RandGen(), genie::NtpWriter::Save(), SelectAnnihilationMode(), SelectInitState(), genie::GHepRecord::SetPrintLevel(), genie::GMCJMonitor::SetRefreshRate(), and genie::GMCJMonitor::Update().

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

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

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

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

  // Get the nucleon decay generator
  const EventRecordVisitorI * mcgen = NeutronOscGenerator();

  // Event loop
  int ievent = 0;
  while (1)
  {
     if(ievent == gOptNev) break;

     LOG("gevgen_nnbar_osc", pNOTICE)
          << " *** Generating event............ " << ievent;

     EventRecord * event = new EventRecord;
     int target = SelectInitState();
     int decay = SelectAnnihilationMode(target);
     Interaction * interaction = Interaction::NOsc(target,decay);
     event->AttachSummary(interaction);

     // Simulate decay     
     mcgen->ProcessEventRecord(event);

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

  return 0;
}

const EventRecordVisitorI * NeutronOscGenerator

(

void

)

Definition at line 323 of file gNNBarOscEvGen.cxx.

References genie::gAbortingInErr, genie::AlgFactory::GetAlgorithm(), genie::AlgFactory::Instance(), LOG, and pFATAL.

Referenced by main().

{
  string sname   = "genie::EventGenerator";
  string sconfig = "NNBarOsc";
  AlgFactory * algf = AlgFactory::Instance();
  const EventRecordVisitorI * mcgen =
     dynamic_cast<const EventRecordVisitorI *> (algf->GetAlgorithm(sname,sconfig));
  if(!mcgen) {
     LOG("gevgen_nnbar_osc", pFATAL)
       << "Couldn't instantiate the neutron oscillation generator";
     gAbortingInErr = true;
     exit(1);
  }
  return mcgen;
}

void PrintSyntax

(

void

)

int SelectAnnihilationMode

(

int

pdg_code

)

Definition at line 244 of file gNNBarOscEvGen.cxx.

References genie::gAbortingInErr, gOptDecayMode, genie::RandomGen::Instance(), genie::kNORandom, LOG, pERROR, pFATAL, genie::RandomGen::RndNum(), and genie::RandomGen::SetSeed().

Referenced by main().

{
  // if the mode is set to 'random' (the default), pick one at random!
  if (gOptDecayMode == kNORandom) {
    int mode;

    std::string pdg_string = std::to_string(static_cast<long long>(pdg_code));
    if (pdg_string.size() != 10) {
      LOG("gevgen_nnbar_osc", pERROR)
        << "Expecting PDG code to be a 10-digit integer; instead, it's the following: " << pdg_string;
      gAbortingInErr = true;
      exit(1);
    }

    // count number of protons & neutrons
    int n_nucleons = std::stoi(pdg_string.substr(6,3)) - 1;
    int n_protons  = std::stoi(pdg_string.substr(3,3));

    // factor proton / neutron ratio into branching ratios
    double proton_frac  = ((double)n_protons) / ((double)n_nucleons);
    double neutron_frac = 1 - proton_frac;

    // set branching ratios, taken from bubble chamber data
    const int n_modes = 16;
    double br [n_modes] = { 0.010, 0.080, 0.100, 0.220,
                            0.360, 0.160, 0.070, 0.020,
                            0.015, 0.065, 0.110, 0.280,
                            0.070, 0.240, 0.100, 0.100 };

    for (int i = 0; i < n_modes; i++) {
      if (i < 7)
        br[i] *= proton_frac;
      else
        br[i] *= neutron_frac;
    }

    // randomly generate a number between 1 and 0
    RandomGen * rnd = RandomGen::Instance();
    rnd->SetSeed(0);
    double p = rnd->RndNum().Rndm();

    // loop through all modes, figure out which one our random number corresponds to
    double threshold = 0;
    for (int i = 0; i < n_modes; i++) {
      threshold += br[i];
      if (p < threshold) {
        // once we've found our mode, return it!
        mode = i + 1;
        return mode;
      }
    }

    // error message, in case the random number selection fails
    LOG("gevgen_nnbar_osc", pFATAL) << "Random selection of final state failed!";
    gAbortingInErr = true;
    exit(1);
  }

  // if specific annihilation mode specified, just use that
  else {
    int mode = (int) gOptDecayMode;
    return mode;
  }
}

int SelectInitState

(

void

)

Definition at line 309 of file gNNBarOscEvGen.cxx.

References genie::gAbortingInErr, gOptTgtMix, LOG, and pERROR.

Referenced by main().

{
  if (gOptTgtMix.size() > 1) {
    LOG("gevgen_nnbar_osc", pERROR)
      << "Target mix not currently supported. You must specify a single target nucleus!";
    gAbortingInErr = true;
    exit(1);
  }

  int pdg_code = gOptTgtMix.begin()->first;

  return pdg_code;
}

Variable Documentation

NNBarOscMode_t gOptDecayMode = kNONull

Definition at line 172 of file gNNBarOscEvGen.cxx.

string gOptEvFilePrefix = kDefOptEvFilePrefix

Definition at line 173 of file gNNBarOscEvGen.cxx.

double gOptGeomDUnits = 0

Definition at line 179 of file gNNBarOscEvGen.cxx.

double gOptGeomLUnits = 0

Definition at line 178 of file gNNBarOscEvGen.cxx.

int gOptNev = 10

Definition at line 171 of file gNNBarOscEvGen.cxx.

long int gOptRanSeed = -1

Definition at line 180 of file gNNBarOscEvGen.cxx.

string gOptRootGeom

Definition at line 176 of file gNNBarOscEvGen.cxx.

string gOptRootGeomTopVol = ""

Definition at line 177 of file gNNBarOscEvGen.cxx.

Long_t gOptRunNu = 1000

Definition at line 170 of file gNNBarOscEvGen.cxx.

map<int,double> gOptTgtMix

Definition at line 175 of file gNNBarOscEvGen.cxx.

bool gOptUsingRootGeom = false

Definition at line 174 of file gNNBarOscEvGen.cxx.

string kDefOptEvFilePrefix = "gntp"

Definition at line 167 of file gNNBarOscEvGen.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 165 of file gNNBarOscEvGen.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 164 of file gNNBarOscEvGen.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 166 of file gNNBarOscEvGen.cxx.