gEvGenDM.cxx File Reference

#include <cstdlib>
#include <cassert>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <TFile.h>
#include <TTree.h>
#include <TSystem.h>
#include <TVector3.h>
#include <TH1.h>
#include <TF1.h>
#include "Framework/Algorithm/AlgConfigPool.h"
#include "Framework/Conventions/XmlParserStatus.h"
#include "Framework/Conventions/GBuild.h"
#include "Framework/Conventions/Controls.h"
#include "Framework/Conventions/Constants.h"
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/EventGen/GFluxI.h"
#include "Framework/EventGen/GEVGDriver.h"
#include "Framework/EventGen/GMCJDriver.h"
#include "Framework/EventGen/GMCJMonitor.h"
#include "Framework/Interaction/Interaction.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Ntuple/NtpWriter.h"
#include "Framework/Ntuple/NtpMCFormat.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Numerical/Spline.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/PrintUtils.h"
#include "Framework/Utils/SystemUtils.h"
#include "Framework/Utils/CmdLnArgParser.h"

Go to the source code of this file.

Functions
void	GetCommandLineArgs (int argc, char **argv)
void	Initialize (void)
void	PrintSyntax (void)
bool	CheckUnitarityLimit (void)
void	GenerateEventsAtFixedInitState (void)
int	main (int argc, char **argv)

Variables
int	kDefOptNevents = 0
NtpMCFormat_t	kDefOptNtpFormat = kNFGHEP
Long_t	kDefOptRunNu = 0
int	gOptNevents
double	gOptDMEnergy
double	gOptDMEnergyRange
double	gOptDMMass
double	gOptZpCoupling
map< int, double >	gOptTgtMix
double	gOptMedRatio
Long_t	gOptRunNu
string	gOptFlux
bool	gOptWeighted
bool	gOptUsingFluxOrTgtMix = false
long int	gOptRanSeed
string	gOptInpXSecFile
string	gOptOutFileName
string	gOptStatFileName

Function Documentation

bool CheckUnitarityLimit

(

void

)

Definition at line 859 of file gEvGenDM.cxx.

{
  // Before generating the events, perform a simple sanity check
  // We estimate the leading divergent piece of the cross-section
  // We make sure it does not exceed the unitarity limit
  double gzp;
  Registry * r = AlgConfigPool::Instance()->CommonList("Param", "BoostedDarkMatter");
  r->Get("ZpCoupling", gzp);
  double gzp4 = TMath::Power(gzp,4);
  double Mzp  = gOptMedRatio * gOptDMMass;
  double Mzp2 = Mzp*Mzp;
  // The leading, forward-dominated piece is the same for both DM models
  double xsec_est = gzp4 / (4. * kPi * Mzp2);
  double ml   = gOptDMMass;
  double ml2  = ml*ml;
  double M    = kNucleonMass;
  double M2   = M*M;
  double Ed   = gOptDMEnergy;
  double Ed2  = Ed*Ed;
  double pcm2 = M2 * (Ed2 - ml2) / (ml2 + M2 + 2.*M*Ed);
  double xsec_lim = kPi / pcm2;
  bool unitary = xsec_lim > xsec_est;
  if (!unitary) {
    LOG("gevgen_dm", pWARN)
      << "Estimated a cross-section " << xsec_est/cm2 << " cm^2";
    LOG("gevgen_dm", pWARN)
      << "Unitarity limit set to " << xsec_lim/cm2 << " cm^2";
  }
  return unitary;
}

void GenerateEventsAtFixedInitState

(

void

)

Definition at line 291 of file gEvGenDM.cxx.

