Back to home page

Project CMSSW displayed by LXR

 
 

    


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   }  // anonymous namespace
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 }  // namespace lhef