simcore::generators::GenieGenerator Class Reference

Class that uses GENIE's GEVGDriver to generator eN interactions. More...

#include <GenieGenerator.h>

Public Member Functions
	GenieGenerator (const std::string &name, const framework::config::Parameters &parameters)
	Constructor.
	~GenieGenerator ()
	Destructor.
void	GeneratePrimaryVertex (G4Event *event) final override
	Generate the primary vertices in the Geant4 event.
void	RecordConfig (const std::string &id, ldmx::RunHeader &rh) final override
	Record the configuration of the primary generator into the run header.
Public Member Functions inherited from simcore::PrimaryGenerator
	PrimaryGenerator (const std::string &name, const framework::config::Parameters &parameters)
	Constructor.
	DECLARE_FACTORY_WITH_WAREHOUSE (PrimaryGenerator, std::shared_ptr< PrimaryGenerator >, const std::string &, const framework::config::Parameters &)
virtual	~PrimaryGenerator ()=default
	Destructor.
std::string	name ()
void	smearBeamspot (G4PrimaryVertex *primary_vertex)
	Apply beam spot smearing to a primary vertex.
bool	useBeamspot () const
	Check if beam spot smearing is enabled for this generator.

Private Member Functions
void	fillConfig (const framework::config::Parameters &)
bool	validateConfig ()
	fill the configuration
void	initializeGENIE ()
	simple validation check on configuration params
void	calculateTotalXS ()
	GENIE initialization.

Private Attributes
std::vector< genie::GEVGDriver >	evg_drivers_
	One GENIE event generator driver per target and convertor to HepMC3GenEvent.
genie::HepMC3Converter	hep_mc3_converter_
double	energy_
std::vector< int >	targets_
std::vector< double >	abundances_
std::vector< double >	position_
std::vector< double >	beam_size_
float	target_thickness_
float	time_
std::vector< double >	direction_
std::string	tune_
std::string	spline_file_
std::string	message_threshold_file_
std::vector< float >	ev_weighting_integral_
size_t	n_events_generated_
std::vector< size_t >	n_events_by_target_
std::vector< float >	xsec_by_target_
float	xsec_total_

Additional Inherited Members
Protected Attributes inherited from simcore::PrimaryGenerator
std::string	name_ {""}
	Name of the PrimaryGenerator.
bool	use_beamspot_ {false}
	Flag denoting whether beam spot smearing is enabled for this generator.
double	beamspot_x_size_ {0}
	Extent of the beamspot in x [mm].
double	beamspot_y_size_ {0}
	Extent of the beamspot in y [mm].
double	beamspot_z_size_ {0}
	Extent of the beamspot in z [mm].

Detailed Description

Class that uses GENIE's GEVGDriver to generator eN interactions.

Definition at line 42 of file GenieGenerator.h.

Constructor & Destructor Documentation

GenieGenerator()

simcore::generators::GenieGenerator::GenieGenerator	(	const std::string &	name,
		const framework::config::Parameters &	parameters )

Constructor.

Parameters

parameters

Parameters used to configure GENIE generator.

Parameters: energy : energy of initial electron (GeV) targets : list of ten-digit 10LZZZAAAI target codes abundances: list of relative abundances for the given targets position : position of interaction from (mm three-vector) time : time to shoot at (ns) direction : direction to shoot in (unitless three-vector) tune : name of GENIE tune seed : seed for random generator

Definition at line 236 of file GenieGenerator.cxx.

    : PrimaryGenerator(name, p) {
  fillConfig(p);
 
  if (!validateConfig())
    EXCEPTION_RAISE("ConfigurationException", "Configuration not valid.");
 
  n_events_generated_ = 0;
  initializeGENIE();
  calculateTotalXS();
}

References calculateTotalXS(), initializeGENIE(), and validateConfig().

~GenieGenerator()

simcore::generators::GenieGenerator::~GenieGenerator

(

)

Destructor.

