ExhumeInterface/src/ExhumeHadronizer.cc

0001 #include "GeneratorInterface/ExhumeInterface/interface/ExhumeHadronizer.h"
0002
0003 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "GeneratorInterface/Core/interface/FortranCallback.h"
0007 #include "GeneratorInterface/Core/interface/FortranInstance.h"
0008 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
0009
0010 #include "HepPID/ParticleIDTranslations.hh"
0011
0012 #include "HepMC/GenEvent.h"
0013 #include "HepMC/PdfInfo.h"
0014 #include "HepMC/HEPEVT_Wrapper.h"
0015 #include "HepMC/IO_HEPEVT.h"
0016
0017 //ExHuME headers
0018 #include "GeneratorInterface/ExhumeInterface/interface/Event.h"
0019 #include "GeneratorInterface/ExhumeInterface/interface/QQ.h"
0020 #include "GeneratorInterface/ExhumeInterface/interface/GG.h"
0021 #include "GeneratorInterface/ExhumeInterface/interface/Higgs.h"
0022 #include "GeneratorInterface/ExhumeInterface/interface/DiPhoton.h"
0023
0024 #include <string>
0025 #include <sstream>
0026
0027 HepMC::IO_HEPEVT exhume_conv;
0028
0029 namespace gen {
0030   extern "C" {
0031   extern struct {
0032     int mstu[200];
0033     double paru[200];
0034     int mstj[200];
0035     double parj[200];
0036   } pydat1_;
0037 #define pydat1 pydat1_
0038
0039   extern struct {
0040     int mstp[200];
0041     double parp[200];
0042     int msti[200];
0043     double pari[200];
0044   } pypars_;
0045 #define pypars pypars_
0046
0047   extern struct {
0048     int mint[400];
0049     double vint[400];
0050   } pyint1_;
0051 #define pyint1 pyint1_
0052   }
0053
0054   extern "C" {
0055   void pylist_(int*);
0056   int pycomp_(int&);
0057   void pygive_(const char*, int);
0058   }
0059 #define pylist pylist_
0060 #define pycomp pycomp_
0061 #define pygive pygive_
0062
0063   inline void call_pylist(int mode) { pylist(&mode); }
0064   inline bool call_pygive(const std::string& line) {
0065     int numWarn = pydat1.mstu[26];  // # warnings
0066     int numErr = pydat1.mstu[22];   // # errors
0067
0068     pygive(line.c_str(), line.length());
0069
0070     return (pydat1.mstu[26] == numWarn) && (pydat1.mstu[22] == numErr);
0071   }
0072
0073   const std::vector<std::string> ExhumeHadronizer::theSharedResources = {edm::SharedResourceNames::kPythia6,
0074                                                                          gen::FortranInstance::kFortranInstance};
0075
0076   ExhumeHadronizer::ExhumeHadronizer(edm::ParameterSet const& pset)
0077       : BaseHadronizer(pset),
0078         pythia6Service_(new Pythia6Service(pset)),
0079         randomEngine_(nullptr),
0080         comEnergy_(pset.getParameter<double>("comEnergy")),
0081         myPSet_(pset),
0082         hepMCVerbosity_(pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false)),
0083         maxEventsToPrint_(pset.getUntrackedParameter<int>("maxEventsToPrint", 0)),
0084         pythiaListVerbosity_(pset.getUntrackedParameter<int>("pythiaPylistVerbosity", 0)),
0085         exhumeEvent_(nullptr) {
0086     convertToPDG_ = false;
0087     if (pset.exists("doPDGConvert")) {
0088       convertToPDG_ = pset.getParameter<bool>("doPDGConvert");
0089     }
0090
0091     //pythia6Hadronizer_ = new Pythia6Hadronizer(pset);
0092   }
0093
0094   ExhumeHadronizer::~ExhumeHadronizer() {
0095     //delete pythia6Hadronizer_;
0096     delete pythia6Service_;
0097     delete exhumeEvent_;
0098     delete exhumeProcess_;
0099   }
0100
0101   void ExhumeHadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) {
0102     pythia6Service_->setRandomEngine(v);
0103     randomEngine_ = v;
0104     if (exhumeEvent_) {
0105       exhumeEvent_->SetRandomEngine(v);
0106     }
0107   }
0108
0109   void ExhumeHadronizer::finalizeEvent() {
0110     //pythia6Hadronizer_->finalizeEvent();
0111
0112     event()->set_signal_process_id(pypars.msti[0]);
0113     event()->set_event_scale(pypars.pari[16]);
0114
0115     HepMC::PdfInfo pdf;
0116     pdf.set_id1(pyint1.mint[14] == 21 ? 0 : pyint1.mint[14]);
0117     pdf.set_id2(pyint1.mint[15] == 21 ? 0 : pyint1.mint[15]);
0118     pdf.set_x1(pyint1.vint[40]);
0119     pdf.set_x2(pyint1.vint[41]);
0120     pdf.set_pdf1(pyint1.vint[38] / pyint1.vint[40]);
0121     pdf.set_pdf2(pyint1.vint[39] / pyint1.vint[41]);
0122     pdf.set_scalePDF(pyint1.vint[50]);
0123
0124     event()->set_pdf_info(pdf);
0125
0126     event()->weights().push_back(pyint1.vint[96]);
0127
0128     // convert particle IDs Py6->PDG, if requested
0129     if (convertToPDG_) {
0130       for (HepMC::GenEvent::particle_iterator part = event()->particles_begin(); part != event()->particles_end();
0131            ++part) {
0132         (*part)->set_pdg_id(HepPID::translatePythiatoPDT((*part)->pdg_id()));
0133       }
0134     }
0135
0136     // service printouts, if requested
0137     //
0138     if (maxEventsToPrint_ > 0) {
0139       --maxEventsToPrint_;
0140       if (pythiaListVerbosity_)
0141         call_pylist(pythiaListVerbosity_);
0142       if (hepMCVerbosity_) {
0143         std::cout << "Event process = " << pypars.msti[0] << std::endl << "----------------------" << std::endl;
0144         event()->print();
0145       }
0146     }
0147
0148     return;
0149   }
0150
0151   bool ExhumeHadronizer::generatePartonsAndHadronize() {
0152     Pythia6Service::InstanceWrapper guard(pythia6Service_);
0153
0154     FortranCallback::getInstance()->resetIterationsPerEvent();
0155
0156     // generate event
0157
0158     exhumeEvent_->Generate();
0159     exhumeProcess_->Hadronise();
0160
0161     event().reset(exhume_conv.read_next_event());
0162
0163     return true;
0164   }
0165
0166   bool ExhumeHadronizer::hadronize() { return false; }
0167
0168   bool ExhumeHadronizer::decay() { return true; }
0169
0170   bool ExhumeHadronizer::residualDecay() { return true; }
0171
0172   bool ExhumeHadronizer::initializeForExternalPartons() { return false; }
0173
0174   bool ExhumeHadronizer::readSettings(int) {
0175     Pythia6Service::InstanceWrapper guard(pythia6Service_);
0176
0177     pythia6Service_->setGeneralParams();
0178
0179     return true;
0180   }
0181
0182   bool ExhumeHadronizer::initializeForInternalPartons() {
0183     Pythia6Service::InstanceWrapper guard(pythia6Service_);
0184
0185     // pythia6Service_->setGeneralParams();
0186
0187     //Exhume Initialization
0188     edm::ParameterSet processPSet = myPSet_.getParameter<edm::ParameterSet>("ExhumeProcess");
0189     std::string processType = processPSet.getParameter<std::string>("ProcessType");
0190     int sigID = -1;
0191     if (processType == "Higgs") {
0192       exhumeProcess_ = new Exhume::Higgs(myPSet_);
0193       int higgsDecay = processPSet.getParameter<int>("HiggsDecay");
0194       (static_cast<Exhume::Higgs*>(exhumeProcess_))->SetHiggsDecay(higgsDecay);
0195       sigID = 100 + higgsDecay;
0196     } else if (processType == "QQ") {
0197       exhumeProcess_ = new Exhume::QQ(myPSet_);
0198       int quarkType = processPSet.getParameter<int>("QuarkType");
0199       double thetaMin = processPSet.getParameter<double>("ThetaMin");
0200       ((Exhume::QQ*)exhumeProcess_)->SetQuarkType(quarkType);
0201       (static_cast<Exhume::QQ*>(exhumeProcess_))->SetThetaMin(thetaMin);
0202       sigID = 200 + quarkType;
0203     } else if (processType == "GG") {
0204       exhumeProcess_ = new Exhume::GG(myPSet_);
0205       double thetaMin = processPSet.getParameter<double>("ThetaMin");
0206       (static_cast<Exhume::GG*>(exhumeProcess_))->SetThetaMin(thetaMin);
0207       sigID = 300;
0208     } else if (processType == "DiPhoton") {
0209       exhumeProcess_ = new Exhume::DiPhoton(myPSet_);
0210       double thetaMin = processPSet.getParameter<double>("ThetaMin");
0211       (static_cast<Exhume::DiPhoton*>(exhumeProcess_))->SetThetaMin(thetaMin);
0212       sigID = 400;
0213     } else {
0214       throw edm::Exception(edm::errors::Configuration, "ExhumeError") << " No valid Exhume Process";
0215     }
0216
0217     pypars.msti[0] = sigID;
0218     exhumeEvent_ = new Exhume::Event(*exhumeProcess_, randomEngine_);
0219
0220     double massRangeLow = processPSet.getParameter<double>("MassRangeLow");
0221     double massRangeHigh = processPSet.getParameter<double>("MassRangeHigh");
0222     exhumeEvent_->SetMassRange(massRangeLow, massRangeHigh);
0223     exhumeEvent_->SetParameterSpace();
0224
0225     return true;
0226   }
0227
0228   bool ExhumeHadronizer::declareStableParticles(const std::vector<int>& _pdg) {
0229     std::vector<int> pdg = _pdg;
0230     //return pythia6Hadronizer_->declareStableParticles(pdg);
0231
0232     for (size_t i = 0; i < pdg.size(); i++) {
0233       int pyCode = pycomp(pdg[i]);
0234       std::ostringstream pyCard;
0235       pyCard << "MDCY(" << pyCode << ",1)=0";
0236       std::cout << pyCard.str() << std::endl;
0237       call_pygive(pyCard.str());
0238     }
0239
0240     return true;
0241   }
0242
0243   bool ExhumeHadronizer::declareSpecialSettings(const std::vector<std::string>&) { return true; }
0244
0245   void ExhumeHadronizer::statistics() {
0246     std::ostringstream footer_str;
0247
0248     double cs = exhumeEvent_->CrossSectionCalculation();
0249     double eff = exhumeEvent_->GetEfficiency();
0250     std::string name = exhumeProcess_->GetName();
0251
0252     footer_str << "\n"
0253                << "   You have just been ExHuMEd."
0254                << "\n"
0255                << "\n";
0256     footer_str << "   The cross section for process " << name << " is " << cs << " fb"
0257                << "\n"
0258                << "\n";
0259     footer_str << "   The efficiency of event generation was " << eff << "%"
0260                << "\n"
0261                << "\n";
0262
0263     edm::LogInfo("") << footer_str.str();
0264
0265     if (!runInfo().internalXSec()) {
0266       runInfo().setInternalXSec(cs);
0267     }
0268
0269     return;
0270   }
0271
0272   const char* ExhumeHadronizer::classname() const { return "gen::ExhumeHadronizer"; }
0273
0274 }  // namespace gen