PartonShowerVeto/src/JetMatchingAlpgen.cc

0001 #include <iostream>
0002 #include <cstring>
0003 #include <vector>
0004 #include <memory>
0005 #include <string>
0006
0007 #include <HepMC/GenEvent.h>
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011
0012 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0013 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0014 #include "GeneratorInterface/PartonShowerVeto/interface/JetMatchingAlpgen.h"
0015
0016 namespace gen {
0017
0018   extern "C" {
0019   // Put here all the functions for interfacing Fortran code.
0020
0021   //    extern void alinit_();
0022   extern void alsetp_();
0023   //    extern void alevnt_();
0024   extern void alveto_(int* ipveto);
0025
0026   extern void dbpart_();
0027   extern void pyupre_();
0028
0029   void alshcd_(char csho[3]);
0030   void alshen_();
0031   }  // extern "C"
0032
0033   // Constructor
0034   JetMatchingAlpgen::JetMatchingAlpgen(const edm::ParameterSet& params)
0035       : JetMatching(params), applyMatching(params.getParameter<bool>("applyMatching")), runInitialized(false) {
0036     ahopts_.etclus = params.getParameter<double>("etMin");
0037     ahopts_.rclus = params.getParameter<double>("drMin");
0038     ahopts_.iexc = params.getParameter<bool>("exclusive");
0039   }
0040
0041   // Destructor
0042   JetMatchingAlpgen::~JetMatchingAlpgen() {}
0043
0044   std::set<std::string> JetMatchingAlpgen::capabilities() const {
0045     std::set<std::string> result;
0046     result.insert("psFinalState");
0047     result.insert("hepevt");
0048     result.insert("pythia6");
0049     return result;
0050   }
0051
0052   // Implements the Alpgen MLM method - use alpg_match.F
0053
0054   void JetMatchingAlpgen::init(const lhef::LHERunInfo* runInfo) {
0055     // Read Alpgen run card stored in the LHERunInfo object.
0056     std::vector<std::string> headerLines = runInfo->findHeader("AlpgenUnwParFile");
0057     if (headerLines.empty())
0058       throw cms::Exception("Generator|PartonShowerVeto") << "In order to use Alpgen jet matching, "
0059                                                             "the input file has to contain the corresponding "
0060                                                             "Alpgen headers."
0061                                                          << std::endl;
0062
0063     // Parse the header using its bultin function.
0064     header.parse(headerLines.begin(), headerLines.end());
0065
0066     // I don't want to print this right now.
0067     //  std::cout << "Alpgen header" << std::endl;
0068     //  std::cout << "========================" << std::endl;
0069     //  std::cout << "\tihrd = " << header.ihrd << std::endl;
0070     //  std::cout << "\tmc = " << header.masses[AlpgenHeader::mc]
0071     //            << ", mb = " << header.masses[AlpgenHeader::mb]
0072     //            << ", mt = " << header.masses[AlpgenHeader::mt]
0073     //            << ", mw = " << header.masses[AlpgenHeader::mw]
0074     //            << ", mz = " << header.masses[AlpgenHeader::mz]
0075     //            << ", mh = " << header.masses[AlpgenHeader::mh]
0076     //            << std::endl;
0077     //  for(std::map<AlpgenHeader::Parameter, double>::const_iterator iter =
0078     //      header.params.begin(); iter != header.params.end(); ++iter)
0079     //      std::cout << "\t" << AlpgenHeader::parameterName(iter->first)
0080     //                << " = " << iter->second << std::endl;
0081     //  std::cout << "\txsec = " << header.xsec
0082     //            << " +-" << header.xsecErr << std::endl;
0083     //  std::cout << "========================" << std::endl;
0084
0085     // Here we pass a few header variables to common block and
0086     // call Alpgen init routine to do the rest.
0087     // The variables passed first are the ones directly that
0088     // need to be set up "manually": IHRD and the masses.
0089     // (ebeam is set just to mimic the original code)
0090     // Afterwards, we pass the full spectrum of Alpgen
0091     // parameters directly into the AHPARS structure, to be
0092     // treated by AHSPAR which is called inside alsetp_().
0093
0094     std::copy(header.masses, header.masses + AlpgenHeader::MASS_MAX, ahppara_.masses);
0095     ahppara_.ihrd = header.ihrd;
0096     ahppara_.ebeam = header.ebeam;
0097
0098     for (std::map<AlpgenHeader::Parameter, double>::const_iterator iter = header.params.begin();
0099          iter != header.params.end();
0100          ++iter) {
0101       if (iter->first <= 0 || iter->first >= (int)AHPARS::nparam - 1)
0102         continue;
0103       ahpars_.parval[(int)iter->first - 1] = iter->second;
0104     }
0105
0106     // Run the rest of the setup.
0107     alsetp_();
0108
0109     // When we reach this point, the run is fully initialized.
0110     runInitialized = true;
0111   }
0112
0113   void JetMatchingAlpgen::beforeHadronisation(const lhef::LHEEvent* event) {
0114     // We can't continue if the run has not been initialized.
0115     if (!runInitialized)
0116       throw cms::Exception("Generator|PartonShowerVeto") << "Run not initialized in JetMatchingAlpgen" << std::endl;
0117
0118     // We are called just after LHEInterface has filled in
0119     // the Fortran common block (and Pythia6 called UPEVNT).
0120
0121     // Possibly not interesting for us.
0122     // (except perhaps for debugging?)
0123     //  pyupre_();
0124     //  dbpart_();
0125     eventInitialized = true;
0126   }
0127
0128   /*
0129 int JetMatchingAlpgen::match(const HepMC::GenEvent *partonLevel,
0130                  const HepMC::GenEvent *finalState,
0131                  bool showeredFinalState)
0132 */
0133   int JetMatchingAlpgen::match(const lhef::LHEEvent* partonLevel, const std::vector<fastjet::PseudoJet>* jetInput) {
0134     /*
0135   if (!showeredFinalState)
0136     throw cms::Exception("Generator|PartonShowerVeto")
0137       << "Alpgen matching expected parton shower "
0138       "final state." << std::endl;
0139 */
0140
0141     if (!runInitialized)
0142       throw cms::Exception("Generator|PartonShowerVeto") << "Run not initialized in JetMatchingAlpgen" << std::endl;
0143
0144     if (!eventInitialized)
0145       throw cms::Exception("Generator|PartonShowerVeto") << "Event not initialized in JetMatchingAlpgen" << std::endl;
0146
0147     // If matching not required (e.g., icckw = 0), don't run the
0148     // FORTRAN veto code.
0149     if (!applyMatching)
0150       return 0;
0151
0152     // Call the Fortran veto code.
0153     int veto = 0;
0154     alveto_(&veto);
0155
0156     eventInitialized = false;
0157
0158     // If event was vetoed, the variable veto will contain the number 1.
0159     // In this case, we must return 1 - that will be used as the return value from UPVETO.
0160     // If event was accepted, the variable veto will contain the number 0.
0161     // In this case, we must return 0 - that will be used as the return value from UPVETO.
0162     return veto ? 1 : 0;
0163   }
0164
0165   void alshcd_(char csho[3]) {
0166     std::memcpy(csho, "PYT", 3);  // or "HER"
0167   }
0168
0169   void alshen_() {}
0170
0171 }  // end namespace gen