Definition at line 249 of file GenieGenerator.cxx.

                                {
  std::cout << "--- GENIE Generation Summary BEGIN ---" << std::endl;
  double total_xsec = 0;
  for (size_t i_t = 0; i_t < targets_.size(); ++i_t) {
    std::cout << "Target=" << targets_[i_t]
              << "\tAbundance=" << abundances_[i_t]
              << "\tXSEC=" << xsec_by_target_[i_t] / genie::units::millibarn
              << " mb" << "\tEvents=" << n_events_by_target_[i_t] << std::endl;
    if (n_events_by_target_[i_t] > 0)
      total_xsec += xsec_by_target_[i_t] * abundances_[i_t];
  }
 
  std::cout << "Total events generated = " << n_events_generated_
            << "\nTotal XSEC = " << total_xsec / genie::units::millibarn
            << " mb" << std::endl;
 
  std::cout << "--- GENIE Generation Summary *END* ---" << std::endl;
}

Member Function Documentation

calculateTotalXS()

void simcore::generators::GenieGenerator::calculateTotalXS ( )

private

GENIE initialization.

Definition at line 192 of file GenieGenerator.cxx.

                                      {
  // initializing...
  xsec_total_ = 0;
  ev_weighting_integral_.resize(targets_.size(), 0.0);
  evg_drivers_.resize(targets_.size());
 
  // calculate the total xsec per target...
  for (size_t i_t = 0; i_t < targets_.size(); ++i_t) {
    genie::InitialState initial_state(targets_[i_t], 11);
    evg_drivers_[i_t].SetEventGeneratorList(
        genie::RunOpt::Instance()->EventGeneratorList());
    evg_drivers_[i_t].SetUnphysEventMask(
        *genie::RunOpt::Instance()->UnphysEventMask());
    evg_drivers_[i_t].Configure(initial_state);
    evg_drivers_[i_t].UseSplines();
 
    // setup the initial election
    TParticle initial_e;
    initial_e.SetPdgCode(11);
    auto elec_i_p = std::sqrt(energy_ * energy_ -
                              initial_e.GetMass() * initial_e.GetMass());
    initial_e.SetMomentum(elec_i_p * direction_[0], elec_i_p * direction_[1],
                          elec_i_p * direction_[2], energy_);
    TLorentzVector e_p4;
    initial_e.Momentum(e_p4);
 
    xsec_by_target_[i_t] = evg_drivers_[i_t].XSecSum(e_p4);
    xsec_total_ += xsec_by_target_[i_t] * abundances_[i_t];
 
    ev_weighting_integral_[i_t] = xsec_total_;  // running sum
 
    // print...
    ldmx_log(debug) << "Target=" << targets_[i_t]
                    << "\tAbundance=" << abundances_[i_t] << "\tXSEC="
                    << xsec_by_target_[i_t] / genie::units::millibarn << "mb";
  }
  ldmx_log(debug) << "Total XSEC = " << xsec_total_ / genie::units::millibarn
                  << " mb";
 
  // renormalize our weighting integral
  for (size_t i_t = 0; i_t < ev_weighting_integral_.size(); ++i_t)
    ev_weighting_integral_[i_t] = ev_weighting_integral_[i_t] / xsec_total_;
}

References evg_drivers_.

Referenced by GenieGenerator().

fillConfig()

void simcore::generators::GenieGenerator::fillConfig ( const framework::config::Parameters & p )

private

Definition at line 62 of file GenieGenerator.cxx.

                                                                  {
  energy_ = p.get<double>("energy");  // * GeV;
 
  targets_ = p.get<std::vector<int> >("targets");
  abundances_ = p.get<std::vector<double> >("abundances");
 
  time_ = p.get<double>("time");                          // * ns;
  position_ = p.get<std::vector<double> >("position");    // mm
  beam_size_ = p.get<std::vector<double> >("beam_size");  // mm
  direction_ = p.get<std::vector<double> >("direction");
  target_thickness_ = p.get<double>("target_thickness");  // mm
 
  tune_ = p.get<std::string>("tune");
  spline_file_ = p.get<std::string>("spline_file");
 
  message_threshold_file_ = p.get<std::string>("message_threshold_file");
}

GeneratePrimaryVertex()

void simcore::generators::GenieGenerator::GeneratePrimaryVertex ( G4Event * event )

finaloverridevirtual

Generate the primary vertices in the Geant4 event.

Parameters

event

The Geant4 event.

Implements simcore::PrimaryGenerator.

