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