Core/plugins/GenXSecAnalyzer.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "TMath.h"
0003 #include <iostream>
0004 #include <iomanip>
0005
0006 // analyzer of a summary information product on filter efficiency for a user specified path
0007 // meant for the generator filter efficiency calculation
0008
0009 // system include files
0010 #include <memory>
0011 #include <vector>
0012 #include <map>
0013
0014 // user include files
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/global/EDAnalyzer.h"
0017
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/Run.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/LuminosityBlock.h"
0022 #include "FWCore/Utilities/interface/InputTag.h"
0023 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoProduct.h"
0026 #include "SimDataFormats/GeneratorProducts/interface/GenFilterInfo.h"
0027 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0028
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030
0031 //
0032 // class declaration
0033 //
0034 namespace gxsec {
0035   struct LumiCache {};
0036   struct RunCache {
0037     RunCache()
0038         : product_(-9999),
0039           filterOnlyEffRun_(0, 0, 0, 0, 0., 0., 0., 0.),
0040           hepMCFilterEffRun_(0, 0, 0, 0, 0., 0., 0., 0.) {}
0041     // for weight before GenFilter and HepMCFilter and before matching
0042     CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable double thisRunWeightPre_ = 0;
0043
0044     // for weight after GenFilter and HepMCFilter and after matching
0045     CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable double thisRunWeight_ = 0;
0046
0047     // GenLumiInfo before HepMCFilter and GenFilter, this is used
0048     // for computation
0049     CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable GenLumiInfoProduct product_;
0050
0051     // statistics from additional generator filter, for computation
0052     // reset for each run
0053     CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable GenFilterInfo filterOnlyEffRun_;
0054
0055     // statistics from HepMC filter, for computation
0056     CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable GenFilterInfo hepMCFilterEffRun_;
0057
0058     // the following vectors all have the same size
0059     // LHE or Pythia/Herwig cross section of previous luminosity block
0060     // vector size = number of processes, used for computation
0061     CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable std::map<int, GenLumiInfoProduct::XSec> previousLumiBlockLHEXSec_;
0062
0063     // LHE or Pythia/Herwig combined cross section of current luminosity block
0064     // updated for each luminosity block, initialized in every run
0065     // used for computation
0066     CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable std::map<int, GenLumiInfoProduct::XSec> currentLumiBlockLHEXSec_;
0067   };
0068 }  // namespace gxsec
0069
0070 class GenXSecAnalyzer
0071     : public edm::global::EDAnalyzer<edm::RunCache<gxsec::RunCache>, edm::LuminosityBlockCache<gxsec::LumiCache>> {
0072 public:
0073   explicit GenXSecAnalyzer(const edm::ParameterSet &);
0074   ~GenXSecAnalyzer() override;
0075
0076 private:
0077   void beginJob() final;
0078   std::shared_ptr<gxsec::RunCache> globalBeginRun(edm::Run const &, edm::EventSetup const &) const final;
0079   std::shared_ptr<gxsec::LumiCache> globalBeginLuminosityBlock(edm::LuminosityBlock const &,
0080                                                                edm::EventSetup const &) const final;
0081   void analyze(edm::StreamID, const edm::Event &, const edm::EventSetup &) const final;
0082   void globalEndLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) const final;
0083   void globalEndRun(edm::Run const &, edm::EventSetup const &) const final;
0084   void endJob() final;
0085   // computation of cross section after matching and before HepcFilter and GenFilter
0086   GenLumiInfoProduct::XSec compute(const GenLumiInfoProduct &) const;
0087   // combination of cross section from different MCs after matching (could be either before or after HepcFilter and GenFilter)
0088   void combine(GenLumiInfoProduct::XSec &, double &, const GenLumiInfoProduct::XSec &, const double &) const;
0089   void combine(double &, double &, double &, const double &, const double &, const double &) const;
0090
0091   edm::EDGetTokenT<GenFilterInfo> genFilterInfoToken_;
0092   edm::EDGetTokenT<GenFilterInfo> hepMCFilterInfoToken_;
0093   edm::EDGetTokenT<GenLumiInfoProduct> genLumiInfoToken_;
0094   edm::EDGetTokenT<LHERunInfoProduct> lheRunInfoToken_;
0095
0096   // ----------member data --------------------------
0097
0098   mutable std::atomic<int> nMCs_;
0099
0100   mutable std::atomic<int> hepidwtup_;
0101
0102   mutable std::mutex mutex_;
0103
0104   // for weight before GenFilter and HepMCFilter and before matching
0105   CMS_THREAD_GUARD(mutex_) mutable double totalWeightPre_;
0106
0107   // for weight after GenFilter and HepMCFilter and after matching
0108   CMS_THREAD_GUARD(mutex_) mutable double totalWeight_;
0109
0110   // combined cross sections before HepMCFilter and GenFilter
0111   CMS_THREAD_GUARD(mutex_) mutable GenLumiInfoProduct::XSec xsecPreFilter_;
0112
0113   // final combined cross sections
0114   CMS_THREAD_GUARD(mutex_) mutable GenLumiInfoProduct::XSec xsec_;
0115
0116   // statistics from additional generator filter, for print-out only
0117   CMS_THREAD_GUARD(mutex_) mutable GenFilterInfo filterOnlyEffStat_;
0118
0119   // statistics from HepMC filter, for print-out only
0120   CMS_THREAD_GUARD(mutex_) mutable GenFilterInfo hepMCFilterEffStat_;
0121
0122   // the vector/map size is the number of LHE processes + 1
0123   // needed only for printouts, not used for computation
0124   // only printed out when combining the same physics process
0125   // uncertainty-averaged cross sections before matching
0126   CMS_THREAD_GUARD(mutex_) mutable std::vector<GenLumiInfoProduct::XSec> xsecBeforeMatching_;
0127   // uncertainty-averaged cross sections after matching
0128   CMS_THREAD_GUARD(mutex_) mutable std::vector<GenLumiInfoProduct::XSec> xsecAfterMatching_;
0129   // statistics from jet matching
0130   CMS_THREAD_GUARD(mutex_) mutable std::map<int, GenFilterInfo> jetMatchEffStat_;
0131 };
0132
0133 GenXSecAnalyzer::GenXSecAnalyzer(const edm::ParameterSet &iConfig)
0134     : nMCs_(0),
0135       hepidwtup_(-9999),
0136       totalWeightPre_(0),
0137       totalWeight_(0),
0138       xsecPreFilter_(-1, -1),
0139       xsec_(-1, -1),
0140       filterOnlyEffStat_(0, 0, 0, 0, 0., 0., 0., 0.),
0141       hepMCFilterEffStat_(0, 0, 0, 0, 0., 0., 0., 0.) {
0142   genFilterInfoToken_ = consumes<GenFilterInfo, edm::InLumi>(edm::InputTag("genFilterEfficiencyProducer", ""));
0143   hepMCFilterInfoToken_ = consumes<GenFilterInfo, edm::InLumi>(edm::InputTag("generator", ""));
0144   genLumiInfoToken_ = consumes<GenLumiInfoProduct, edm::InLumi>(edm::InputTag("generator", ""));
0145   lheRunInfoToken_ = consumes<LHERunInfoProduct, edm::InRun>(edm::InputTag("externalLHEProducer", ""));
0146 }
0147
0148 GenXSecAnalyzer::~GenXSecAnalyzer() {}
0149
0150 void GenXSecAnalyzer::beginJob() {}
0151
0152 std::shared_ptr<gxsec::RunCache> GenXSecAnalyzer::globalBeginRun(edm::Run const &iRun, edm::EventSetup const &) const {
0153   // initialization for every different physics MC
0154
0155   nMCs_++;
0156
0157   {
0158     std::lock_guard l{mutex_};
0159     xsecBeforeMatching_.clear();
0160     xsecAfterMatching_.clear();
0161     jetMatchEffStat_.clear();
0162   }
0163   return std::make_shared<gxsec::RunCache>();
0164 }
0165
0166 std::shared_ptr<gxsec::LumiCache> GenXSecAnalyzer::globalBeginLuminosityBlock(edm::LuminosityBlock const &iLumi,
0167                                                                               edm::EventSetup const &) const {
0168   return std::shared_ptr<gxsec::LumiCache>();
0169 }
0170
0171 void GenXSecAnalyzer::analyze(edm::StreamID, const edm::Event &, const edm::EventSetup &) const {}
0172
0173 void GenXSecAnalyzer::globalEndLuminosityBlock(edm::LuminosityBlock const &iLumi, edm::EventSetup const &) const {
0174   edm::Handle<GenLumiInfoProduct> genLumiInfo;
0175   iLumi.getByToken(genLumiInfoToken_, genLumiInfo);
0176   if (!genLumiInfo.isValid())
0177     return;
0178   hepidwtup_ = genLumiInfo->getHEPIDWTUP();
0179
0180   std::vector<GenLumiInfoProduct::ProcessInfo> const &theProcesses = genLumiInfo->getProcessInfos();
0181
0182   unsigned int theProcesses_size = theProcesses.size();
0183
0184   // if it's a pure parton-shower generator, check there should be only one element in thisProcessInfos
0185   // the error of lheXSec is -1
0186   if (hepidwtup_ == -1) {
0187     if (theProcesses_size != 1) {
0188       edm::LogError("GenXSecAnalyzer::endLuminosityBlock") << "Pure parton shower has thisProcessInfos size!=1";
0189       return;
0190     }
0191   }
0192
0193   for (unsigned int ip = 0; ip < theProcesses_size; ip++) {
0194     if (theProcesses[ip].lheXSec().value() < 0) {
0195       edm::LogError("GenXSecAnalyzer::endLuminosityBlock")
0196           << "cross section of process " << ip << " value = " << theProcesses[ip].lheXSec().value();
0197       return;
0198     }
0199   }
0200
0201   auto runC = runCache(iLumi.getRun().index());
0202   {
0203     std::lock_guard g{mutex_};
0204     runC->product_.mergeProduct(*genLumiInfo);
0205   }
0206   edm::Handle<GenFilterInfo> genFilter;
0207   iLumi.getByToken(genFilterInfoToken_, genFilter);
0208
0209   if (genFilter.isValid()) {
0210     std::lock_guard g{mutex_};
0211     filterOnlyEffStat_.mergeProduct(*genFilter);
0212     runC->filterOnlyEffRun_.mergeProduct(*genFilter);
0213     runC->thisRunWeight_ += genFilter->sumPassWeights();
0214   }
0215
0216   edm::Handle<GenFilterInfo> hepMCFilter;
0217   iLumi.getByToken(hepMCFilterInfoToken_, hepMCFilter);
0218
0219   if (hepMCFilter.isValid()) {
0220     std::lock_guard g{mutex_};
0221     hepMCFilterEffStat_.mergeProduct(*hepMCFilter);
0222     runC->hepMCFilterEffRun_.mergeProduct(*hepMCFilter);
0223   }
0224
0225   std::lock_guard g{mutex_};
0226   // doing generic summing for jet matching statistics
0227   // and computation of combined LHE information
0228   for (unsigned int ip = 0; ip < theProcesses_size; ip++) {
0229     int id = theProcesses[ip].process();
0230     GenFilterInfo &x = jetMatchEffStat_[id];
0231     GenLumiInfoProduct::XSec &y = runC->currentLumiBlockLHEXSec_[id];
0232     GenLumiInfoProduct::FinalStat temp_killed = theProcesses[ip].killed();
0233     GenLumiInfoProduct::FinalStat temp_selected = theProcesses[ip].selected();
0234     double passw = temp_killed.sum();
0235     double passw2 = temp_killed.sum2();
0236     double totalw = temp_selected.sum();
0237     double totalw2 = temp_selected.sum2();
0238     GenFilterInfo tempInfo(theProcesses[ip].nPassPos(),
0239                            theProcesses[ip].nPassNeg(),
0240                            theProcesses[ip].nTotalPos(),
0241                            theProcesses[ip].nTotalNeg(),
0242                            passw,
0243                            passw2,
0244                            totalw,
0245                            totalw2);
0246
0247     // matching statistics for all processes
0248     jetMatchEffStat_[10000].mergeProduct(tempInfo);
0249     double currentValue = theProcesses[ip].lheXSec().value();
0250     double currentError = theProcesses[ip].lheXSec().error();
0251
0252     // this process ID has occurred before
0253     auto &thisRunWeightPre = runC->thisRunWeightPre_;
0254     if (y.value() > 0) {
0255       x.mergeProduct(tempInfo);
0256       double previousValue = runC->previousLumiBlockLHEXSec_[id].value();
0257
0258       if (currentValue != previousValue)  // transition of cross section
0259       {
0260         double xsec = y.value();
0261         double err = y.error();
0262         combine(xsec, err, thisRunWeightPre, currentValue, currentError, totalw);
0263         y = GenLumiInfoProduct::XSec(xsec, err);
0264       } else  // LHE cross section is the same as previous lumiblock
0265         thisRunWeightPre += totalw;
0266
0267     }
0268     // this process ID has never occurred before
0269     else {
0270       x = tempInfo;
0271       y = theProcesses[ip].lheXSec();
0272       thisRunWeightPre += totalw;
0273     }
0274
0275     runC->previousLumiBlockLHEXSec_[id] = theProcesses[ip].lheXSec();
0276   }  // end
0277
0278   return;
0279 }
0280
0281 void GenXSecAnalyzer::globalEndRun(edm::Run const &iRun, edm::EventSetup const &) const {
0282   //xsection before matching
0283   edm::Handle<LHERunInfoProduct> run;
0284
0285   if (iRun.getByToken(lheRunInfoToken_, run)) {
0286     const lhef::HEPRUP thisHeprup = run->heprup();
0287
0288     for (unsigned int iSize = 0; iSize < thisHeprup.XSECUP.size(); iSize++) {
0289       std::cout << std::setw(14) << std::fixed << thisHeprup.XSECUP[iSize] << std::setw(14) << std::fixed
0290                 << thisHeprup.XERRUP[iSize] << std::setw(14) << std::fixed << thisHeprup.XMAXUP[iSize] << std::setw(14)
0291                 << std::fixed << thisHeprup.LPRUP[iSize] << std::endl;
0292     }
0293     std::cout << " " << std::endl;
0294   }
0295
0296   auto runC = runCache(iRun.index());
0297   std::lock_guard l{mutex_};
0298
0299   // compute cross section for this run first
0300   // set the correct combined LHE+filter cross sections
0301   unsigned int i = 0;
0302   std::vector<GenLumiInfoProduct::ProcessInfo> newInfos;
0303   for (std::map<int, GenLumiInfoProduct::XSec>::const_iterator iter = runC->currentLumiBlockLHEXSec_.begin();
0304        iter != runC->currentLumiBlockLHEXSec_.end();
0305        ++iter, i++) {
0306     GenLumiInfoProduct::ProcessInfo temp = runC->product_.getProcessInfos()[i];
0307     temp.setLheXSec(iter->second.value(), iter->second.error());
0308     newInfos.push_back(temp);
0309   }
0310   runC->product_.setProcessInfo(newInfos);
0311
0312   const GenLumiInfoProduct::XSec thisRunXSecPre = compute(runC->product_);
0313   // xsection after matching before filters
0314   combine(xsecPreFilter_, totalWeightPre_, thisRunXSecPre, runC->thisRunWeightPre_);
0315
0316   double thisHepFilterEff = 1;
0317   double thisHepFilterErr = 0;
0318
0319   if (runC->hepMCFilterEffRun_.sumWeights2() > 0) {
0320     thisHepFilterEff = runC->hepMCFilterEffRun_.filterEfficiency(hepidwtup_);
0321     thisHepFilterErr = runC->hepMCFilterEffRun_.filterEfficiencyError(hepidwtup_);
0322     if (thisHepFilterEff < 0) {
0323       thisHepFilterEff = 1;
0324       thisHepFilterErr = 0;
0325     }
0326   }
0327
0328   double thisGenFilterEff = 1;
0329   double thisGenFilterErr = 0;
0330
0331   if (runC->filterOnlyEffRun_.sumWeights2() > 0) {
0332     thisGenFilterEff = runC->filterOnlyEffRun_.filterEfficiency(hepidwtup_);
0333     thisGenFilterErr = runC->filterOnlyEffRun_.filterEfficiencyError(hepidwtup_);
0334     if (thisGenFilterEff < 0) {
0335       thisGenFilterEff = 1;
0336       thisGenFilterErr = 0;
0337     }
0338   }
0339   double thisXsec = thisRunXSecPre.value() > 0 ? thisHepFilterEff * thisGenFilterEff * thisRunXSecPre.value() : 0;
0340   double thisErr =
0341       thisRunXSecPre.value() > 0
0342           ? thisXsec * sqrt(pow(TMath::Max(thisRunXSecPre.error(), (double)0) / thisRunXSecPre.value(), 2) +
0343                             pow(thisHepFilterErr / thisHepFilterEff, 2) + pow(thisGenFilterErr / thisGenFilterEff, 2))
0344           : 0;
0345   const GenLumiInfoProduct::XSec thisRunXSec = GenLumiInfoProduct::XSec(thisXsec, thisErr);
0346   combine(xsec_, totalWeight_, thisRunXSec, runC->thisRunWeight_);
0347 }
0348
0349 void GenXSecAnalyzer::combine(double &finalValue,
0350                               double &finalError,
0351                               double &finalWeight,
0352                               const double &currentValue,
0353                               const double &currentError,
0354                               const double &currentWeight) const {
0355   if (finalValue <= 0) {
0356     finalValue = currentValue;
0357     finalError = currentError;
0358     finalWeight += currentWeight;
0359   } else {
0360     double wgt1 = (finalError <= 0 || currentError <= 0) ? finalWeight : 1 / (finalError * finalError);
0361     double wgt2 = (finalError <= 0 || currentError <= 0) ? currentWeight : 1 / (currentError * currentError);
0362     double xsec = (wgt1 * finalValue + wgt2 * currentValue) / (wgt1 + wgt2);
0363     double err = (finalError <= 0 || currentError <= 0) ? 0 : 1.0 / std::sqrt(wgt1 + wgt2);
0364     finalValue = xsec;
0365     finalError = err;
0366     finalWeight += currentWeight;
0367   }
0368   return;
0369 }
0370
0371 void GenXSecAnalyzer::combine(GenLumiInfoProduct::XSec &finalXSec,
0372                               double &totalw,
0373                               const GenLumiInfoProduct::XSec &thisRunXSec,
0374                               const double &thisw) const {
0375   double value = finalXSec.value();
0376   double error = finalXSec.error();
0377   double thisValue = thisRunXSec.value();
0378   double thisError = thisRunXSec.error();
0379   combine(value, error, totalw, thisValue, thisError, thisw);
0380   finalXSec = GenLumiInfoProduct::XSec(value, error);
0381   return;
0382 }
0383
0384 GenLumiInfoProduct::XSec GenXSecAnalyzer::compute(const GenLumiInfoProduct &iLumiInfo) const {
0385   // sum of cross sections and errors over different processes
0386   double sigSelSum = 0.0;
0387   double err2SelSum = 0.0;
0388
0389   std::vector<GenLumiInfoProduct::XSec> tempVector_before;
0390   std::vector<GenLumiInfoProduct::XSec> tempVector_after;
0391
0392   // loop over different processes for each sample
0393   unsigned int vectorSize = iLumiInfo.getProcessInfos().size();
0394   for (unsigned int ip = 0; ip < vectorSize; ip++) {
0395     GenLumiInfoProduct::ProcessInfo proc = iLumiInfo.getProcessInfos()[ip];
0396     double hepxsec_value = proc.lheXSec().value();
0397     double hepxsec_error = proc.lheXSec().error() <= 0 ? 0 : proc.lheXSec().error();
0398     tempVector_before.push_back(GenLumiInfoProduct::XSec(hepxsec_value, hepxsec_error));
0399
0400     sigSelSum += hepxsec_value;
0401     err2SelSum += hepxsec_error * hepxsec_error;
0402
0403     // skips computation if jet matching efficiency=0
0404     if (proc.killed().n() < 1) {
0405       tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
0406       continue;
0407     }
0408
0409     // computing jet matching efficiency for this process
0410     double fracAcc = 0;
0411     double ntotal = proc.nTotalPos() - proc.nTotalNeg();
0412     double npass = proc.nPassPos() - proc.nPassNeg();
0413     switch (hepidwtup_) {
0414       case 3:
0415       case -3:
0416         fracAcc = ntotal > 0 ? npass / ntotal : -1;
0417         break;
0418       default:
0419         fracAcc = proc.selected().sum() > 0 ? proc.killed().sum() / proc.selected().sum() : -1;
0420         break;
0421     }
0422
0423     if (fracAcc <= 0) {
0424       tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
0425       continue;
0426     }
0427
0428     // cross section after matching for this particular process
0429     double sigmaFin = hepxsec_value * fracAcc;
0430
0431     // computing error on jet matching efficiency
0432     double relErr = 1.0;
0433     double efferr2 = 0;
0434     switch (hepidwtup_) {
0435       case 3:
0436       case -3: {
0437         double ntotal_pos = proc.nTotalPos();
0438         double effp = ntotal_pos > 0 ? (double)proc.nPassPos() / ntotal_pos : 0;
0439         double effp_err2 = ntotal_pos > 0 ? (1 - effp) * effp / ntotal_pos : 0;
0440
0441         double ntotal_neg = proc.nTotalNeg();
0442         double effn = ntotal_neg > 0 ? (double)proc.nPassNeg() / ntotal_neg : 0;
0443         double effn_err2 = ntotal_neg > 0 ? (1 - effn) * effn / ntotal_neg : 0;
0444
0445         efferr2 = ntotal > 0
0446                       ? (ntotal_pos * ntotal_pos * effp_err2 + ntotal_neg * ntotal_neg * effn_err2) / ntotal / ntotal
0447                       : 0;
0448         break;
0449       }
0450       default: {
0451         double denominator = pow(proc.selected().sum(), 4);
0452         double passw = proc.killed().sum();
0453         double passw2 = proc.killed().sum2();
0454         double failw = proc.selected().sum() - passw;
0455         double failw2 = proc.selected().sum2() - passw2;
0456         double numerator = (passw2 * failw * failw + failw2 * passw * passw);
0457
0458         efferr2 = denominator > 0 ? numerator / denominator : 0;
0459         break;
0460       }
0461     }
0462     double delta2Veto = efferr2 / fracAcc / fracAcc;
0463
0464     // computing total error on cross section after matching efficiency
0465
0466     double sigma2Sum, sigma2Err;
0467     sigma2Sum = hepxsec_value * hepxsec_value;
0468     sigma2Err = hepxsec_error * hepxsec_error;
0469
0470     double delta2Sum = delta2Veto + sigma2Err / sigma2Sum;
0471     relErr = (delta2Sum > 0.0 ? std::sqrt(delta2Sum) : 0.0);
0472     double deltaFin = sigmaFin * relErr;
0473
0474     tempVector_after.push_back(GenLumiInfoProduct::XSec(sigmaFin, deltaFin));
0475
0476   }  // end of loop over different processes
0477   tempVector_before.push_back(GenLumiInfoProduct::XSec(sigSelSum, sqrt(err2SelSum)));
0478
0479   double total_matcheff = jetMatchEffStat_[10000].filterEfficiency(hepidwtup_);
0480   double total_matcherr = jetMatchEffStat_[10000].filterEfficiencyError(hepidwtup_);
0481
0482   double xsec_after = sigSelSum * total_matcheff;
0483   double xsecerr_after = (total_matcheff > 0 && sigSelSum > 0)
0484                              ? xsec_after * sqrt(err2SelSum / sigSelSum / sigSelSum +
0485                                                  total_matcherr * total_matcherr / total_matcheff / total_matcheff)
0486                              : 0;
0487
0488   GenLumiInfoProduct::XSec result(xsec_after, xsecerr_after);
0489   tempVector_after.push_back(result);
0490
0491   xsecBeforeMatching_ = tempVector_before;
0492   xsecAfterMatching_ = tempVector_after;
0493
0494   return result;
0495 }
0496
0497 void GenXSecAnalyzer::endJob() {
0498   edm::LogPrint("GenXSecAnalyzer") << "\n"
0499                                    << "------------------------------------"
0500                                    << "\n"
0501                                    << "GenXsecAnalyzer:"
0502                                    << "\n"
0503                                    << "------------------------------------";
0504
0505   if (jetMatchEffStat_.empty()) {
0506     edm::LogPrint("GenXSecAnalyzer") << "------------------------------------"
0507                                      << "\n"
0508                                      << "Cross-section summary not available"
0509                                      << "\n"
0510                                      << "------------------------------------";
0511     return;
0512   }
0513
0514   // fraction of negative weights
0515   double final_fract_neg_w = 0;
0516   double final_fract_neg_w_unc = 0;
0517
0518   // below print out is only for combination of same physics MC samples and ME+Pythia MCs
0519
0520   if (nMCs_ == 1 && hepidwtup_ != -1) {
0521     edm::LogPrint("GenXSecAnalyzer")
0522         << "-----------------------------------------------------------------------------------------------------------"
0523            "--------------------------------------------------------------- \n"
0524         << "Overall cross-section summary \n"
0525         << "-----------------------------------------------------------------------------------------------------------"
0526            "---------------------------------------------------------------";
0527     edm::LogPrint("GenXSecAnalyzer") << "Process\t\txsec_before [pb]\t\tpassed\tnposw\tnnegw\ttried\tnposw\tnnegw "
0528                                         "\txsec_match [pb]\t\t\taccepted [%]\t event_eff [%]";
0529
0530     const unsigned sizeOfInfos = jetMatchEffStat_.size();
0531     const unsigned last = sizeOfInfos - 1;
0532     std::string *title = new std::string[sizeOfInfos];
0533     unsigned int i = 0;
0534     double jetmatch_eff = 0;
0535     double jetmatch_err = 0;
0536     double matching_eff = 1;
0537     double matching_efferr = 1;
0538
0539     for (std::map<int, GenFilterInfo>::const_iterator iter = jetMatchEffStat_.begin(); iter != jetMatchEffStat_.end();
0540          ++iter, i++) {
0541       GenFilterInfo thisJetMatchStat = iter->second;
0542       GenFilterInfo thisEventEffStat =
0543           GenFilterInfo(thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
0544                         0,
0545                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents(),
0546                         0,
0547                         thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
0548                         thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
0549                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents(),
0550                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents());
0551
0552       jetmatch_eff = thisJetMatchStat.filterEfficiency(hepidwtup_);
0553       jetmatch_err = thisJetMatchStat.filterEfficiencyError(hepidwtup_);
0554
0555       if (i == last) {
0556         title[i] = "Total";
0557
0558         edm::LogPrint("GenXSecAnalyzer")
0559             << "-------------------------------------------------------------------------------------------------------"
0560                "------------------------------------------------------------------- ";
0561
0562         // fill negative fraction of negative weights and uncertainty after matching
0563         final_fract_neg_w = thisEventEffStat.numEventsPassed() > 0
0564                                 ? thisJetMatchStat.numPassNegativeEvents() / thisEventEffStat.numEventsPassed()
0565                                 : 0;
0566         final_fract_neg_w_unc =
0567             thisJetMatchStat.numPassNegativeEvents() > 0
0568                 ? final_fract_neg_w * final_fract_neg_w / thisEventEffStat.numEventsPassed() *
0569                       sqrt(thisJetMatchStat.numPassPositiveEvents() * thisJetMatchStat.numPassPositiveEvents() /
0570                                thisJetMatchStat.numPassNegativeEvents() +
0571                            thisJetMatchStat.numPassPositiveEvents())
0572                 : 0;
0573       } else {
0574         title[i] = Form("%d", i);
0575       }
0576
0577       edm::LogPrint("GenXSecAnalyzer") << title[i] << "\t\t" << std::scientific << std::setprecision(3)
0578                                        << xsecBeforeMatching_[i].value() << " +/- " << xsecBeforeMatching_[i].error()
0579                                        << "\t\t" << thisEventEffStat.numEventsPassed() << "\t"
0580                                        << thisJetMatchStat.numPassPositiveEvents() << "\t"
0581                                        << thisJetMatchStat.numPassNegativeEvents() << "\t"
0582                                        << thisEventEffStat.numEventsTotal() << "\t"
0583                                        << thisJetMatchStat.numTotalPositiveEvents() << "\t"
0584                                        << thisJetMatchStat.numTotalNegativeEvents() << "\t" << std::scientific
0585                                        << std::setprecision(3) << xsecAfterMatching_[i].value() << " +/- "
0586                                        << xsecAfterMatching_[i].error() << "\t\t" << std::fixed << std::setprecision(1)
0587                                        << (jetmatch_eff * 100) << " +/- " << (jetmatch_err * 100) << "\t" << std::fixed
0588                                        << std::setprecision(1) << (thisEventEffStat.filterEfficiency(+3) * 100)
0589                                        << " +/- " << (thisEventEffStat.filterEfficiencyError(+3) * 100);
0590
0591       matching_eff = thisEventEffStat.filterEfficiency(+3);
0592       matching_efferr = thisEventEffStat.filterEfficiencyError(+3);
0593     }
0594     delete[] title;
0595
0596     edm::LogPrint("GenXSecAnalyzer")
0597         << "-----------------------------------------------------------------------------------------------------------"
0598            "---------------------------------------------------------------";
0599
0600     edm::LogPrint("GenXSecAnalyzer") << "Before matching: total cross section = " << std::scientific
0601                                      << std::setprecision(3) << xsecBeforeMatching_[last].value() << " +- "
0602                                      << xsecBeforeMatching_[last].error() << " pb";
0603
0604     edm::LogPrint("GenXSecAnalyzer") << "After matching: total cross section = " << std::scientific
0605                                      << std::setprecision(3) << xsecAfterMatching_[last].value() << " +- "
0606                                      << xsecAfterMatching_[last].error() << " pb";
0607
0608     edm::LogPrint("GenXSecAnalyzer") << "Matching efficiency = " << std::fixed << std::setprecision(1) << matching_eff
0609                                      << " +/- " << matching_efferr << "   [TO BE USED IN MCM]";
0610
0611   } else if (hepidwtup_ == -1)
0612     edm::LogPrint("GenXSecAnalyzer") << "Before Filter: total cross section = " << std::scientific
0613                                      << std::setprecision(3) << xsecPreFilter_.value() << " +- "
0614                                      << xsecPreFilter_.error() << " pb";
0615
0616   // hepMC filter efficiency
0617   double hepMCFilter_eff = 1.0;
0618   double hepMCFilter_err = 0.0;
0619   if (hepMCFilterEffStat_.sumWeights2() > 0) {
0620     hepMCFilter_eff = hepMCFilterEffStat_.filterEfficiency(-1);
0621     hepMCFilter_err = hepMCFilterEffStat_.filterEfficiencyError(-1);
0622     edm::LogPrint("GenXSecAnalyzer") << "HepMC filter efficiency (taking into account weights)= "
0623                                      << "(" << hepMCFilterEffStat_.sumPassWeights() << ")"
0624                                      << " / "
0625                                      << "(" << hepMCFilterEffStat_.sumWeights() << ")"
0626                                      << " = " << std::scientific << std::setprecision(3) << hepMCFilter_eff << " +- "
0627                                      << hepMCFilter_err;
0628
0629     double hepMCFilter_event_total =
0630         hepMCFilterEffStat_.numTotalPositiveEvents() + hepMCFilterEffStat_.numTotalNegativeEvents();
0631     double hepMCFilter_event_pass =
0632         hepMCFilterEffStat_.numPassPositiveEvents() + hepMCFilterEffStat_.numPassNegativeEvents();
0633     double hepMCFilter_event_eff = hepMCFilter_event_total > 0 ? hepMCFilter_event_pass / hepMCFilter_event_total : 0;
0634     double hepMCFilter_event_err =
0635         hepMCFilter_event_total > 0
0636             ? sqrt((1 - hepMCFilter_event_eff) * hepMCFilter_event_eff / hepMCFilter_event_total)
0637             : -1;
0638     edm::LogPrint("GenXSecAnalyzer") << "HepMC filter efficiency (event-level)= "
0639                                      << "(" << hepMCFilter_event_pass << ")"
0640                                      << " / "
0641                                      << "(" << hepMCFilter_event_total << ")"
0642                                      << " = " << std::scientific << std::setprecision(3) << hepMCFilter_event_eff
0643                                      << " +- " << hepMCFilter_event_err;
0644   }
0645
0646   // gen-particle filter efficiency
0647   if (filterOnlyEffStat_.sumWeights2() > 0) {
0648     double filterOnly_eff = filterOnlyEffStat_.filterEfficiency(-1);
0649     double filterOnly_err = filterOnlyEffStat_.filterEfficiencyError(-1);
0650
0651     edm::LogPrint("GenXSecAnalyzer") << "Filter efficiency (taking into account weights)= "
0652                                      << "(" << filterOnlyEffStat_.sumPassWeights() << ")"
0653                                      << " / "
0654                                      << "(" << filterOnlyEffStat_.sumWeights() << ")"
0655                                      << " = " << std::scientific << std::setprecision(3) << filterOnly_eff << " +- "
0656                                      << filterOnly_err;
0657
0658     double filterOnly_event_total =
0659         filterOnlyEffStat_.numTotalPositiveEvents() + filterOnlyEffStat_.numTotalNegativeEvents();
0660     double filterOnly_event_pass =
0661         filterOnlyEffStat_.numPassPositiveEvents() + filterOnlyEffStat_.numPassNegativeEvents();
0662     double filterOnly_event_eff = filterOnly_event_total > 0 ? filterOnly_event_pass / filterOnly_event_total : 0;
0663     double filterOnly_event_err = filterOnly_event_total > 0
0664                                       ? sqrt((1 - filterOnly_event_eff) * filterOnly_event_eff / filterOnly_event_total)
0665                                       : -1;
0666     edm::LogPrint("GenXSecAnalyzer") << "Filter efficiency (event-level)= "
0667                                      << "(" << filterOnly_event_pass << ")"
0668                                      << " / "
0669                                      << "(" << filterOnly_event_total << ")"
0670                                      << " = " << std::scientific << std::setprecision(3) << filterOnly_event_eff
0671                                      << " +- " << filterOnly_event_err << "    [TO BE USED IN MCM]";
0672
0673     // fill negative fraction of negative weights and uncertainty after filter
0674     final_fract_neg_w =
0675         filterOnly_event_pass > 0 ? filterOnlyEffStat_.numPassNegativeEvents() / (filterOnly_event_pass) : 0;
0676     final_fract_neg_w_unc =
0677         filterOnlyEffStat_.numPassNegativeEvents() > 0
0678             ? final_fract_neg_w * final_fract_neg_w / filterOnly_event_pass *
0679                   sqrt(filterOnlyEffStat_.numPassPositiveEvents() * filterOnlyEffStat_.numPassPositiveEvents() /
0680                            filterOnlyEffStat_.numPassNegativeEvents() +
0681                        filterOnlyEffStat_.numPassPositiveEvents())
0682             : 0;
0683   }
0684
0685   edm::LogPrint("GenXSecAnalyzer") << "\nAfter filter: final cross section = " << std::scientific
0686                                    << std::setprecision(3) << xsec_.value() << " +- " << xsec_.error() << " pb";
0687
0688   edm::LogPrint("GenXSecAnalyzer") << "After filter: final fraction of events with negative weights = "
0689                                    << std::scientific << std::setprecision(3) << final_fract_neg_w << " +- "
0690                                    << final_fract_neg_w_unc;
0691
0692   // L=[N*(1-2f)^2]/s
0693   double lumi_1M_evts =
0694       xsec_.value() > 0 ? 1e6 * (1 - 2 * final_fract_neg_w) * (1 - 2 * final_fract_neg_w) / xsec_.value() / 1e3 : 0;
0695   double lumi_1M_evts_unc =
0696       xsec_.value() > 0 ? (1 - 2 * final_fract_neg_w) * lumi_1M_evts *
0697                               sqrt(1e-6 + 16 * pow(final_fract_neg_w_unc, 2) / pow(1 - 2 * final_fract_neg_w, 2) +
0698                                    pow(xsec_.error() / xsec_.value(), 2))
0699                         : 0;
0700   edm::LogPrint("GenXSecAnalyzer") << "After filter: final equivalent lumi for 1M events (1/fb) = " << std::scientific
0701                                    << std::setprecision(3) << lumi_1M_evts << " +- " << lumi_1M_evts_unc;
0702 }
0703
0704 DEFINE_FWK_MODULE(GenXSecAnalyzer);