Definition at line 268 of file GenieGenerator.cxx.

                                                         {
  if (n_events_generated_ == 0) {
    // set random seed
    // have to do this here since seeds aren't properly set until we know the
    // run number
    auto seed = G4Random::getTheEngine()->getSeed();
    ldmx_log(debug) << "Initializing GENIE with seed " << seed;
    genie::utils::app_init::RandGen(seed);
  }
 
  auto nucl_target_i = 0;
 
  if (targets_.size() > 0) {
    double rand_uniform = G4Random::getTheGenerator()->flat();
 
    nucl_target_i = std::distance(
        ev_weighting_integral_.begin(),
        std::lower_bound(ev_weighting_integral_.begin(),
                         ev_weighting_integral_.end(), rand_uniform));
 
    ldmx_log(debug) << "Random number = " << rand_uniform << ", target picked "
                    << targets_.at(nucl_target_i);
  }
 
  auto x_pos = position_[0] +
               (G4Random::getTheGenerator()->flat() - 0.5) * beam_size_[0];
  auto y_pos = position_[1] +
               (G4Random::getTheGenerator()->flat() - 0.5) * beam_size_[1];
  auto z_pos = position_[2] +
               (G4Random::getTheGenerator()->flat() - 0.5) * target_thickness_;
 
  ldmx_log(debug) << "Generating interaction at (x_,y_,z_)=" << "(" << x_pos
                  << "," << y_pos << "," << z_pos << ")";
 
  // setup the initial election
  TParticle initial_e;
  initial_e.SetPdgCode(11);
  float elec_i_p =
      std::sqrt(energy_ * energy_ - initial_e.GetMass() * initial_e.GetMass());
  initial_e.SetMomentum(elec_i_p * direction_[0], elec_i_p * direction_[1],
                        elec_i_p * direction_[2], energy_);
  TLorentzVector e_p4;
  initial_e.Momentum(e_p4);
 
  ldmx_log(debug) << "Generating interation with (px,py,pz,e)=" << "("
                  << e_p4.Px() << "," << e_p4.Py() << "," << e_p4.Pz() << ","
                  << e_p4.E() << ")";
 
  n_events_by_target_[nucl_target_i] += 1;
 
  // GENIE magic — use the pre-configured driver for this target
  genie::EventRecord* genie_event = NULL;
  while (!genie_event)
    genie_event = evg_drivers_[nucl_target_i].GenerateEvent(e_p4);
 
  auto ev_info = new UserEventInformation;
  auto hepmc3_genie = hep_mc3_converter_.ConvertToHepMC3(*genie_event);
  ldmx::HepMC3GenEvent hepmc3_ldmx_genie;
  hepmc3_genie->write_data(hepmc3_ldmx_genie);
  ev_info->addHepMC3GenEvent(hepmc3_ldmx_genie);
  event->SetUserInformation(ev_info);
 
  // setup the primary vertex now
 
  G4PrimaryVertex* vertex = new G4PrimaryVertex();
  vertex->SetPosition(x_pos, y_pos, z_pos);
  vertex->SetWeight(genie_event->Weight());
 
  // loop over the entries and add to the G4Event
  int n_entries = genie_event->GetEntries();
 
  ldmx_log(debug) << "---------- " << "Generated Event "
                  << n_events_generated_ + 1 << " ----------";
 
  for (int i_p = 0; i_p < n_entries; ++i_p) {
    genie::GHepParticle* p = (genie::GHepParticle*)(*genie_event)[i_p];
 
    // make sure it's a final state particle
    if (p->Status() != 1) continue;
 
    ldmx_log(debug) << "\tAdding particle " << p->Pdg() << " with status "
                    << p->Status() << " energy " << p->E() << " ...";
 
    G4PrimaryParticle* primary = new G4PrimaryParticle();
    primary->SetPDGcode(p->Pdg());
 
    // this sets the 4momentum. But may not respect masses in G4...
    // primary->Set4Momentum(p->Px() * CLHEP::GeV, p->Py() * CLHEP::GeV,
    //                       p->Pz() * CLHEP::GeV, p->E() * CLHEP::GeV);
 
    // for now, do this to set the masses to be the G4 ones, through the PDG
    // code
    primary->SetMomentum(p->Px() * CLHEP::GeV, p->Py() * CLHEP::GeV,
                         p->Pz() * CLHEP::GeV);
 
    primary->SetProperTime(time_ * CLHEP::ns);
 
    UserPrimaryParticleInformation* primary_info =
        new UserPrimaryParticleInformation();
    primary_info->setHepEvtStatus(1);
    primary->SetUserInformation(primary_info);
 
    vertex->SetPrimary(primary);
  }
 
  // add the vertex to the event
  event->AddPrimaryVertex(vertex);
 
  // Apply beam spot smearing if configured for this generator
  if (useBeamspot()) {
    smearBeamspot(vertex);
  }
 
  ++n_events_generated_;
  delete genie_event;
}