{
  int dark_matter = kPdgDarkMatter;
  int target   = gOptTgtMix.begin()->first;
  double Ed    = gOptDMEnergy;
  double Md    = gOptDMMass;
  double pd    = TMath::Sqrt(Ed*Ed - Md*Md);
  assert(pd>=0.);
  TLorentzVector dm_p4(0.,0.,pd,Ed); // px,py,pz,E (GeV)

  // Create init state
  InitialState init_state(target, dark_matter);

  bool unitary = CheckUnitarityLimit();
  if (!unitary) {
    LOG("gevgen_dm", pFATAL)
      << "Cross-section risks exceeding unitarity limit - Exiting";
    exit(1);
  }


  // Create/config event generation driver
  GEVGDriver evg_driver;
  evg_driver.SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
  evg_driver.SetUnphysEventMask(*RunOpt::Instance()->UnphysEventMask());
  evg_driver.Configure(init_state);

  // Initialize an Ntuple Writer
  NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu, gOptRanSeed);

  // If an output file name has been specified... use it
  if (!gOptOutFileName.empty()){
    ntpw.CustomizeFilename(gOptOutFileName);
  }
  ntpw.Initialize();


  // Create an MC Job Monitor
  GMCJMonitor mcjmonitor(gOptRunNu);
  mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());

  // If a status file name has been given... use it
  if (!gOptStatFileName.empty()){
    mcjmonitor.CustomizeFilename(gOptStatFileName);
  }


  LOG("gevgen_dm", pNOTICE)
    << "\n ** Will generate " << gOptNevents << " events for \n"
    << init_state << " at Ev = " << Ed << " GeV";

  // Generate events / print the GHEP record / add it to the ntuple
  int ievent = 0;
  while (ievent < gOptNevents) {
    LOG("gevgen_dm", pNOTICE)
      << " *** Generating event............ " << ievent;

    // generate a single event
    EventRecord * event = evg_driver.GenerateEvent(dm_p4);

    if(!event) {
      LOG("gevgen_dm", pNOTICE)
        << "Last attempt failed. Re-trying....";
      continue;
    }

    LOG("gevgen_dm", pNOTICE)
      << "Generated Event GHEP Record: " << *event;

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

  // Save the generated MC events
  ntpw.Save();
}

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 586 of file gEvGenDM.cxx.

