JetMatchingMGFastJet.cc

Go to the documentation of this file.

#include "GeneratorInterface/PartonShowerVeto/interface/JetMatchingMGFastJet.h"
#include "FWCore/Utilities/interface/Exception.h"
#include <iomanip>
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
#include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
 
#include <boost/algorithm/string/trim.hpp>
 
namespace gen {
 
  extern "C" {
 
#define PARAMLEN 20
  namespace {
    struct Param {
      Param(const std::string &str) {
        int len = std::min(PARAMLEN, (int)str.length());
        std::memcpy(value, str.c_str(), len);
        std::memset(value + len, ' ', PARAMLEN - len);
      }
 
      char value[PARAMLEN];
    };
  }  // namespace
 
  extern struct UPPRIV {
    int lnhin, lnhout;
    int mscal, ievnt;
    int ickkw, iscale;
  } uppriv_;
 
  extern struct MEMAIN {
    double etcjet, rclmax, etaclmax, qcut, showerkt, clfact;
    int maxjets, minjets, iexcfile, ktsche;
    int mektsc, nexcres, excres[30];
    int nqmatch, excproc, iexcproc[1000], iexcval[1000];
    bool nosingrad, jetprocs;
  } memain_;
 
  extern struct OUTTREE { int flag; } outtree_;
  }
 
  template <typename T>
  T JetMatchingMGFastJet::parseParameter(const std::string &value) {
    std::istringstream ss(value);
    T result;
    ss >> result;
    return result;
  }
 
  template <>
  std::string JetMatchingMGFastJet::parseParameter(const std::string &value) {
    std::string result;
    if (!result.empty() && result[0] == '\'')
      result = result.substr(1);
    if (!result.empty() && result[result.length() - 1] == '\'')
      result.resize(result.length() - 1);
    return result;
  }
 
  template <>
  bool JetMatchingMGFastJet::parseParameter(const std::string &value_) {
    std::string value(value_);
    std::transform(value.begin(), value.end(), value.begin(), (int (*)(int))std::toupper);
    return value == "T" || value == "Y" || value == "True" || value == "1" || value == ".TRUE.";
  }
 
  template <typename T>
  T JetMatchingMGFastJet::getParameter(const std::map<std::string, std::string> &params,
                                       const std::string &var,
                                       const T &defValue) {
    std::map<std::string, std::string>::const_iterator pos = params.find(var);
    if (pos == params.end())
      return defValue;
    return parseParameter<T>(pos->second);
  }
 
  template <typename T>
  T JetMatchingMGFastJet::getParameter(const std::string &var, const T &defValue) const {
    return getParameter(mgParams, var, defValue);
  }
 
  static std::map<std::string, std::string> parseHeader(const std::vector<std::string> &header) {
    std::map<std::string, std::string> params;
 
    for (std::vector<std::string>::const_iterator iter = header.begin(); iter != header.end(); ++iter) {
      std::string line = *iter;
      if (line.empty() || line[0] == '#')
        continue;
 
      std::string::size_type pos = line.find('!');
      if (pos != std::string::npos)
        line.resize(pos);
 
      pos = line.find('=');
      if (pos == std::string::npos)
        continue;
 
      std::string var = boost::algorithm::trim_copy(line.substr(pos + 1));
      std::string value = boost::algorithm::trim_copy(line.substr(0, pos));
 
      params[var] = value;
    }
 
    return params;
  }
 
  template <typename T>
  void JetMatchingMGFastJet::updateOrDie(const std::map<std::string, std::string> &params,
                                         T &param,
                                         const std::string &name) {
    if (param < 0) {
      param = getParameter(params, name, param);
    }
    if (param < 0)
      throw cms::Exception("Generator|PartonShowerVeto") << "The MGParamCMS header does not specify the jet "
                                                            "matching parameter \""
                                                         << name
                                                         << "\", but it "
                                                            "is requested by the CMSSW configuration."
                                                         << std::endl;
  }
 