References evg_drivers_, simcore::UserPrimaryParticleInformation::setHepEvtStatus(), simcore::PrimaryGenerator::smearBeamspot(), and simcore::PrimaryGenerator::useBeamspot().

initializeGENIE()

void simcore::generators::GenieGenerator::initializeGENIE ( )

private

simple validation check on configuration params

Definition at line 166 of file GenieGenerator.cxx.

                                     {
  // initialize some RunOpt by hacking the command line interface
  {
    char* in_arr[3] = {const_cast<char*>(""),
                       const_cast<char*>("--event-generator-list"),
                       const_cast<char*>("EM")};
    genie::RunOpt::Instance()->ReadFromCommandLine(3, in_arr);
  }
 
  // set message thresholds
  genie::utils::app_init::MesgThresholds(message_threshold_file_);
 
  // set tune info
  genie::RunOpt::Instance()->SetTuneName(tune_);
  if (!genie::RunOpt::Instance()->Tune()) {
    EXCEPTION_RAISE("ConfigurationException", "No TuneId in RunOption.");
  }
  genie::RunOpt::Instance()->BuildTune();
 
  // give it the splint file and require it
  genie::utils::app_init::XSecTable(spline_file_, true);
 
  // set GHEP print level (needed?)
  genie::GHepRecord::SetPrintLevel(0);
}

Referenced by GenieGenerator().

RecordConfig()

void simcore::generators::GenieGenerator::RecordConfig ( const std::string & id, ldmx::RunHeader & rh )

finaloverridevirtual

Record the configuration of the primary generator into the run header.

Note

you must include the id number in each entry into the run header just in case there are other generators

Implements simcore::PrimaryGenerator.

Definition at line 385 of file GenieGenerator.cxx.

                                                                        {
  rh.setStringParameter(id + " Class", "simcore::generators::GenieGenerator");
  rh.setStringParameter(id + "GenieTune", tune_);
}

References ldmx::RunHeader::setStringParameter().

validateConfig()

bool simcore::generators::GenieGenerator::validateConfig ( )

private

fill the configuration

Definition at line 80 of file GenieGenerator.cxx.

                                    {
  bool ret = true;
 
  if (targets_.size() == 0 || abundances_.size() == 0) {
    ldmx_log(error) << "targets and/or abundances sizes are zero." << "  "
                    << targets_.size() << ", " << abundances_.size();
    ret = false;
  }
  if (targets_.size() != abundances_.size()) {
    ldmx_log(error) << "targets and abundances sizes unequal." << "  "
                    << targets_.size() << " != " << abundances_.size();
    ret = false;
  }
 
  if (position_.size() != 3 || direction_.size() != 3) {
    ldmx_log(error) << "position and/or direction sizes are not 3." << "  "
                    << position_.size() << ", " << direction_.size();
    ret = false;
  }
 
  if (target_thickness_ < 0) {
    ldmx_log(warn) << "target thickness cannot be less than 0! (thickness="
                   << target_thickness_ << "). Taking absolute value.";
    target_thickness_ = std::abs(target_thickness_);
  }
 
  if (beam_size_.size() != 2) {
    if (beam_size_.size() == 0) {
      ldmx_log(info) << "beam size not set. Using zero.";
      beam_size_.resize(2);
      beam_size_[0] = 0.0;
      beam_size_[1] = 0.0;
    } else {
      ldmx_log(error) << "beam size is set, but does not have size 2." << " "
                      << beam_size_.size();
      ret = false;
    }
  } else if (beam_size_[0] < 0 || beam_size_[1] < 0) {
    ldmx_log(warn) << "Beam size set as negative value? " << "("
                   << beam_size_[0] << "," << beam_size_[1] << ")"
                   << ". Changing to positive.";
    beam_size_[0] = std::abs(beam_size_[0]);
    beam_size_[1] = std::abs(beam_size_[1]);
  }
 
  // normalize abundances
  float abundance_sum = 0;
  for (auto a : abundances_) {
    abundance_sum += a;
  }
 
  if (std::abs(abundance_sum) < 1e-6) {
    ldmx_log(error) << "abundances list sums to zero? " << abundance_sum;
    ret = false;
  }
 
  if (std::abs(abundance_sum - 1.0) > 2e-2) {
    ldmx_log(info) << "abundances list sums is not unity (" << abundance_sum
                   << " instead.) Will renormalize abundances to unity!";
  }
 
  for (size_t i_a = 0; i_a < abundances_.size(); ++i_a) {
    abundances_[i_a] = abundances_[i_a] / abundance_sum;
 
    ldmx_log(debug) << "Target=" << targets_[i_a]
                    << ", Abundance=" << abundances_[i_a];
  }
 
  float dir_total_sq = 0;
  for (auto d : direction_) dir_total_sq += d * d;
 
  if (dir_total_sq < 1e-6) {
    ldmx_log(error) << "direction vector is zero or negative? " << "("
                    << direction_[0] << "," << direction_[1] << ","
                    << direction_[2] << ")";
    ret = false;
  }
  for (size_t i_d = 0; i_d < direction_.size(); ++i_d)
    direction_[i_d] = direction_[i_d] / std::sqrt(dir_total_sq);
 
  xsec_by_target_.resize(targets_.size(), -999.);
  n_events_by_target_.resize(targets_.size(), 0);
 
  return ret;
}

