BeamHaloProducer.cc

Go to the documentation of this file.

#include <iostream>
#include <ctime>

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Run.h"
#include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"
#include "FWCore/Utilities/interface/Exception.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"

#include "GeneratorInterface/BeamHaloGenerator/interface/BeamHaloProducer.h"
#include "GeneratorInterface/BeamHaloGenerator/interface/PYR.h"

using namespace edm;
using namespace std;

#include "HepMC/IO_HEPEVT.h"
#include "HepMC/HEPEVT_Wrapper.h"
// #include "HepMC/ConvertHEPEVT.h"
// #include "HepMC/CBhepevt.h"
#include "HepMC/WeightContainer.h"

#define KI_BHG_INIT ki_bhg_init_
extern "C" {
void KI_BHG_INIT(long& seed);
}

#define BHSETPARAM bhsetparam_
extern "C" {
void BHSETPARAM(int* iparam, float* fparam, const char* cparam, int length);
}

#define KI_BHG_FILL ki_bhg_fill_
extern "C" {
void KI_BHG_FILL(int& iret, float& weight);
}

#define KI_BHG_STAT ki_bhg_stat_
extern "C" {
void KI_BHG_STAT(int& iret);
}

// HepMC::ConvertHEPEVT conv;
//include "HepMC/HEPEVT_Wrapper.h"
static HepMC::HEPEVT_Wrapper wrapper;
static HepMC::IO_HEPEVT conv;

BeamHaloProducer::~BeamHaloProducer() {
  int iret = 0;
  call_ki_bhg_stat(iret);
}

BeamHaloProducer::BeamHaloProducer(const ParameterSet& pset) : evt(nullptr), isInitialized_(false) {
  int iparam[8];
  float fparam[4];
  std::string cparam;
  // -- from bhgctrl.inc
  iparam[0] = pset.getUntrackedParameter<int>("GENMOD");
  iparam[1] = pset.getUntrackedParameter<int>("LHC_B1");
  iparam[2] = pset.getUntrackedParameter<int>("LHC_B2");
  iparam[3] = pset.getUntrackedParameter<int>("IW_MUO");
  iparam[4] = pset.getUntrackedParameter<int>("IW_HAD");
  iparam[5] = 9999999;
  iparam[6] = pset.getUntrackedParameter<int>("OFFSET", 0);
  iparam[7] = pset.getUntrackedParameter<int>("shift_bx");

  fparam[0] = (float)pset.getUntrackedParameter<double>("EG_MIN");
  fparam[1] = (float)pset.getUntrackedParameter<double>("EG_MAX");

  fparam[2] = (float)pset.getUntrackedParameter<double>("BXNS");
  fparam[3] = (float)pset.getUntrackedParameter<double>("W0", 1.0);

  cparam = pset.getUntrackedParameter<std::string>("G3FNAME", "input.txt");
  call_bh_set_parameters(iparam, fparam, cparam);

  produces<HepMCProduct>("unsmeared");
  produces<GenEventInfoProduct>();
  produces<GenRunInfoProduct, Transition::EndRun>();

  usesResource("BeamHaloProducer");

  cout << "BeamHaloProducer: starting event generation ... " << endl;
}

void BeamHaloProducer::clear() {}

void BeamHaloProducer::setRandomEngine(CLHEP::HepRandomEngine* v) { _BeamHalo_randomEngine = v; }

void BeamHaloProducer::beginLuminosityBlock(LuminosityBlock const& lumi, EventSetup const&) {
  if (!isInitialized_) {
    isInitialized_ = true;
    RandomEngineSentry<BeamHaloProducer> randomEngineSentry(this, lumi.index());

    // -- initialisation
    long seed = 1;  // This seed is not actually used
    call_ki_bhg_init(seed);
  }
}

void BeamHaloProducer::produce(Event& e, const EventSetup& es) {
  RandomEngineSentry<BeamHaloProducer> randomEngineSentry(this, e.streamID());

  // cout << "in produce " << endl;

  //        unique_ptr<HepMCProduct> bare_product(new HepMCProduct());

  // cout << "apres autoptr " << endl;

  int iret = 0;
  float weight = 0;
  call_ki_bhg_fill(iret, weight);

  // Throw an exception if call_ki_bhg_fill(...) fails.  Use the EventCorruption
  // exception since it maps onto SkipEvent which is what we want to do here.

  if (iret < 0)
    throw edm::Exception(edm::errors::EventCorruption)
        << "BeamHaloProducer: function call_ki_bhg_fill returned " << iret << endl;

  // cout << "apres fortran " << endl;

  // HepMC::GenEvent* evt = conv.getGenEventfromHEPEVT();
  //    HepMC::GenEvent* evt = conv.read_next_event();  seems to be broken (?)
  evt = new HepMC::GenEvent();

  for (int theindex = 1; theindex <= wrapper.number_entries(); theindex++) {
    HepMC::GenVertex* Vtx = new HepMC::GenVertex(
        HepMC::FourVector(wrapper.x(theindex), wrapper.y(theindex), wrapper.z(theindex), wrapper.t(theindex)));
    HepMC::FourVector p(wrapper.px(theindex), wrapper.py(theindex), wrapper.pz(theindex), wrapper.e(theindex));
    HepMC::GenParticle* Part = new HepMC::GenParticle(p, wrapper.id(theindex), wrapper.status(theindex));
    Vtx->add_particle_out(Part);
    evt->add_vertex(Vtx);
  }

  evt->set_event_number(e.id().event());

  HepMC::WeightContainer& weights = evt->weights();
  weights.push_back(weight);
  //    evt->print();
  std::unique_ptr<HepMCProduct> CMProduct(new HepMCProduct());
  if (evt)
    CMProduct->addHepMCData(evt);
  e.put(std::move(CMProduct), "unsmeared");

  unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(evt));
  e.put(std::move(genEventInfo));
}

void BeamHaloProducer::endRunProduce(Run& run, const EventSetup& es) {
  // just create an empty product
  // to keep the EventContent definitions happy
  // later on we might put the info into the run info that this is a PGun
  unique_ptr<GenRunInfoProduct> genRunInfo(new GenRunInfoProduct());
  run.put(std::move(genRunInfo));
}

bool BeamHaloProducer::call_bh_set_parameters(int* ival, float* fval, const std::string cval_string) {
  BHSETPARAM(ival, fval, cval_string.c_str(), cval_string.length());
  return true;
}

bool BeamHaloProducer::call_ki_bhg_init(long& seed) {
  KI_BHG_INIT(seed);
  return true;
}

bool BeamHaloProducer::call_ki_bhg_fill(int& iret, float& weight) {
  KI_BHG_FILL(iret, weight);
  return true;
}

bool BeamHaloProducer::call_ki_bhg_stat(int& iret) {
  KI_BHG_STAT(iret);
  return true;
}