  JetMatchingMGFastJet::JetMatchingMGFastJet(const edm::ParameterSet &params)
      : JetMatching(params), runInitialized(false), fJetFinder(nullptr), fIsInit(false) {
    std::string mode = params.getParameter<std::string>("mode");
    if (mode == "inclusive") {
      soup = false;
      exclusive = false;
    } else if (mode == "exclusive") {
      soup = false;
      exclusive = true;
    } else if (mode == "auto")
      soup = true;
    else
      throw cms::Exception("Generator|LHEInterface") << "MGFastJet jet matching scheme requires \"mode\" "
                                                        "parameter to be set to either \"inclusive\", "
                                                        "\"exclusive\" or \"auto\"."
                                                     << std::endl;
 
    memain_.etcjet = 0.;
    memain_.rclmax = 0.0;
    memain_.clfact = 0.0;
    memain_.ktsche = 0.0;
    memain_.etaclmax = params.getParameter<double>("MEMAIN_etaclmax");
    memain_.qcut = params.getParameter<double>("MEMAIN_qcut");
    memain_.minjets = params.getParameter<int>("MEMAIN_minjets");
    memain_.maxjets = params.getParameter<int>("MEMAIN_maxjets");
    memain_.showerkt = params.getParameter<double>("MEMAIN_showerkt");
    memain_.nqmatch = params.getParameter<int>("MEMAIN_nqmatch");
    outtree_.flag = params.getParameter<int>("outTree_flag");
    std::string list_excres = params.getParameter<std::string>("MEMAIN_excres");
    std::vector<std::string> elems;
    std::stringstream ss(list_excres);
    std::string item;
    int index = 0;
    while (std::getline(ss, item, ',')) {
      elems.push_back(item);
      memain_.excres[index] = std::atoi(item.c_str());
      index++;
    }
    memain_.nexcres = index;
 
    fDJROutFlag = params.getParameter<int>("outTree_flag");
  }
 
  double JetMatchingMGFastJet::getJetEtaMax() const { return memain_.etaclmax; }
 
  bool JetMatchingMGFastJet::initAfterBeams() {
    if (fIsInit)
      return true;
 
    //
    doMerge = uppriv_.ickkw;
    // doMerge = true;
    qCut = memain_.qcut;  //
    nQmatch = memain_.nqmatch;
    clFact = 1.;  // Default value
                  // NOTE: ME2pythia seems to default to 1.5 - need to check !!!
                  // in general, needs to read key ALPSFACT from LHE file - fix CMSSW code !!!
    nJetMin = memain_.minjets;
    nJetMax = memain_.maxjets;
 
    etaJetMax = memain_.etaclmax;
 
    coneRadius = 1.0;
    jetAlgoPower = 1;  //   this is the kT algorithm !!!
 
    // Matching procedure
    //
    qCutSq = pow(qCut, 2);
    // this should be something like memaev_.iexc
    fExcLocal = true;
 
    // If not merging, then done (?????)
    //
    // if (!doMerge) return true;
 
    // Initialise chosen jet algorithm.
    //
    fJetFinder = new fastjet::JetDefinition(fastjet::kt_algorithm, coneRadius);
    fClusJets.clear();
    fPtSortedJets.clear();
 
    fIsInit = true;
 
    return true;
  }
 
  void JetMatchingMGFastJet::beforeHadronisation(const lhef::LHEEvent *lhee) {
    if (!runInitialized)
      throw cms::Exception("Generator|PartonShowerVeto") << "Run not initialized in JetMatchingMGFastJet" << std::endl;
 
    for (int i = 0; i < 3; i++) {
      typeIdx[i].clear();
    }
 
    // Sort original process final state into light/heavy jets and 'other'.
    // Criteria:
    //   1 <= ID <= 5 and massless, or ID == 21 --> light jet (typeIdx[0])
    //   4 <= ID <= 6 and massive               --> heavy jet (typeIdx[1])
    //   All else                               --> other     (typeIdx[2])
    //
    const lhef::HEPEUP &hepeup = *lhee->getHEPEUP();
    int idx = 2;
    for (int i = 0; i < hepeup.NUP; i++) {
      if (hepeup.ISTUP[i] < 0)
        continue;
      if (hepeup.MOTHUP[i].first != 1 && hepeup.MOTHUP[i].second != 2)
        continue;  // this way we skip entries that come
                   // from resonance decays;
                   // we only take those that descent
                   // directly from "incoming partons"
      idx = 2;
      if (hepeup.IDUP[i] == ID_GLUON || (std::abs(hepeup.IDUP[i]) <= nQmatch))  // light jet
                                                                                // light jet
        idx = 0;
      else if (std::abs(hepeup.IDUP[i]) > nQmatch && std::abs(hepeup.IDUP[i]) <= ID_TOP)  // heavy jet
        idx = 1;
      // Store
      typeIdx[idx].push_back(i);
    }
 
    // NOTE: In principle, I should use exclusive, inclusive, or soup !!!
    // should be like this:
    if (soup) {
      int NPartons = typeIdx[0].size();
      fExcLocal = (NPartons < nJetMax);
    } else
      fExcLocal = exclusive;
 
    return;
  }
 
