BeamHaloGenerator/src/BeamHaloProducer.cc

0001 #include <iostream>
0002 #include <ctime>
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/Run.h"
0006 #include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008
0009 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0010 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0011 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0012
0013 #include "GeneratorInterface/BeamHaloGenerator/interface/BeamHaloProducer.h"
0014 #include "GeneratorInterface/BeamHaloGenerator/interface/PYR.h"
0015
0016 using namespace edm;
0017 using namespace std;
0018
0019 #include "HepMC/IO_HEPEVT.h"
0020 #include "HepMC/HEPEVT_Wrapper.h"
0021 // #include "HepMC/ConvertHEPEVT.h"
0022 // #include "HepMC/CBhepevt.h"
0023 #include "HepMC/WeightContainer.h"
0024
0025 #define KI_BHG_INIT ki_bhg_init_
0026 extern "C" {
0027 void KI_BHG_INIT(long& seed);
0028 }
0029
0030 #define BHSETPARAM bhsetparam_
0031 extern "C" {
0032 void BHSETPARAM(int* iparam, float* fparam, const char* cparam, int length);
0033 }
0034
0035 #define KI_BHG_FILL ki_bhg_fill_
0036 extern "C" {
0037 void KI_BHG_FILL(int& iret, float& weight);
0038 }
0039
0040 #define KI_BHG_STAT ki_bhg_stat_
0041 extern "C" {
0042 void KI_BHG_STAT(int& iret);
0043 }
0044
0045 // HepMC::ConvertHEPEVT conv;
0046 //include "HepMC/HEPEVT_Wrapper.h"
0047 static HepMC::HEPEVT_Wrapper wrapper;
0048 static HepMC::IO_HEPEVT conv;
0049
0050 BeamHaloProducer::~BeamHaloProducer() {
0051   int iret = 0;
0052   call_ki_bhg_stat(iret);
0053 }
0054
0055 BeamHaloProducer::BeamHaloProducer(const ParameterSet& pset) : evt(nullptr), isInitialized_(false) {
0056   int iparam[8];
0057   float fparam[4];
0058   std::string cparam;
0059   // -- from bhgctrl.inc
0060   iparam[0] = pset.getUntrackedParameter<int>("GENMOD");
0061   iparam[1] = pset.getUntrackedParameter<int>("LHC_B1");
0062   iparam[2] = pset.getUntrackedParameter<int>("LHC_B2");
0063   iparam[3] = pset.getUntrackedParameter<int>("IW_MUO");
0064   iparam[4] = pset.getUntrackedParameter<int>("IW_HAD");
0065   iparam[5] = 9999999;
0066   iparam[6] = pset.getUntrackedParameter<int>("OFFSET", 0);
0067   iparam[7] = pset.getUntrackedParameter<int>("shift_bx");
0068
0069   fparam[0] = (float)pset.getUntrackedParameter<double>("EG_MIN");
0070   fparam[1] = (float)pset.getUntrackedParameter<double>("EG_MAX");
0071
0072   fparam[2] = (float)pset.getUntrackedParameter<double>("BXNS");
0073   fparam[3] = (float)pset.getUntrackedParameter<double>("W0", 1.0);
0074
0075   cparam = pset.getUntrackedParameter<std::string>("G3FNAME", "input.txt");
0076   call_bh_set_parameters(iparam, fparam, cparam);
0077
0078   produces<HepMCProduct>("unsmeared");
0079   produces<GenEventInfoProduct>();
0080   produces<GenRunInfoProduct, Transition::EndRun>();
0081
0082   usesResource("BeamHaloProducer");
0083
0084   cout << "BeamHaloProducer: starting event generation ... " << endl;
0085 }
0086
0087 void BeamHaloProducer::clear() {}
0088
0089 void BeamHaloProducer::setRandomEngine(CLHEP::HepRandomEngine* v) { _BeamHalo_randomEngine = v; }
0090
0091 void BeamHaloProducer::beginLuminosityBlock(LuminosityBlock const& lumi, EventSetup const&) {
0092   if (!isInitialized_) {
0093     isInitialized_ = true;
0094     RandomEngineSentry<BeamHaloProducer> randomEngineSentry(this, lumi.index());
0095
0096     // -- initialisation
0097     long seed = 1;  // This seed is not actually used
0098     call_ki_bhg_init(seed);
0099   }
0100 }
0101
0102 void BeamHaloProducer::produce(Event& e, const EventSetup& es) {
0103   RandomEngineSentry<BeamHaloProducer> randomEngineSentry(this, e.streamID());
0104
0105   // cout << "in produce " << endl;
0106
0107   //        unique_ptr<HepMCProduct> bare_product(new HepMCProduct());
0108
0109   // cout << "apres autoptr " << endl;
0110
0111   int iret = 0;
0112   float weight = 0;
0113   call_ki_bhg_fill(iret, weight);
0114
0115   // Throw an exception if call_ki_bhg_fill(...) fails.  Use the EventCorruption
0116   // exception since it maps onto SkipEvent which is what we want to do here.
0117
0118   if (iret < 0)
0119     throw edm::Exception(edm::errors::EventCorruption)
0120         << "BeamHaloProducer: function call_ki_bhg_fill returned " << iret << endl;
0121
0122   // cout << "apres fortran " << endl;
0123
0124   // HepMC::GenEvent* evt = conv.getGenEventfromHEPEVT();
0125   //    HepMC::GenEvent* evt = conv.read_next_event();  seems to be broken (?)
0126   evt = new HepMC::GenEvent();
0127
0128   for (int theindex = 1; theindex <= wrapper.number_entries(); theindex++) {
0129     HepMC::GenVertex* Vtx = new HepMC::GenVertex(
0130         HepMC::FourVector(wrapper.x(theindex), wrapper.y(theindex), wrapper.z(theindex), wrapper.t(theindex)));
0131     HepMC::FourVector p(wrapper.px(theindex), wrapper.py(theindex), wrapper.pz(theindex), wrapper.e(theindex));
0132     HepMC::GenParticle* Part = new HepMC::GenParticle(p, wrapper.id(theindex), wrapper.status(theindex));
0133     Vtx->add_particle_out(Part);
0134     evt->add_vertex(Vtx);
0135   }
0136
0137   evt->set_event_number(e.id().event());
0138
0139   HepMC::WeightContainer& weights = evt->weights();
0140   weights.push_back(weight);
0141   //    evt->print();
0142   std::unique_ptr<HepMCProduct> CMProduct(new HepMCProduct());
0143   if (evt)
0144     CMProduct->addHepMCData(evt);
0145   e.put(std::move(CMProduct), "unsmeared");
0146
0147   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(evt));
0148   e.put(std::move(genEventInfo));
0149 }
0150
0151 void BeamHaloProducer::endRunProduce(Run& run, const EventSetup& es) {
0152   // just create an empty product
0153   // to keep the EventContent definitions happy
0154   // later on we might put the info into the run info that this is a PGun
0155   unique_ptr<GenRunInfoProduct> genRunInfo(new GenRunInfoProduct());
0156   run.put(std::move(genRunInfo));
0157 }
0158
0159 bool BeamHaloProducer::call_bh_set_parameters(int* ival, float* fval, const std::string cval_string) {
0160   BHSETPARAM(ival, fval, cval_string.c_str(), cval_string.length());
0161   return true;
0162 }
0163
0164 bool BeamHaloProducer::call_ki_bhg_init(long& seed) {
0165   KI_BHG_INIT(seed);
0166   return true;
0167 }
0168
0169 bool BeamHaloProducer::call_ki_bhg_fill(int& iret, float& weight) {
0170   KI_BHG_FILL(iret, weight);
0171   return true;
0172 }
0173
0174 bool BeamHaloProducer::call_ki_bhg_stat(int& iret) {
0175   KI_BHG_STAT(iret);
0176   return true;
0177 }