{
  LOG("gevgen_dm", pINFO) << "Parsing command line arguments";

  // Common run options. Set defaults and read.
  RunOpt::Instance()->EnableBareXSecPreCalc(true);
  RunOpt::Instance()->ReadFromCommandLine(argc,argv);

  // Parse run options for this app

  CmdLnArgParser parser(argc,argv);

  // help?
  bool help = parser.OptionExists('h');
  if(help) {
      PrintSyntax();
      exit(0);
  }

  if ( ! parser.OptionExists("tune") ) {
    LOG("gevgen_dm", pFATAL) << "No Dark Matter tune selected, please select one ";
    LOG("gevgen_dm", pFATAL) << "Exiting ";
    exit( 0 ) ;
  }

  // number of events
  if( parser.OptionExists('n') ) {
    LOG("gevgen_dm", pINFO) << "Reading number of events to generate";
    gOptNevents = parser.ArgAsInt('n');
  } else {
    LOG("gevgen_dm", pINFO)
       << "Unspecified number of events to generate - Using default";
    gOptNevents = kDefOptNevents;
  }

  // run number
  if( parser.OptionExists('r') ) {
    LOG("gevgen_dm", pINFO) << "Reading MC run number";
    gOptRunNu = parser.ArgAsLong('r');
  } else {
    LOG("gevgen_dm", pINFO) << "Unspecified run number - Using default";
    gOptRunNu = kDefOptRunNu;
  }

  // Output file name
  if( parser.OptionExists('o') ) {
    LOG("gevgen_dm", pINFO) << "Reading output file name";
    gOptOutFileName = parser.ArgAsString('o');

    gOptStatFileName = gOptOutFileName;
    // strip the output file format and replace with .status
    if (gOptOutFileName.find_last_of(".") != string::npos)
      gOptStatFileName =
        gOptStatFileName.substr(0, gOptOutFileName.find_last_of("."));
    gOptStatFileName .append(".status");
  }

  // flux functional form
  bool using_flux = false;
  if( parser.OptionExists('f') ) {
    LOG("gevgen_dm", pINFO) << "Reading flux function";
    gOptFlux = parser.ArgAsString('f');
    using_flux = true;
  }

  if(parser.OptionExists('s')) {
    LOG("gevgen_dm", pWARN)
      << "-s option no longer available. Please read the revised code documentation";
    gAbortingInErr = true;
    exit(1);
  }


  // generate weighted events option (only relevant if using a flux)
  gOptWeighted = parser.OptionExists('w');

  // dark matter energy
  if( parser.OptionExists('e') ) {
    LOG("gevgen_dm", pINFO) << "Reading dark matter energy";
    string dme = parser.ArgAsString('e');

    // is it just a value or a range (comma separated set of values)
    if(dme.find(",") != string::npos) {
       // split the comma separated list
       vector<string> nurange = utils::str::Split(dme, ",");
       assert(nurange.size() == 2);
       double emin = atof(nurange[0].c_str());
       double emax = atof(nurange[1].c_str());
       assert(emax>emin && emin>=0);
       gOptDMEnergy      = emin;
       gOptDMEnergyRange = emax-emin;
       if(!using_flux) {
          LOG("gevgen_dm", pWARN)
             << "No flux was specified but an energy range was input!";
          LOG("gevgen_dm", pWARN)
             << "Events will be generated at fixed E = " << gOptDMEnergy << " GeV";
          gOptDMEnergyRange = -1;
       }
    } else {
       gOptDMEnergy       = atof(dme.c_str());
       gOptDMEnergyRange = -1;
    }
  } else {
    LOG("gevgen_dm", pFATAL) << "Unspecified dark matter energy - Exiting";
    PrintSyntax();
    exit(1);
  }

  // dark matter mass
  if( parser.OptionExists('m') ) {
    LOG("gevgen_dm", pINFO) << "Reading dark matter mass";
    gOptDMMass = parser.ArgAsDouble('m');
  } else {
    LOG("gevgen_dm", pFATAL) << "Unspecified dark matter mass - Exiting";
    PrintSyntax();
    exit(1);
  }

  // mediator coupling
  if( parser.OptionExists('g') ) {
    LOG("gevgen_dm", pINFO) << "Reading mediator coupling";
    gOptZpCoupling = parser.ArgAsDouble('g');
  } else {
    LOG("gevgen_dm", pINFO) << "Unspecified mediator coupling - Using value from config file";
    gOptZpCoupling = -1.;
  }

  // target mix (their PDG codes with their corresponding weights)
  bool using_tgtmix = false;
  if( parser.OptionExists('t') ) {
    LOG("gevgen_dm", pINFO) << "Reading target mix";
    string stgtmix = parser.ArgAsString('t');
    gOptTgtMix.clear();
    vector<string> tgtmix = utils::str::Split(stgtmix,",");
    if(tgtmix.size()==1) {
         int    pdg = atoi(tgtmix[0].c_str());
         double wgt = 1.0;
         gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
    } else {
      using_tgtmix = true;
      vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
      for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
         string tgt_with_wgt = *tgtmix_iter;
         string::size_type open_bracket  = tgt_with_wgt.find("[");
         string::size_type close_bracket = tgt_with_wgt.find("]");
         string::size_type ibeg = 0;
         string::size_type iend = open_bracket;
         string::size_type jbeg = open_bracket+1;
         string::size_type jend = close_bracket-1;
         int    pdg = atoi(tgt_with_wgt.substr(ibeg,iend).c_str());
         double wgt = atof(tgt_with_wgt.substr(jbeg,jend).c_str());
         LOG("Main", pNOTICE)
            << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
         gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
      }//tgtmix_iter
    }//>1

  } else {
    LOG("gevgen_dm", pFATAL) << "Unspecified target PDG code - Exiting";
    PrintSyntax();
    exit(1);
  }

  // mediator mass ratio
  if( parser.OptionExists('z') ) {
    LOG("gevgen_dm", pINFO) << "Reading mediator mass ratio";
    gOptMedRatio = parser.ArgAsDouble('z');
  } else {
    LOG("gevgen_dm", pINFO) << "Unspecified mediator mass ratio - Using default";
    gOptMedRatio = 0.5;
  }

  gOptUsingFluxOrTgtMix = using_flux || using_tgtmix;

  // random number seed
  if( parser.OptionExists("seed") ) {
    LOG("gevgen_dm", pINFO) << "Reading random number seed";
    gOptRanSeed = parser.ArgAsLong("seed");
  } else {
    LOG("gevgen_dm", pINFO) << "Unspecified random number seed - Using default";
    gOptRanSeed = -1;
  }

  // input cross-section file
  if( parser.OptionExists("cross-sections") ) {
    LOG("gevgen_dm", pINFO) << "Reading cross-section file";
    gOptInpXSecFile = parser.ArgAsString("cross-sections");
  } else {
    LOG("gevgen_dm", pINFO) << "Unspecified cross-section file";
    gOptInpXSecFile = "";
  }

  //
  // print-out the command line options
  //
  LOG("gevgen_dm", pNOTICE)
     << "\n"
     << utils::print::PrintFramedMesg("gevgen_dm job configuration");
  LOG("gevgen_dm", pNOTICE)
     << "MC Run Number: " << gOptRunNu;
  if(gOptRanSeed != -1) {
     LOG("gevgen_dm", pNOTICE)
       << "Random number seed: " << gOptRanSeed;
  } else {
     LOG("gevgen_dm", pNOTICE)
       << "Random number seed was not set, using default";
  }
  LOG("gevgen_dm", pNOTICE)
       << "Number of events requested: " << gOptNevents;
  if(gOptInpXSecFile.size() > 0) {
     LOG("gevgen_dm", pNOTICE)
       << "Using cross-section splines read from: " << gOptInpXSecFile;
  } else {
     LOG("gevgen_dm", pNOTICE)
       << "No input cross-section spline file";
  }
  LOG("gevgen_dm", pNOTICE)
       << "Flux: " << gOptFlux;
  LOG("gevgen_dm", pNOTICE)
       << "Generate weighted events? " << gOptWeighted;
  if(gOptDMEnergyRange>0) {
     LOG("gevgen_dm", pNOTICE)
        << "Dark matter energy: ["
        << gOptDMEnergy << ", " << gOptDMEnergy+gOptDMEnergyRange << "]";
  } else {
     LOG("gevgen_dm", pNOTICE)
        << "Dark matter energy: " << gOptDMEnergy;
  }
  LOG("gevgen_dm", pNOTICE)
      << "Dark matter mass: " << gOptDMMass;
  LOG("gevgen_dm", pNOTICE)
      << "Target code (PDG) & weight fraction (in case of multiple targets): ";
  LOG("gevgen_dm", pNOTICE)
      << "Mediator mass ratio: " << gOptMedRatio;
  map<int,double>::const_iterator iter;
  for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
      int    tgtpdgc = iter->first;
      double wgt     = iter->second;
      LOG("gevgen_dm", pNOTICE)
          << " >> " <<  tgtpdgc << " (weight fraction = " << wgt << ")";
  }
  LOG("gevgen_dm", pNOTICE) << "\n";

  LOG("gevgen_dm", pNOTICE) << *RunOpt::Instance();

}