  void JetMatchingMGFastJet::init(const lhef::LHERunInfo *runInfo) {
    if (fIsInit)
      return;
 
    // read MGFastJet run card
 
    std::map<std::string, std::string> parameters;
 
    std::vector<std::string> header = runInfo->findHeader("MGRunCard");
    if (header.empty())
      throw cms::Exception("Generator|PartonShowerVeto") << "In order to use MGFastJet jet matching, "
                                                            "the input file has to contain the corresponding "
                                                            "MadGraph headers."
                                                         << std::endl;
 
    mgParams = parseHeader(header);
 
    // set variables in common block
 
    std::vector<Param> params;
    std::vector<Param> values;
    for (std::map<std::string, std::string>::const_iterator iter = mgParams.begin(); iter != mgParams.end(); ++iter) {
      params.push_back(" " + iter->first);
      values.push_back(iter->second);
    }
 
    // set MG matching parameters
 
    uppriv_.ickkw = getParameter<int>("ickkw", 0);
    memain_.mektsc = getParameter<int>("ktscheme", 0);
 
    header = runInfo->findHeader("MGParamCMS");
 
    std::map<std::string, std::string> mgInfoCMS = parseHeader(header);
 
    for (std::map<std::string, std::string>::const_iterator iter = mgInfoCMS.begin(); iter != mgInfoCMS.end(); ++iter) {
      std::cout << "mgInfoCMS: " << iter->first << " " << iter->second << std::endl;
    }
 
    updateOrDie(mgInfoCMS, memain_.etaclmax, "etaclmax");
    updateOrDie(mgInfoCMS, memain_.qcut, "qcut");
    updateOrDie(mgInfoCMS, memain_.minjets, "minjets");
    updateOrDie(mgInfoCMS, memain_.maxjets, "maxjets");
    updateOrDie(mgInfoCMS, memain_.showerkt, "showerkt");
    updateOrDie(mgInfoCMS, memain_.nqmatch, "nqmatch");
 
    runInitialized = true;
 
    initAfterBeams();
 
    fIsInit = true;
 
    return;
  }
 
