File indexing completed on 2021-02-14 13:07:21
0001 #include "GeneratorInterface/PartonShowerVeto/interface/JetMatchingMGFastJet.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include <iomanip>
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
0007 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0008 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0009
0010 #include <boost/algorithm/string/trim.hpp>
0011
0012 namespace gen {
0013
0014 extern "C" {
0015
0016 #define PARAMLEN 20
0017 namespace {
0018 struct Param {
0019 Param(const std::string &str) {
0020 int len = std::min(PARAMLEN, (int)str.length());
0021 std::memcpy(value, str.c_str(), len);
0022 std::memset(value + len, ' ', PARAMLEN - len);
0023 }
0024
0025 char value[PARAMLEN];
0026 };
0027 }
0028
0029 extern struct UPPRIV {
0030 int lnhin, lnhout;
0031 int mscal, ievnt;
0032 int ickkw, iscale;
0033 } uppriv_;
0034
0035 extern struct MEMAIN {
0036 double etcjet, rclmax, etaclmax, qcut, showerkt, clfact;
0037 int maxjets, minjets, iexcfile, ktsche;
0038 int mektsc, nexcres, excres[30];
0039 int nqmatch, excproc, iexcproc[1000], iexcval[1000];
0040 bool nosingrad, jetprocs;
0041 } memain_;
0042
0043 extern struct OUTTREE { int flag; } outtree_;
0044 }
0045
0046 template <typename T>
0047 T JetMatchingMGFastJet::parseParameter(const std::string &value) {
0048 std::istringstream ss(value);
0049 T result;
0050 ss >> result;
0051 return result;
0052 }
0053
0054 template <>
0055 std::string JetMatchingMGFastJet::parseParameter(const std::string &value) {
0056 std::string result;
0057 if (!result.empty() && result[0] == '\'')
0058 result = result.substr(1);
0059 if (!result.empty() && result[result.length() - 1] == '\'')
0060 result.resize(result.length() - 1);
0061 return result;
0062 }
0063
0064 template <>
0065 bool JetMatchingMGFastJet::parseParameter(const std::string &value_) {
0066 std::string value(value_);
0067 std::transform(value.begin(), value.end(), value.begin(), (int (*)(int))std::toupper);
0068 return value == "T" || value == "Y" || value == "True" || value == "1" || value == ".TRUE.";
0069 }
0070
0071 template <typename T>
0072 T JetMatchingMGFastJet::getParameter(const std::map<std::string, std::string> ¶ms,
0073 const std::string &var,
0074 const T &defValue) {
0075 std::map<std::string, std::string>::const_iterator pos = params.find(var);
0076 if (pos == params.end())
0077 return defValue;
0078 return parseParameter<T>(pos->second);
0079 }
0080
0081 template <typename T>
0082 T JetMatchingMGFastJet::getParameter(const std::string &var, const T &defValue) const {
0083 return getParameter(mgParams, var, defValue);
0084 }
0085
0086 static std::map<std::string, std::string> parseHeader(const std::vector<std::string> &header) {
0087 std::map<std::string, std::string> params;
0088
0089 for (std::vector<std::string>::const_iterator iter = header.begin(); iter != header.end(); ++iter) {
0090 std::string line = *iter;
0091 if (line.empty() || line[0] == '#')
0092 continue;
0093
0094 std::string::size_type pos = line.find('!');
0095 if (pos != std::string::npos)
0096 line.resize(pos);
0097
0098 pos = line.find('=');
0099 if (pos == std::string::npos)
0100 continue;
0101
0102 std::string var = boost::algorithm::trim_copy(line.substr(pos + 1));
0103 std::string value = boost::algorithm::trim_copy(line.substr(0, pos));
0104
0105 params[var] = value;
0106 }
0107
0108 return params;
0109 }
0110
0111 template <typename T>
0112 void JetMatchingMGFastJet::updateOrDie(const std::map<std::string, std::string> ¶ms,
0113 T ¶m,
0114 const std::string &name) {
0115 if (param < 0) {
0116 param = getParameter(params, name, param);
0117 }
0118 if (param < 0)
0119 throw cms::Exception("Generator|PartonShowerVeto") << "The MGParamCMS header does not specify the jet "
0120 "matching parameter \""
0121 << name
0122 << "\", but it "
0123 "is requested by the CMSSW configuration."
0124 << std::endl;
0125 }
0126
0127 JetMatchingMGFastJet::JetMatchingMGFastJet(const edm::ParameterSet ¶ms)
0128 : JetMatching(params), runInitialized(false), fJetFinder(nullptr), fIsInit(false) {
0129 std::string mode = params.getParameter<std::string>("mode");
0130 if (mode == "inclusive") {
0131 soup = false;
0132 exclusive = false;
0133 } else if (mode == "exclusive") {
0134 soup = false;
0135 exclusive = true;
0136 } else if (mode == "auto")
0137 soup = true;
0138 else
0139 throw cms::Exception("Generator|LHEInterface") << "MGFastJet jet matching scheme requires \"mode\" "
0140 "parameter to be set to either \"inclusive\", "
0141 "\"exclusive\" or \"auto\"."
0142 << std::endl;
0143
0144 memain_.etcjet = 0.;
0145 memain_.rclmax = 0.0;
0146 memain_.clfact = 0.0;
0147 memain_.ktsche = 0.0;
0148 memain_.etaclmax = params.getParameter<double>("MEMAIN_etaclmax");
0149 memain_.qcut = params.getParameter<double>("MEMAIN_qcut");
0150 memain_.minjets = params.getParameter<int>("MEMAIN_minjets");
0151 memain_.maxjets = params.getParameter<int>("MEMAIN_maxjets");
0152 memain_.showerkt = params.getParameter<double>("MEMAIN_showerkt");
0153 memain_.nqmatch = params.getParameter<int>("MEMAIN_nqmatch");
0154 outtree_.flag = params.getParameter<int>("outTree_flag");
0155 std::string list_excres = params.getParameter<std::string>("MEMAIN_excres");
0156 std::vector<std::string> elems;
0157 std::stringstream ss(list_excres);
0158 std::string item;
0159 int index = 0;
0160 while (std::getline(ss, item, ',')) {
0161 elems.push_back(item);
0162 memain_.excres[index] = std::atoi(item.c_str());
0163 index++;
0164 }
0165 memain_.nexcres = index;
0166
0167 fDJROutFlag = params.getParameter<int>("outTree_flag");
0168 }
0169
0170 double JetMatchingMGFastJet::getJetEtaMax() const { return memain_.etaclmax; }
0171
0172 bool JetMatchingMGFastJet::initAfterBeams() {
0173 if (fIsInit)
0174 return true;
0175
0176
0177 doMerge = uppriv_.ickkw;
0178
0179 qCut = memain_.qcut;
0180 nQmatch = memain_.nqmatch;
0181 clFact = 1.;
0182
0183
0184 nJetMin = memain_.minjets;
0185 nJetMax = memain_.maxjets;
0186
0187 etaJetMax = memain_.etaclmax;
0188
0189 coneRadius = 1.0;
0190 jetAlgoPower = 1;
0191
0192
0193
0194 qCutSq = pow(qCut, 2);
0195
0196 fExcLocal = true;
0197
0198
0199
0200
0201
0202
0203
0204 fJetFinder = new fastjet::JetDefinition(fastjet::kt_algorithm, coneRadius);
0205 fClusJets.clear();
0206 fPtSortedJets.clear();
0207
0208 fIsInit = true;
0209
0210 return true;
0211 }
0212
0213 void JetMatchingMGFastJet::beforeHadronisation(const lhef::LHEEvent *lhee) {
0214 if (!runInitialized)
0215 throw cms::Exception("Generator|PartonShowerVeto") << "Run not initialized in JetMatchingMGFastJet" << std::endl;
0216
0217 for (int i = 0; i < 3; i++) {
0218 typeIdx[i].clear();
0219 }
0220
0221
0222
0223
0224
0225
0226
0227 const lhef::HEPEUP &hepeup = *lhee->getHEPEUP();
0228 int idx = 2;
0229 for (int i = 0; i < hepeup.NUP; i++) {
0230 if (hepeup.ISTUP[i] < 0)
0231 continue;
0232 if (hepeup.MOTHUP[i].first != 1 && hepeup.MOTHUP[i].second != 2)
0233 continue;
0234
0235
0236
0237 idx = 2;
0238 if (hepeup.IDUP[i] == ID_GLUON || (std::abs(hepeup.IDUP[i]) <= nQmatch))
0239
0240 idx = 0;
0241 else if (std::abs(hepeup.IDUP[i]) > nQmatch && std::abs(hepeup.IDUP[i]) <= ID_TOP)
0242 idx = 1;
0243
0244 typeIdx[idx].push_back(i);
0245 }
0246
0247
0248
0249 if (soup) {
0250 int NPartons = typeIdx[0].size();
0251 fExcLocal = (NPartons < nJetMax);
0252 } else
0253 fExcLocal = exclusive;
0254
0255 return;
0256 }
0257
0258 void JetMatchingMGFastJet::init(const lhef::LHERunInfo *runInfo) {
0259 if (fIsInit)
0260 return;
0261
0262
0263
0264 std::map<std::string, std::string> parameters;
0265
0266 std::vector<std::string> header = runInfo->findHeader("MGRunCard");
0267 if (header.empty())
0268 throw cms::Exception("Generator|PartonShowerVeto") << "In order to use MGFastJet jet matching, "
0269 "the input file has to contain the corresponding "
0270 "MadGraph headers."
0271 << std::endl;
0272
0273 mgParams = parseHeader(header);
0274
0275
0276
0277 std::vector<Param> params;
0278 std::vector<Param> values;
0279 for (std::map<std::string, std::string>::const_iterator iter = mgParams.begin(); iter != mgParams.end(); ++iter) {
0280 params.push_back(" " + iter->first);
0281 values.push_back(iter->second);
0282 }
0283
0284
0285
0286 uppriv_.ickkw = getParameter<int>("ickkw", 0);
0287 memain_.mektsc = getParameter<int>("ktscheme", 0);
0288
0289 header = runInfo->findHeader("MGParamCMS");
0290
0291 std::map<std::string, std::string> mgInfoCMS = parseHeader(header);
0292
0293 for (std::map<std::string, std::string>::const_iterator iter = mgInfoCMS.begin(); iter != mgInfoCMS.end(); ++iter) {
0294 std::cout << "mgInfoCMS: " << iter->first << " " << iter->second << std::endl;
0295 }
0296
0297 updateOrDie(mgInfoCMS, memain_.etaclmax, "etaclmax");
0298 updateOrDie(mgInfoCMS, memain_.qcut, "qcut");
0299 updateOrDie(mgInfoCMS, memain_.minjets, "minjets");
0300 updateOrDie(mgInfoCMS, memain_.maxjets, "maxjets");
0301 updateOrDie(mgInfoCMS, memain_.showerkt, "showerkt");
0302 updateOrDie(mgInfoCMS, memain_.nqmatch, "nqmatch");
0303
0304 runInitialized = true;
0305
0306 initAfterBeams();
0307
0308 fIsInit = true;
0309
0310 return;
0311 }
0312
0313 int JetMatchingMGFastJet::match(const lhef::LHEEvent *partonLevel, const std::vector<fastjet::PseudoJet> *jetInput) {
0314
0315
0316 int NPartons = typeIdx[0].size();
0317
0318 fClusJets.clear();
0319 fPtSortedJets.clear();
0320
0321 int ClusSeqNJets = 0;
0322
0323 fastjet::ClusterSequence ClusSequence(*jetInput, *fJetFinder);
0324
0325 if (fExcLocal) {
0326 fClusJets = ClusSequence.exclusive_jets(qCutSq);
0327 } else {
0328 fClusJets = ClusSequence.inclusive_jets(qCut);
0329 }
0330
0331 ClusSeqNJets = fClusJets.size();
0332
0333 if (ClusSeqNJets < NPartons)
0334 return LESS_JETS;
0335
0336 double localQcutSq = qCutSq;
0337
0338 if (fExcLocal)
0339 {
0340 if (ClusSeqNJets > NPartons)
0341 return MORE_JETS;
0342 } else
0343 {
0344 fPtSortedJets = fastjet::sorted_by_pt(*jetInput);
0345 localQcutSq = std::max(qCutSq, fPtSortedJets[0].pt2());
0346 fClusJets = ClusSequence.exclusive_jets(NPartons);
0347 ClusSeqNJets = NPartons;
0348 }
0349
0350 if (clFact != 0)
0351 localQcutSq *= pow(clFact, 2);
0352
0353 std::vector<fastjet::PseudoJet> MatchingInput;
0354
0355 std::vector<bool> jetAssigned;
0356 jetAssigned.assign(fClusJets.size(), false);
0357
0358 int counter = 0;
0359
0360 const lhef::HEPEUP &hepeup = *partonLevel->getHEPEUP();
0361
0362 while (counter < NPartons) {
0363 MatchingInput.clear();
0364
0365 for (int i = 0; i < ClusSeqNJets; i++) {
0366 if (jetAssigned[i])
0367 continue;
0368 if (i == NPartons)
0369 break;
0370
0371
0372
0373
0374 fastjet::PseudoJet exjet = fClusJets[i];
0375 MatchingInput.push_back(fastjet::PseudoJet(exjet.px(), exjet.py(), exjet.pz(), exjet.e()));
0376 MatchingInput.back().set_user_index(i);
0377 }
0378
0379 int idx = typeIdx[0][counter];
0380 MatchingInput.push_back(
0381 fastjet::PseudoJet(hepeup.PUP[idx][0], hepeup.PUP[idx][1], hepeup.PUP[idx][2], hepeup.PUP[idx][3]));
0382
0383
0384
0385
0386 int NBefore = MatchingInput.size();
0387
0388
0389
0390
0391
0392
0393
0394
0395 fastjet::ClusterSequence ClusSeq(MatchingInput, *fJetFinder);
0396
0397 const std::vector<fastjet::PseudoJet> &output = ClusSeq.jets();
0398 int NClusJets = output.size() - NBefore;
0399
0400
0401
0402
0403
0404
0405
0406 if (NClusJets < 1) {
0407 return UNMATCHED_PARTON;
0408 }
0409
0410
0411
0412 if (NClusJets >= NBefore) {
0413 return MORE_JETS;
0414 }
0415
0416
0417
0418
0419
0420
0421
0422
0423 bool matched = false;
0424 const std::vector<fastjet::ClusterSequence::history_element> &history = ClusSeq.history();
0425
0426
0427 for (unsigned int i = NBefore; i < output.size(); i++) {
0428 int hidx = output[i].cluster_sequence_history_index();
0429 double dNext = history[hidx].dij;
0430 if (dNext < localQcutSq) {
0431
0432
0433
0434
0435 int parent1 = history[hidx].parent1;
0436 int parent2 = history[hidx].parent2;
0437 if (parent1 < 0 || parent2 < 0)
0438 break;
0439
0440
0441
0442 int pidx = MatchingInput[parent1].user_index();
0443 jetAssigned[pidx] = true;
0444 matched = true;
0445 break;
0446 }
0447 }
0448 if (!matched) {
0449 return UNMATCHED_PARTON;
0450 }
0451
0452 counter++;
0453 }
0454
0455
0456
0457
0458
0459
0460 if (fDJROutFlag > 0) {
0461 std::vector<double> MergingScale;
0462 MergingScale.clear();
0463 for (int nj = 0; nj < 4; nj++) {
0464 double dmscale2 = ClusSequence.exclusive_dmerge(nj);
0465 double dmscale = sqrt(dmscale2);
0466 MergingScale.push_back(dmscale);
0467 }
0468 fDJROutput.open("events.tree", std::ios_base::app);
0469 double dNJets = (double)NPartons;
0470 fDJROutput << " " << dNJets << " " << MergingScale[0] << " " << MergingScale[1] << " " << MergingScale[2] << " "
0471 << MergingScale[3] << std::endl;
0472 fDJROutput.close();
0473 }
0474
0475 return NONE;
0476 }
0477
0478 }