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