  int JetMatchingMGFastJet::match(const lhef::LHEEvent *partonLevel, const std::vector<fastjet::PseudoJet> *jetInput) {
    // Number of hard partons
    //
    int NPartons = typeIdx[0].size();
 
    fClusJets.clear();
    fPtSortedJets.clear();
 
    int ClusSeqNJets = 0;
 
    fastjet::ClusterSequence ClusSequence(*jetInput, *fJetFinder);
 
    if (fExcLocal) {
      fClusJets = ClusSequence.exclusive_jets(qCutSq);
    } else {
      fClusJets = ClusSequence.inclusive_jets(qCut);
    }
 
    ClusSeqNJets = fClusJets.size();
 
    if (ClusSeqNJets < NPartons)
      return LESS_JETS;
 
    double localQcutSq = qCutSq;
 
    if (fExcLocal)  // exclusive
    {
      if (ClusSeqNJets > NPartons)
        return MORE_JETS;
    } else  // inclusive
    {
      fPtSortedJets = fastjet::sorted_by_pt(*jetInput);
      localQcutSq = std::max(qCutSq, fPtSortedJets[0].pt2());
      fClusJets = ClusSequence.exclusive_jets(NPartons);  // override
      ClusSeqNJets = NPartons;
    }
 
    if (clFact != 0)
      localQcutSq *= pow(clFact, 2);
 
    std::vector<fastjet::PseudoJet> MatchingInput;
 
    std::vector<bool> jetAssigned;
    jetAssigned.assign(fClusJets.size(), false);
 
    int counter = 0;
 
    const lhef::HEPEUP &hepeup = *partonLevel->getHEPEUP();
 
    while (counter < NPartons) {
      MatchingInput.clear();
 
      for (int i = 0; i < ClusSeqNJets; i++) {
        if (jetAssigned[i])
          continue;
        if (i == NPartons)
          break;
        //
        // this looks "awkward" but this way we do NOT pass in cluster_hist_index
        // which would actually indicate position in "history" of the "master" ClusSeq
        //
        fastjet::PseudoJet exjet = fClusJets[i];
        MatchingInput.push_back(fastjet::PseudoJet(exjet.px(), exjet.py(), exjet.pz(), exjet.e()));
        MatchingInput.back().set_user_index(i);
      }
 
      int idx = typeIdx[0][counter];
      MatchingInput.push_back(
          fastjet::PseudoJet(hepeup.PUP[idx][0], hepeup.PUP[idx][1], hepeup.PUP[idx][2], hepeup.PUP[idx][3]));
 
      //
      // in principle, one can use ClusterSequence::n_particles()
      //
      int NBefore = MatchingInput.size();
 
      // Create new clustering object - which includes the 1st clustering run !!!
      // NOTE-1: it better be a new object for each try, or history will be "too crowded".
      // NOTE-2: when created, the it ALWAYS makes at least 1 clustering step, so in most
      //         cases at least 1 new jet object will be added; although in some cases
      //         no jet is added, if the system doesn't cluster to anything (for example,
      //         input jet(s) and parton(s) are going opposite directions)
      //
      fastjet::ClusterSequence ClusSeq(MatchingInput, *fJetFinder);
 
      const std::vector<fastjet::PseudoJet> &output = ClusSeq.jets();
      int NClusJets = output.size() - NBefore;
 
      //
      // JVY - I think this is the right idea:
      // at least 1 (one) new jet needs to be added
      // however, need to double check details and refine, especially for inclusive mode
      //
      //
      if (NClusJets < 1) {
        return UNMATCHED_PARTON;
      }
      //
      // very unlikely case but let's do it just to be safe
      //
      if (NClusJets >= NBefore) {
        return MORE_JETS;
      }
 
      // Now browse history and see how close the clustering distance
      //
      // NOTE: Remember, there maybe more than one new jet in the list (for example,
      //       for process=2,3,...);
      //       in this case we take the ** first ** one that satisfies the distance/cut,
      //       which is ** typically ** also the best one
      //
      bool matched = false;
      const std::vector<fastjet::ClusterSequence::history_element> &history = ClusSeq.history();
 
      // for ( unsigned int i=nBefore; i<history.size(); i++ )
      for (unsigned int i = NBefore; i < output.size(); i++) {
        int hidx = output[i].cluster_sequence_history_index();
        double dNext = history[hidx].dij;
        if (dNext < localQcutSq) {
          //
          // the way we form input, parent1 is always jet, and parent2 can be any,
          // but if it's a good match/cluster, it should point at the parton
          //
          int parent1 = history[hidx].parent1;
          int parent2 = history[hidx].parent2;
          if (parent1 < 0 || parent2 < 0)
            break;  // bad cluster, no match
          //
          // pull up jet's "global" index
          //
          int pidx = MatchingInput[parent1].user_index();
          jetAssigned[pidx] = true;
          matched = true;
          break;
        }
      }
      if (!matched) {
        return UNMATCHED_PARTON;
      }
 
      counter++;
    }
 
    // Now save some info for DJR analysis (if requested).
    // This mimics what is done in ME2pythia.f
    // Basically, NJets and matching scale for these 4 cases:
    // 1->0, 2->1, 3->2, and 4->3
    //
    if (fDJROutFlag > 0) {
      std::vector<double> MergingScale;
      MergingScale.clear();
      for (int nj = 0; nj < 4; nj++) {
        double dmscale2 = ClusSequence.exclusive_dmerge(nj);
        double dmscale = sqrt(dmscale2);
        MergingScale.push_back(dmscale);
      }
      fDJROutput.open("events.tree", std::ios_base::app);
      double dNJets = (double)NPartons;
      fDJROutput << " " << dNJets << " " << MergingScale[0] << " " << MergingScale[1] << " " << MergingScale[2] << " "
                 << MergingScale[3] << std::endl;
      fDJROutput.close();
    }
 
    return NONE;
  }
 
}  // namespace gen