File indexing completed on 2023-03-17 11:04:49
0001 #include <algorithm>
0002 #include <iostream>
0003 #include <iomanip>
0004 #include <string>
0005 #include <cctype>
0006 #include <vector>
0007 #include <memory>
0008 #include <cmath>
0009 #include <cstring>
0010
0011 #include <xercesc/dom/DOM.hpp>
0012 #include <xercesc/parsers/XercesDOMParser.hpp>
0013 #include <xercesc/sax/HandlerBase.hpp>
0014
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 #include "SimDataFormats/GeneratorProducts/interface/LesHouches.h"
0019
0020 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0021
0022 #include "XMLUtils.h"
0023
0024 XERCES_CPP_NAMESPACE_USE
0025
0026 static int skipWhitespace(std::istream &in) {
0027 int ch;
0028 do {
0029 ch = in.get();
0030 } while (std::isspace(ch));
0031 if (ch != std::istream::traits_type::eof())
0032 in.putback(ch);
0033 return ch;
0034 }
0035
0036 namespace lhef {
0037
0038 LHERunInfo::LHERunInfo(std::istream &in) {
0039 in >> heprup.IDBMUP.first >> heprup.IDBMUP.second >> heprup.EBMUP.first >> heprup.EBMUP.second >>
0040 heprup.PDFGUP.first >> heprup.PDFGUP.second >> heprup.PDFSUP.first >> heprup.PDFSUP.second >> heprup.IDWTUP >>
0041 heprup.NPRUP;
0042 if (!in.good())
0043 throw cms::Exception("InvalidFormat") << "Les Houches file contained invalid"
0044 " header in init section."
0045 << std::endl;
0046
0047 heprup.resize();
0048
0049 for (int i = 0; i < heprup.NPRUP; i++) {
0050 in >> heprup.XSECUP[i] >> heprup.XERRUP[i] >> heprup.XMAXUP[i] >> heprup.LPRUP[i];
0051 if (!in.good())
0052 throw cms::Exception("InvalidFormat") << "Les Houches file contained invalid data"
0053 " in header payload line "
0054 << (i + 1) << "." << std::endl;
0055 }
0056
0057 while (skipWhitespace(in) == '#') {
0058 std::string line;
0059 std::getline(in, line);
0060 comments.push_back(line + "\n");
0061 }
0062
0063 if (!in.eof())
0064 edm::LogInfo("Generator|LHEInterface")
0065 << "Les Houches file contained spurious"
0066 " content after the regular data (this is normal for LHEv3 files for now)."
0067 << std::endl;
0068
0069 init();
0070 }
0071
0072 LHERunInfo::LHERunInfo(const HEPRUP &heprup) : heprup(heprup) { init(); }
0073
0074 LHERunInfo::LHERunInfo(const HEPRUP &heprup,
0075 const std::vector<LHERunInfoProduct::Header> &headers,
0076 const std::vector<std::string> &comments)
0077 : heprup(heprup) {
0078 std::copy(headers.begin(), headers.end(), std::back_inserter(this->headers));
0079 std::copy(comments.begin(), comments.end(), std::back_inserter(this->comments));
0080
0081 init();
0082 }
0083
0084 LHERunInfo::LHERunInfo(const LHERunInfoProduct &product) : heprup(product.heprup()) {
0085 std::copy(product.headers_begin(), product.headers_end(), std::back_inserter(headers));
0086 std::copy(product.comments_begin(), product.comments_end(), std::back_inserter(comments));
0087
0088 init();
0089 }
0090
0091 LHERunInfo::~LHERunInfo() {}
0092
0093 void LHERunInfo::init() {
0094 for (int i = 0; i < heprup.NPRUP; i++) {
0095 Process proc;
0096 proc.setProcess(heprup.LPRUP[i]);
0097 proc.setHepRupIndex((unsigned int)i);
0098 processes.push_back(proc);
0099 }
0100
0101 std::sort(processes.begin(), processes.end());
0102 }
0103
0104 void LHERunInfo::initLumi() {
0105 processesLumi.clear();
0106 for (int i = 0; i < heprup.NPRUP; i++) {
0107 Process proc;
0108 proc.setProcess(heprup.LPRUP[i]);
0109 proc.setHepRupIndex((unsigned int)i);
0110 proc.setLHEXSec(heprup.XSECUP[i], heprup.XERRUP[i]);
0111 processesLumi.push_back(proc);
0112 }
0113
0114 std::sort(processesLumi.begin(), processesLumi.end());
0115 }
0116
0117 bool LHERunInfo::operator==(const LHERunInfo &other) const { return heprup == other.heprup; }
0118
0119 void LHERunInfo::count(int process, CountMode mode, double eventWeight, double brWeight, double matchWeight) {
0120 std::vector<Process>::iterator proc = std::lower_bound(processes.begin(), processes.end(), process);
0121 if (proc == processes.end() || proc->process() != process)
0122 return;
0123
0124 std::vector<Process>::iterator procLumi = std::lower_bound(processesLumi.begin(), processesLumi.end(), process);
0125 if (procLumi == processesLumi.end() || procLumi->process() != process)
0126 return;
0127
0128 switch (mode) {
0129 case kAccepted:
0130 proc->addAcceptedBr(eventWeight * brWeight * matchWeight);
0131 proc->addAccepted(eventWeight * matchWeight);
0132 procLumi->addAcceptedBr(eventWeight * brWeight * matchWeight);
0133 procLumi->addAccepted(eventWeight * matchWeight);
0134 [[fallthrough]];
0135 case kKilled:
0136 proc->addKilled(eventWeight * matchWeight);
0137 procLumi->addKilled(eventWeight * matchWeight);
0138 if (eventWeight > 0) {
0139 proc->addNPassPos();
0140 procLumi->addNPassPos();
0141 } else {
0142 proc->addNPassNeg();
0143 procLumi->addNPassNeg();
0144 }
0145 [[fallthrough]];
0146 case kSelected:
0147 proc->addSelected(eventWeight);
0148 procLumi->addSelected(eventWeight);
0149 if (eventWeight > 0) {
0150 proc->addNTotalPos();
0151 procLumi->addNTotalPos();
0152 } else {
0153 proc->addNTotalNeg();
0154 procLumi->addNTotalNeg();
0155 }
0156 [[fallthrough]];
0157 case kTried:
0158 proc->addTried(eventWeight);
0159 procLumi->addTried(eventWeight);
0160 }
0161 }
0162
0163 LHERunInfo::XSec LHERunInfo::xsec() const {
0164 double sigBrSum = 0.0;
0165 double errBr2Sum = 0.0;
0166 int idwtup = heprup.IDWTUP;
0167 for (std::vector<Process>::const_iterator proc = processes.begin(); proc != processes.end(); ++proc) {
0168 unsigned int idx = proc->heprupIndex();
0169
0170 if (!proc->killed().n())
0171 continue;
0172
0173 double sigma2Sum, sigma2Err;
0174 sigma2Sum = heprup.XSECUP[idx] * heprup.XSECUP[idx];
0175 sigma2Err = heprup.XERRUP[idx] * heprup.XERRUP[idx];
0176
0177 double sigmaAvg = heprup.XSECUP[idx];
0178
0179 double fracAcc = 0;
0180 double ntotal = proc->nTotalPos() - proc->nTotalNeg();
0181 double npass = proc->nPassPos() - proc->nPassNeg();
0182 switch (idwtup) {
0183 case 3:
0184 case -3:
0185 fracAcc = ntotal > 0 ? npass / ntotal : -1;
0186 break;
0187 default:
0188 fracAcc = proc->selected().sum() > 0 ? proc->killed().sum() / proc->selected().sum() : -1;
0189 break;
0190 }
0191
0192 if (fracAcc <= 0)
0193 continue;
0194
0195 double fracBr = proc->accepted().sum() > 0.0 ? proc->acceptedBr().sum() / proc->accepted().sum() : 1;
0196 double sigmaFin = sigmaAvg * fracAcc;
0197 double sigmaFinBr = sigmaFin * fracBr;
0198
0199 double relErr = 1.0;
0200
0201 double efferr2 = 0;
0202 switch (idwtup) {
0203 case 3:
0204 case -3: {
0205 double ntotal_pos = proc->nTotalPos();
0206 double effp = ntotal_pos > 0 ? (double)proc->nPassPos() / ntotal_pos : 0;
0207 double effp_err2 = ntotal_pos > 0 ? (1 - effp) * effp / ntotal_pos : 0;
0208
0209 double ntotal_neg = proc->nTotalNeg();
0210 double effn = ntotal_neg > 0 ? (double)proc->nPassNeg() / ntotal_neg : 0;
0211 double effn_err2 = ntotal_neg > 0 ? (1 - effn) * effn / ntotal_neg : 0;
0212
0213 efferr2 = ntotal > 0
0214 ? (ntotal_pos * ntotal_pos * effp_err2 + ntotal_neg * ntotal_neg * effn_err2) / ntotal / ntotal
0215 : 0;
0216 break;
0217 }
0218 default: {
0219 double denominator = pow(proc->selected().sum(), 4);
0220 double passw = proc->killed().sum();
0221 double passw2 = proc->killed().sum2();
0222 double failw = proc->selected().sum() - passw;
0223 double failw2 = proc->selected().sum2() - passw2;
0224 double numerator = (passw2 * failw * failw + failw2 * passw * passw);
0225
0226 efferr2 = denominator > 0 ? numerator / denominator : 0;
0227 break;
0228 }
0229 }
0230 double delta2Veto = efferr2 / fracAcc / fracAcc;
0231 double delta2Sum = delta2Veto + sigma2Err / sigma2Sum;
0232 relErr = (delta2Sum > 0.0 ? std::sqrt(delta2Sum) : 0.0);
0233
0234 double deltaFinBr = sigmaFinBr * relErr;
0235
0236 sigBrSum += sigmaFinBr;
0237 errBr2Sum += deltaFinBr * deltaFinBr;
0238 }
0239
0240 XSec result(sigBrSum, std::sqrt(errBr2Sum));
0241
0242 return result;
0243 }
0244
0245 void LHERunInfo::statistics() const {
0246 double sigSelSum = 0.0;
0247 double sigSum = 0.0;
0248 double sigBrSum = 0.0;
0249 double errSel2Sum = 0.0;
0250 double errBr2Sum = 0.0;
0251 double errMatch2Sum = 0.0;
0252 unsigned long nAccepted = 0;
0253 unsigned long nTried = 0;
0254 unsigned long nAccepted_pos = 0;
0255 unsigned long nTried_pos = 0;
0256 unsigned long nAccepted_neg = 0;
0257 unsigned long nTried_neg = 0;
0258 int idwtup = heprup.IDWTUP;
0259
0260 LogDebug("LHERunInfo") << " statistics";
0261 LogDebug("LHERunInfo") << "Process and cross-section statistics";
0262 LogDebug("LHERunInfo") << "------------------------------------";
0263 LogDebug("LHERunInfo") << "Process\t\txsec_before [pb]\t\tpassed\tnposw\tnnegw\ttried\tnposw\tnnegw \txsec_match "
0264 "[pb]\t\t\taccepted [%]\t event_eff [%]";
0265
0266 for (std::vector<Process>::const_iterator proc = processes.begin(); proc != processes.end(); ++proc) {
0267 unsigned int idx = proc->heprupIndex();
0268
0269 if (!proc->selected().n()) {
0270 LogDebug("LHERunInfo") << proc->process() << "\t0\t0\tn/a\t\t\tn/a";
0271 continue;
0272 }
0273
0274 double sigma2Sum, sigma2Err;
0275 sigma2Sum = heprup.XSECUP[idx] * heprup.XSECUP[idx];
0276 sigma2Err = heprup.XERRUP[idx] * heprup.XERRUP[idx];
0277
0278 double sigmaAvg = heprup.XSECUP[idx];
0279
0280 double fracAcc = 0;
0281 double ntotal = proc->nTotalPos() - proc->nTotalNeg();
0282 double npass = proc->nPassPos() - proc->nPassNeg();
0283 switch (idwtup) {
0284 case 3:
0285 case -3:
0286 fracAcc = ntotal > 0 ? npass / ntotal : -1;
0287 break;
0288 default:
0289 fracAcc = proc->selected().sum() > 0 ? proc->killed().sum() / proc->selected().sum() : -1;
0290 break;
0291 }
0292
0293 double fracBr = proc->accepted().sum() > 0.0 ? proc->acceptedBr().sum() / proc->accepted().sum() : 1;
0294 double sigmaFin = fracAcc > 0 ? sigmaAvg * fracAcc : 0;
0295 double sigmaFinBr = sigmaFin * fracBr;
0296
0297 double relErr = 1.0;
0298 double relAccErr = 1.0;
0299 double efferr2 = 0;
0300
0301 if (proc->killed().n() > 0 && fracAcc > 0) {
0302 switch (idwtup) {
0303 case 3:
0304 case -3: {
0305 double ntotal_pos = proc->nTotalPos();
0306 double effp = ntotal_pos > 0 ? (double)proc->nPassPos() / ntotal_pos : 0;
0307 double effp_err2 = ntotal_pos > 0 ? (1 - effp) * effp / ntotal_pos : 0;
0308
0309 double ntotal_neg = proc->nTotalNeg();
0310 double effn = ntotal_neg > 0 ? (double)proc->nPassNeg() / ntotal_neg : 0;
0311 double effn_err2 = ntotal_neg > 0 ? (1 - effn) * effn / ntotal_neg : 0;
0312
0313 efferr2 = ntotal > 0 ? (ntotal_pos * ntotal_pos * effp_err2 + ntotal_neg * ntotal_neg * effn_err2) /
0314 ntotal / ntotal
0315 : 0;
0316 break;
0317 }
0318 default: {
0319 double denominator = pow(proc->selected().sum(), 4);
0320 double passw = proc->killed().sum();
0321 double passw2 = proc->killed().sum2();
0322 double failw = proc->selected().sum() - passw;
0323 double failw2 = proc->selected().sum2() - passw2;
0324 double numerator = (passw2 * failw * failw + failw2 * passw * passw);
0325
0326 efferr2 = denominator > 0 ? numerator / denominator : 0;
0327 break;
0328 }
0329 }
0330 double delta2Veto = efferr2 / fracAcc / fracAcc;
0331 double delta2Sum = delta2Veto + sigma2Err / sigma2Sum;
0332 relErr = (delta2Sum > 0.0 ? std::sqrt(delta2Sum) : 0.0);
0333 relAccErr = (delta2Veto > 0.0 ? std::sqrt(delta2Veto) : 0.0);
0334 }
0335 double deltaFinBr = sigmaFinBr * relErr;
0336
0337 double ntotal_proc = proc->nTotalPos() + proc->nTotalNeg();
0338 double event_eff_proc = ntotal_proc > 0 ? (double)(proc->nPassPos() + proc->nPassNeg()) / ntotal_proc : -1;
0339 double event_eff_err_proc = ntotal_proc > 0 ? std::sqrt((1 - event_eff_proc) * event_eff_proc / ntotal_proc) : -1;
0340
0341 LogDebug("LHERunInfo") << proc->process() << "\t\t" << std::scientific << std::setprecision(3)
0342 << heprup.XSECUP[proc->heprupIndex()] << " +/- " << heprup.XERRUP[proc->heprupIndex()]
0343 << "\t\t" << proc->accepted().n() << "\t" << proc->nPassPos() << "\t" << proc->nPassNeg()
0344 << "\t" << proc->tried().n() << "\t" << proc->nTotalPos() << "\t" << proc->nTotalNeg()
0345 << "\t" << std::scientific << std::setprecision(3) << sigmaFinBr << " +/- " << deltaFinBr
0346 << "\t\t" << std::fixed << std::setprecision(1) << (fracAcc * 100) << " +/- "
0347 << (std::sqrt(efferr2) * 100) << "\t" << std::fixed << std::setprecision(1)
0348 << (event_eff_proc * 100) << " +/- " << (event_eff_err_proc * 100);
0349
0350 nAccepted += proc->accepted().n();
0351 nTried += proc->tried().n();
0352 nAccepted_pos += proc->nPassPos();
0353 nTried_pos += proc->nTotalPos();
0354 nAccepted_neg += proc->nPassNeg();
0355 nTried_neg += proc->nTotalNeg();
0356 sigSelSum += sigmaAvg;
0357 sigSum += sigmaFin;
0358 sigBrSum += sigmaFinBr;
0359 errSel2Sum += sigma2Err;
0360 errBr2Sum += deltaFinBr * deltaFinBr;
0361 errMatch2Sum += sigmaFin * relAccErr * sigmaFin * relAccErr;
0362 }
0363
0364 double ntotal_all = (nTried_pos + nTried_neg);
0365 double event_eff_all = ntotal_all > 0 ? (double)(nAccepted_pos + nAccepted_neg) / ntotal_all : -1;
0366 double event_eff_err_all = ntotal_all > 0 ? std::sqrt((1 - event_eff_all) * event_eff_all / ntotal_all) : -1;
0367
0368 LogDebug("LHERunInfo") << "Total\t\t" << std::scientific << std::setprecision(3) << sigSelSum << " +/- "
0369 << std::sqrt(errSel2Sum) << "\t\t" << nAccepted << "\t" << nAccepted_pos << "\t"
0370 << nAccepted_neg << "\t" << nTried << "\t" << nTried_pos << "\t" << nTried_neg << "\t"
0371 << std::scientific << std::setprecision(3) << sigBrSum << " +/- " << std::sqrt(errBr2Sum)
0372 << "\t\t" << std::fixed << std::setprecision(1) << (sigSum / sigSelSum * 100) << " +/- "
0373 << (std::sqrt(errMatch2Sum) / sigSelSum * 100) << "\t" << std::fixed << std::setprecision(1)
0374 << (event_eff_all * 100) << " +/- " << (event_eff_err_all * 100);
0375 }
0376
0377 LHERunInfo::Header::Header() : xmlDoc(nullptr) {}
0378
0379 LHERunInfo::Header::Header(const std::string &tag) : LHERunInfoProduct::Header(tag), xmlDoc(nullptr) {}
0380
0381 LHERunInfo::Header::Header(const Header &orig) : LHERunInfoProduct::Header(orig), xmlDoc(nullptr) {}
0382
0383 LHERunInfo::Header::Header(const LHERunInfoProduct::Header &orig)
0384 : LHERunInfoProduct::Header(orig), xmlDoc(nullptr) {}
0385
0386 LHERunInfo::Header::~Header() {
0387 if (xmlDoc)
0388 xmlDoc->release();
0389 }
0390
0391 static void fillLines(std::vector<std::string> &lines, const char *data, int len = -1) {
0392 const char *end = len >= 0 ? (data + len) : nullptr;
0393 while (*data && (!end || data < end)) {
0394 std::size_t len = std::strcspn(data, "\r\n");
0395 if (end && data + len > end)
0396 len = end - data;
0397 if (data[len] == '\r' && data[len + 1] == '\n')
0398 len += 2;
0399 else if (data[len])
0400 len++;
0401 lines.push_back(std::string(data, len));
0402 data += len;
0403 }
0404 }
0405
0406 static std::vector<std::string> domToLines(const DOMNode *node) {
0407 std::vector<std::string> result;
0408 DOMImplementation *impl = DOMImplementationRegistry::getDOMImplementation(XMLUniStr("Core"));
0409 std::unique_ptr<DOMLSSerializer> writer(((DOMImplementationLS *)(impl))->createLSSerializer());
0410
0411 std::unique_ptr<DOMLSOutput> outputDesc(((DOMImplementationLS *)impl)->createLSOutput());
0412 assert(outputDesc.get());
0413 outputDesc->setEncoding(XMLUniStr("UTF-8"));
0414
0415 XMLSimpleStr buffer(writer->writeToString(node));
0416
0417 const char *p = std::strchr((const char *)buffer, '>') + 1;
0418 const char *q = std::strrchr(p, '<');
0419 fillLines(result, p, q - p);
0420
0421 return result;
0422 }
0423
0424 std::vector<std::string> LHERunInfo::findHeader(const std::string &tag) const {
0425 const LHERunInfo::Header *header = nullptr;
0426 for (std::vector<Header>::const_iterator iter = headers.begin(); iter != headers.end(); ++iter) {
0427 if (iter->tag() == tag)
0428 return std::vector<std::string>(iter->begin(), iter->end());
0429 if (iter->tag() == "header")
0430 header = &*iter;
0431 }
0432
0433 if (!header)
0434 return std::vector<std::string>();
0435
0436 const DOMNode *root = header->getXMLNode();
0437 if (!root)
0438 return std::vector<std::string>();
0439
0440 for (const DOMNode *iter = root->getFirstChild(); iter; iter = iter->getNextSibling()) {
0441 if (iter->getNodeType() != DOMNode::ELEMENT_NODE)
0442 continue;
0443 if (tag == (const char *)XMLSimpleStr(iter->getNodeName()))
0444 return domToLines(iter);
0445 }
0446
0447 return std::vector<std::string>();
0448 }
0449
0450 namespace {
0451 class HeaderReader : public CBInputStream::Reader {
0452 public:
0453 HeaderReader(const LHERunInfo::Header *header) : header(header), mode(kHeader), iter(header->begin()) {}
0454
0455 const std::string &data() override {
0456 switch (mode) {
0457 case kHeader:
0458 tmp = "<" + header->tag() + ">";
0459 mode = kBody;
0460 break;
0461 case kBody:
0462 if (iter != header->end())
0463 return *iter++;
0464 tmp = "</" + header->tag() + ">";
0465 mode = kFooter;
0466 break;
0467 case kFooter:
0468 tmp.clear();
0469 }
0470
0471 return tmp;
0472 }
0473
0474 private:
0475 enum Mode { kHeader, kBody, kFooter };
0476
0477 const LHERunInfo::Header *header;
0478 Mode mode;
0479 LHERunInfo::Header::const_iterator iter;
0480 std::string tmp;
0481 };
0482 }
0483
0484 const DOMNode *LHERunInfo::Header::getXMLNode() const {
0485 if (tag().empty())
0486 return nullptr;
0487
0488 if (!xmlDoc) {
0489 XercesDOMParser parser;
0490 parser.setValidationScheme(XercesDOMParser::Val_Auto);
0491 parser.setDoNamespaces(false);
0492 parser.setDoSchema(false);
0493 parser.setValidationSchemaFullChecking(false);
0494
0495 HandlerBase errHandler;
0496 parser.setErrorHandler(&errHandler);
0497 parser.setCreateEntityReferenceNodes(false);
0498
0499 try {
0500 std::unique_ptr<CBInputStream::Reader> reader(new HeaderReader(this));
0501 CBInputSource source(reader);
0502 parser.parse(source);
0503 xmlDoc = parser.adoptDocument();
0504 } catch (const XMLException &e) {
0505 throw cms::Exception("Generator|LHEInterface")
0506 << "XML parser reported DOM error no. " << (unsigned long)e.getCode() << ": "
0507 << XMLSimpleStr(e.getMessage()) << "." << std::endl;
0508 } catch (const SAXException &e) {
0509 throw cms::Exception("Generator|LHEInterface")
0510 << "XML parser reported: " << XMLSimpleStr(e.getMessage()) << "." << std::endl;
0511 }
0512 }
0513
0514 return xmlDoc->getDocumentElement();
0515 }
0516
0517 std::pair<int, int> LHERunInfo::pdfSetTranslation() const {
0518 int pdfA = -1, pdfB = -1;
0519
0520 if (heprup.PDFGUP.first >= 0) {
0521 pdfA = heprup.PDFSUP.first;
0522 }
0523
0524 if (heprup.PDFGUP.second >= 0) {
0525 pdfB = heprup.PDFSUP.second;
0526 }
0527
0528 return std::make_pair(pdfA, pdfB);
0529 }
0530
0531 const bool operator==(const LHERunInfo::Process &lhs, const LHERunInfo::Process &rhs) {
0532 return (lhs.process() == rhs.process());
0533 }
0534
0535 const bool operator<(const LHERunInfo::Process &lhs, const LHERunInfo::Process &rhs) {
0536 return (lhs.process() < rhs.process());
0537 }
0538
0539 }