void Initialize

(

void

)

Definition at line 271 of file gEvGenDM.cxx.

{

  if ( ! RunOpt::Instance()->Tune() ) {
    LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
    exit(-1);
  }
  RunOpt::Instance()->BuildTune();

  // Initialization of random number generators, cross-section table,
  // messenger thresholds, cache file
  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
  utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
  utils::app_init::RandGen(gOptRanSeed);
  utils::app_init::XSecTable(gOptInpXSecFile, false);

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

int main	(	int	argc,
		char **	argv
	)

Definition at line 235 of file gEvGenDM.cxx.

{
  GetCommandLineArgs(argc,argv);
  PDGLibrary::Instance()->AddDarkMatter(gOptDMMass,gOptMedRatio);
  if (gOptZpCoupling > 0.) {
      Registry * r = AlgConfigPool::Instance()->CommonList("Param", "BoostedDarkMatter");
      r->UnLock();
      r->Set("ZpCoupling", gOptZpCoupling);
      r->Lock();
  }
  Initialize();


  // throw on NaNs and Infs...
#if defined(HAVE_FENV_H) && defined(HAVE_FEENABLEEXCEPT)
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
#endif
  //
  // Generate dark matter events
  //

  if(gOptUsingFluxOrTgtMix) {
#ifdef __CAN_GENERATE_EVENTS_USING_A_FLUX_OR_TGTMIX__
        GenerateEventsUsingFluxOrTgtMix();
#else
  LOG("gevgen_dm", pERROR)
    << "\n   To be able to generate dark matter events from a flux and/or a target mix"
    << "\n   you need to add the following config options at your GENIE installation:"
    << "\n   --enable-flux-drivers  --enable-geom-drivers \n" ;
#endif
  } else {
     GenerateEventsAtFixedInitState();
  }
  return 0;
}

void PrintSyntax

(

void

)

Definition at line 833 of file gEvGenDM.cxx.

{
  LOG("gevgen_dm", pNOTICE)
    << "\n\n" << "Syntax:" << "\n"
    << "\n      gevgen_dm [-h]"
    << "\n                [-r run#]"
    << "\n                -n nev"
    << "\n                -e energy (or energy range) "
    << "\n                -m mass"
    << "\n                -t target_pdg "
    << "\n                [-g zp_coupling]"
    << "\n                [-z med_ratio]"
    << "\n                [-f flux_description]"
    << "\n                [-o outfile_name]"
    << "\n                [-w]"
    << "\n                [--seed random_number_seed]"
    << "\n                [--cross-sections xml_file]"
    << "\n                [--event-generator-list list_name]"
    << "\n                [--message-thresholds xml_file]"
    << "\n                [--unphysical-event-mask mask]"
    << "\n                [--event-record-print-level level]"
    << "\n                [--mc-job-status-refresh-rate  rate]"
    << "\n                [--cache-file root_file]"
    << "\n";
}

Variable Documentation

double gOptDMEnergy

Definition at line 219 of file gEvGenDM.cxx.

double gOptDMEnergyRange

Definition at line 220 of file gEvGenDM.cxx.

double gOptDMMass

Definition at line 221 of file gEvGenDM.cxx.

string gOptFlux

Definition at line 226 of file gEvGenDM.cxx.

string gOptInpXSecFile

Definition at line 230 of file gEvGenDM.cxx.

double gOptMedRatio

Definition at line 224 of file gEvGenDM.cxx.

int gOptNevents

Definition at line 218 of file gEvGenDM.cxx.

string gOptOutFileName

Definition at line 231 of file gEvGenDM.cxx.

long int gOptRanSeed

Definition at line 229 of file gEvGenDM.cxx.

Long_t gOptRunNu

Definition at line 225 of file gEvGenDM.cxx.

string gOptStatFileName

Definition at line 232 of file gEvGenDM.cxx.

map<int,double> gOptTgtMix

Definition at line 223 of file gEvGenDM.cxx.

bool gOptUsingFluxOrTgtMix = false

Definition at line 228 of file gEvGenDM.cxx.

bool gOptWeighted

Definition at line 227 of file gEvGenDM.cxx.

double gOptZpCoupling

Definition at line 222 of file gEvGenDM.cxx.

int kDefOptNevents = 0

Definition at line 213 of file gEvGenDM.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 214 of file gEvGenDM.cxx.

Long_t kDefOptRunNu = 0

Definition at line 215 of file gEvGenDM.cxx.