Referenced by GenieGenerator().

Member Data Documentation

abundances_

std::vector<double> simcore::generators::GenieGenerator::abundances_

private

Definition at line 85 of file GenieGenerator.h.

beam_size_

std::vector<double> simcore::generators::GenieGenerator::beam_size_

private

Definition at line 87 of file GenieGenerator.h.

direction_

std::vector<double> simcore::generators::GenieGenerator::direction_

private

Definition at line 90 of file GenieGenerator.h.

energy_

double simcore::generators::GenieGenerator::energy_

private

Definition at line 83 of file GenieGenerator.h.

ev_weighting_integral_

std::vector<float> simcore::generators::GenieGenerator::ev_weighting_integral_

private

Definition at line 97 of file GenieGenerator.h.

evg_drivers_

std::vector<genie::GEVGDriver> simcore::generators::GenieGenerator::evg_drivers_

private

One GENIE event generator driver per target and convertor to HepMC3GenEvent.

Each driver is configured once for its target during initialization.

Definition at line 80 of file GenieGenerator.h.

Referenced by calculateTotalXS(), and GeneratePrimaryVertex().

hep_mc3_converter_

genie::HepMC3Converter simcore::generators::GenieGenerator::hep_mc3_converter_

private

Definition at line 81 of file GenieGenerator.h.

message_threshold_file_

std::string simcore::generators::GenieGenerator::message_threshold_file_

private

Definition at line 95 of file GenieGenerator.h.

n_events_by_target_

std::vector<size_t> simcore::generators::GenieGenerator::n_events_by_target_

private

Definition at line 99 of file GenieGenerator.h.

n_events_generated_

size_t simcore::generators::GenieGenerator::n_events_generated_

private

Definition at line 98 of file GenieGenerator.h.

position_

std::vector<double> simcore::generators::GenieGenerator::position_

private

Definition at line 86 of file GenieGenerator.h.

spline_file_

std::string simcore::generators::GenieGenerator::spline_file_

private

Definition at line 93 of file GenieGenerator.h.

target_thickness_

float simcore::generators::GenieGenerator::target_thickness_

private

Definition at line 88 of file GenieGenerator.h.

targets_

std::vector<int> simcore::generators::GenieGenerator::targets_

private

Definition at line 84 of file GenieGenerator.h.

time_

float simcore::generators::GenieGenerator::time_

private

Definition at line 89 of file GenieGenerator.h.

tune_

std::string simcore::generators::GenieGenerator::tune_

private

Definition at line 92 of file GenieGenerator.h.

xsec_by_target_

std::vector<float> simcore::generators::GenieGenerator::xsec_by_target_

private

Definition at line 100 of file GenieGenerator.h.

xsec_total_

float simcore::generators::GenieGenerator::xsec_total_

private

Definition at line 102 of file GenieGenerator.h.

---

The documentation for this class was generated from the following files:

• SimCore/include/SimCore/Generators/GenieGenerator.h
• SimCore/src/SimCore/Generators/GenieGenerator.cxx