Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 13:29:38 

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   double sigSum = 0.0;
0389   double err2Sum = 0.0;
0390 
0391   std::vector<GenLumiInfoProduct::XSec> tempVector_before;
0392   std::vector<GenLumiInfoProduct::XSec> tempVector_after;
0393 
0394   // loop over different processes for each sample
0395   unsigned int vectorSize = iLumiInfo.getProcessInfos().size();
0396   for (unsigned int ip = 0; ip < vectorSize; ip++) {
0397     GenLumiInfoProduct::ProcessInfo proc = iLumiInfo.getProcessInfos()[ip];
0398     double hepxsec_value = proc.lheXSec().value();
0399     double hepxsec_error = proc.lheXSec().error() <= 0 ? 0 : proc.lheXSec().error();
0400     tempVector_before.push_back(GenLumiInfoProduct::XSec(hepxsec_value, hepxsec_error));
0401 
0402     sigSelSum += hepxsec_value;
0403     err2SelSum += hepxsec_error * hepxsec_error;
0404 
0405     // skips computation if jet matching efficiency=0
0406     if (proc.killed().n() < 1) {
0407       tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
0408       continue;
0409     }
0410 
0411     // computing jet matching efficiency for this process
0412     double fracAcc = 0;
0413     double ntotal = proc.nTotalPos() - proc.nTotalNeg();
0414     double npass = proc.nPassPos() - proc.nPassNeg();
0415     switch (hepidwtup_) {
0416       case 3:
0417       case -3:
0418         fracAcc = ntotal > 0 ? npass / ntotal : -1;
0419         break;
0420       default:
0421         fracAcc = proc.selected().sum() > 0 ? proc.killed().sum() / proc.selected().sum() : -1;
0422         break;
0423     }
0424 
0425     if (fracAcc <= 0) {
0426       tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
0427       continue;
0428     }
0429 
0430     // cross section after matching for this particular process
0431     double sigmaFin = hepxsec_value * fracAcc;
0432 
0433     // computing error on jet matching efficiency
0434     double relErr = 1.0;
0435     double efferr2 = 0;
0436     switch (hepidwtup_) {
0437       case 3:
0438       case -3: {
0439         double ntotal_pos = proc.nTotalPos();
0440         double effp = ntotal_pos > 0 ? (double)proc.nPassPos() / ntotal_pos : 0;
0441         double effp_err2 = ntotal_pos > 0 ? (1 - effp) * effp / ntotal_pos : 0;
0442 
0443         double ntotal_neg = proc.nTotalNeg();
0444         double effn = ntotal_neg > 0 ? (double)proc.nPassNeg() / ntotal_neg : 0;
0445         double effn_err2 = ntotal_neg > 0 ? (1 - effn) * effn / ntotal_neg : 0;
0446 
0447         efferr2 = ntotal > 0
0448                       ? (ntotal_pos * ntotal_pos * effp_err2 + ntotal_neg * ntotal_neg * effn_err2) / ntotal / ntotal
0449                       : 0;
0450         break;
0451       }
0452       default: {
0453         double denominator = pow(proc.selected().sum(), 4);
0454         double passw = proc.killed().sum();
0455         double passw2 = proc.killed().sum2();
0456         double failw = proc.selected().sum() - passw;
0457         double failw2 = proc.selected().sum2() - passw2;
0458         double numerator = (passw2 * failw * failw + failw2 * passw * passw);
0459 
0460         efferr2 = denominator > 0 ? numerator / denominator : 0;
0461         break;
0462       }
0463     }
0464     double delta2Veto = efferr2 / fracAcc / fracAcc;
0465 
0466     // computing total error on cross section after matching efficiency
0467 
0468     double sigma2Sum, sigma2Err;
0469     sigma2Sum = hepxsec_value * hepxsec_value;
0470     sigma2Err = hepxsec_error * hepxsec_error;
0471 
0472     double delta2Sum = delta2Veto + sigma2Err / sigma2Sum;
0473     relErr = (delta2Sum > 0.0 ? std::sqrt(delta2Sum) : 0.0);
0474     double deltaFin = sigmaFin * relErr;
0475 
0476     tempVector_after.push_back(GenLumiInfoProduct::XSec(sigmaFin, deltaFin));
0477 
0478     // sum of cross sections and errors over different processes
0479     sigSum += sigmaFin;
0480     err2Sum += deltaFin * deltaFin;
0481 
0482   }  // end of loop over different processes
0483   tempVector_before.push_back(GenLumiInfoProduct::XSec(sigSelSum, sqrt(err2SelSum)));
0484 
0485   double total_matcheff = jetMatchEffStat_[10000].filterEfficiency(hepidwtup_);
0486   double total_matcherr = jetMatchEffStat_[10000].filterEfficiencyError(hepidwtup_);
0487 
0488   double xsec_after = sigSelSum * total_matcheff;
0489   double xsecerr_after = (total_matcheff > 0 && sigSelSum > 0)
0490                              ? xsec_after * sqrt(err2SelSum / sigSelSum / sigSelSum +
0491                                                  total_matcherr * total_matcherr / total_matcheff / total_matcheff)
0492                              : 0;
0493 
0494   GenLumiInfoProduct::XSec result(xsec_after, xsecerr_after);
0495   tempVector_after.push_back(result);
0496 
0497   xsecBeforeMatching_ = tempVector_before;
0498   xsecAfterMatching_ = tempVector_after;
0499 
0500   return result;
0501 }
0502 
0503 void GenXSecAnalyzer::endJob() {
0504   edm::LogPrint("GenXSecAnalyzer") << "\n"
0505                                    << "------------------------------------"
0506                                    << "\n"
0507                                    << "GenXsecAnalyzer:"
0508                                    << "\n"
0509                                    << "------------------------------------";
0510 
0511   if (jetMatchEffStat_.empty()) {
0512     edm::LogPrint("GenXSecAnalyzer") << "------------------------------------"
0513                                      << "\n"
0514                                      << "Cross-section summary not available"
0515                                      << "\n"
0516                                      << "------------------------------------";
0517     return;
0518   }
0519 
0520   // fraction of negative weights
0521   double final_fract_neg_w = 0;
0522   double final_fract_neg_w_unc = 0;
0523 
0524   // below print out is only for combination of same physics MC samples and ME+Pythia MCs
0525 
0526   if (nMCs_ == 1 && hepidwtup_ != -1) {
0527     edm::LogPrint("GenXSecAnalyzer")
0528         << "-----------------------------------------------------------------------------------------------------------"
0529            "--------------------------------------------------------------- \n"
0530         << "Overall cross-section summary \n"
0531         << "-----------------------------------------------------------------------------------------------------------"
0532            "---------------------------------------------------------------";
0533     edm::LogPrint("GenXSecAnalyzer") << "Process\t\txsec_before [pb]\t\tpassed\tnposw\tnnegw\ttried\tnposw\tnnegw "
0534                                         "\txsec_match [pb]\t\t\taccepted [%]\t event_eff [%]";
0535 
0536     const unsigned sizeOfInfos = jetMatchEffStat_.size();
0537     const unsigned last = sizeOfInfos - 1;
0538     std::string *title = new std::string[sizeOfInfos];
0539     unsigned int i = 0;
0540     double jetmatch_eff = 0;
0541     double jetmatch_err = 0;
0542     double matching_eff = 1;
0543     double matching_efferr = 1;
0544 
0545     for (std::map<int, GenFilterInfo>::const_iterator iter = jetMatchEffStat_.begin(); iter != jetMatchEffStat_.end();
0546          ++iter, i++) {
0547       GenFilterInfo thisJetMatchStat = iter->second;
0548       GenFilterInfo thisEventEffStat =
0549           GenFilterInfo(thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
0550                         0,
0551                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents(),
0552                         0,
0553                         thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
0554                         thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
0555                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents(),
0556                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents());
0557 
0558       jetmatch_eff = thisJetMatchStat.filterEfficiency(hepidwtup_);
0559       jetmatch_err = thisJetMatchStat.filterEfficiencyError(hepidwtup_);
0560 
0561       if (i == last) {
0562         title[i] = "Total";
0563 
0564         edm::LogPrint("GenXSecAnalyzer")
0565             << "-------------------------------------------------------------------------------------------------------"
0566                "------------------------------------------------------------------- ";
0567 
0568         // fill negative fraction of negative weights and uncertainty after matching
0569         final_fract_neg_w = thisEventEffStat.numEventsPassed() > 0
0570                                 ? thisJetMatchStat.numPassNegativeEvents() / thisEventEffStat.numEventsPassed()
0571                                 : 0;
0572         final_fract_neg_w_unc =
0573             thisJetMatchStat.numPassNegativeEvents() > 0
0574                 ? final_fract_neg_w * final_fract_neg_w / thisEventEffStat.numEventsPassed() *
0575                       sqrt(thisJetMatchStat.numPassPositiveEvents() * thisJetMatchStat.numPassPositiveEvents() /
0576                                thisJetMatchStat.numPassNegativeEvents() +
0577                            thisJetMatchStat.numPassPositiveEvents())
0578                 : 0;
0579       } else {
0580         title[i] = Form("%d", i);
0581       }
0582 
0583       edm::LogPrint("GenXSecAnalyzer") << title[i] << "\t\t" << std::scientific << std::setprecision(3)
0584                                        << xsecBeforeMatching_[i].value() << " +/- " << xsecBeforeMatching_[i].error()
0585                                        << "\t\t" << thisEventEffStat.numEventsPassed() << "\t"
0586                                        << thisJetMatchStat.numPassPositiveEvents() << "\t"
0587                                        << thisJetMatchStat.numPassNegativeEvents() << "\t"
0588                                        << thisEventEffStat.numEventsTotal() << "\t"
0589                                        << thisJetMatchStat.numTotalPositiveEvents() << "\t"
0590                                        << thisJetMatchStat.numTotalNegativeEvents() << "\t" << std::scientific
0591                                        << std::setprecision(3) << xsecAfterMatching_[i].value() << " +/- "
0592                                        << xsecAfterMatching_[i].error() << "\t\t" << std::fixed << std::setprecision(1)
0593                                        << (jetmatch_eff * 100) << " +/- " << (jetmatch_err * 100) << "\t" << std::fixed
0594                                        << std::setprecision(1) << (thisEventEffStat.filterEfficiency(+3) * 100)
0595                                        << " +/- " << (thisEventEffStat.filterEfficiencyError(+3) * 100);
0596 
0597       matching_eff = thisEventEffStat.filterEfficiency(+3);
0598       matching_efferr = thisEventEffStat.filterEfficiencyError(+3);
0599     }
0600     delete[] title;
0601 
0602     edm::LogPrint("GenXSecAnalyzer")
0603         << "-----------------------------------------------------------------------------------------------------------"
0604            "---------------------------------------------------------------";
0605 
0606     edm::LogPrint("GenXSecAnalyzer") << "Before matching: total cross section = " << std::scientific
0607                                      << std::setprecision(3) << xsecBeforeMatching_[last].value() << " +- "
0608                                      << xsecBeforeMatching_[last].error() << " pb";
0609 
0610     edm::LogPrint("GenXSecAnalyzer") << "After matching: total cross section = " << std::scientific
0611                                      << std::setprecision(3) << xsecAfterMatching_[last].value() << " +- "
0612                                      << xsecAfterMatching_[last].error() << " pb";
0613 
0614     edm::LogPrint("GenXSecAnalyzer") << "Matching efficiency = " << std::fixed << std::setprecision(1) << matching_eff
0615                                      << " +/- " << matching_efferr << "   [TO BE USED IN MCM]";
0616 
0617   } else if (hepidwtup_ == -1)
0618     edm::LogPrint("GenXSecAnalyzer") << "Before Filter: total cross section = " << std::scientific
0619                                      << std::setprecision(3) << xsecPreFilter_.value() << " +- "
0620                                      << xsecPreFilter_.error() << " pb";
0621 
0622   // hepMC filter efficiency
0623   double hepMCFilter_eff = 1.0;
0624   double hepMCFilter_err = 0.0;
0625   if (hepMCFilterEffStat_.sumWeights2() > 0) {
0626     hepMCFilter_eff = hepMCFilterEffStat_.filterEfficiency(-1);
0627     hepMCFilter_err = hepMCFilterEffStat_.filterEfficiencyError(-1);
0628     edm::LogPrint("GenXSecAnalyzer") << "HepMC filter efficiency (taking into account weights)= "
0629                                      << "(" << hepMCFilterEffStat_.sumPassWeights() << ")"
0630                                      << " / "
0631                                      << "(" << hepMCFilterEffStat_.sumWeights() << ")"
0632                                      << " = " << std::scientific << std::setprecision(3) << hepMCFilter_eff << " +- "
0633                                      << hepMCFilter_err;
0634 
0635     double hepMCFilter_event_total =
0636         hepMCFilterEffStat_.numTotalPositiveEvents() + hepMCFilterEffStat_.numTotalNegativeEvents();
0637     double hepMCFilter_event_pass =
0638         hepMCFilterEffStat_.numPassPositiveEvents() + hepMCFilterEffStat_.numPassNegativeEvents();
0639     double hepMCFilter_event_eff = hepMCFilter_event_total > 0 ? hepMCFilter_event_pass / hepMCFilter_event_total : 0;
0640     double hepMCFilter_event_err =
0641         hepMCFilter_event_total > 0
0642             ? sqrt((1 - hepMCFilter_event_eff) * hepMCFilter_event_eff / hepMCFilter_event_total)
0643             : -1;
0644     edm::LogPrint("GenXSecAnalyzer") << "HepMC filter efficiency (event-level)= "
0645                                      << "(" << hepMCFilter_event_pass << ")"
0646                                      << " / "
0647                                      << "(" << hepMCFilter_event_total << ")"
0648                                      << " = " << std::scientific << std::setprecision(3) << hepMCFilter_event_eff
0649                                      << " +- " << hepMCFilter_event_err;
0650   }
0651 
0652   // gen-particle filter efficiency
0653   if (filterOnlyEffStat_.sumWeights2() > 0) {
0654     double filterOnly_eff = filterOnlyEffStat_.filterEfficiency(-1);
0655     double filterOnly_err = filterOnlyEffStat_.filterEfficiencyError(-1);
0656 
0657     edm::LogPrint("GenXSecAnalyzer") << "Filter efficiency (taking into account weights)= "
0658                                      << "(" << filterOnlyEffStat_.sumPassWeights() << ")"
0659                                      << " / "
0660                                      << "(" << filterOnlyEffStat_.sumWeights() << ")"
0661                                      << " = " << std::scientific << std::setprecision(3) << filterOnly_eff << " +- "
0662                                      << filterOnly_err;
0663 
0664     double filterOnly_event_total =
0665         filterOnlyEffStat_.numTotalPositiveEvents() + filterOnlyEffStat_.numTotalNegativeEvents();
0666     double filterOnly_event_pass =
0667         filterOnlyEffStat_.numPassPositiveEvents() + filterOnlyEffStat_.numPassNegativeEvents();
0668     double filterOnly_event_eff = filterOnly_event_total > 0 ? filterOnly_event_pass / filterOnly_event_total : 0;
0669     double filterOnly_event_err = filterOnly_event_total > 0
0670                                       ? sqrt((1 - filterOnly_event_eff) * filterOnly_event_eff / filterOnly_event_total)
0671                                       : -1;
0672     edm::LogPrint("GenXSecAnalyzer") << "Filter efficiency (event-level)= "
0673                                      << "(" << filterOnly_event_pass << ")"
0674                                      << " / "
0675                                      << "(" << filterOnly_event_total << ")"
0676                                      << " = " << std::scientific << std::setprecision(3) << filterOnly_event_eff
0677                                      << " +- " << filterOnly_event_err << "    [TO BE USED IN MCM]";
0678 
0679     // fill negative fraction of negative weights and uncertainty after filter
0680     final_fract_neg_w =
0681         filterOnly_event_pass > 0 ? filterOnlyEffStat_.numPassNegativeEvents() / (filterOnly_event_pass) : 0;
0682     final_fract_neg_w_unc =
0683         filterOnlyEffStat_.numPassNegativeEvents() > 0
0684             ? final_fract_neg_w * final_fract_neg_w / filterOnly_event_pass *
0685                   sqrt(filterOnlyEffStat_.numPassPositiveEvents() * filterOnlyEffStat_.numPassPositiveEvents() /
0686                            filterOnlyEffStat_.numPassNegativeEvents() +
0687                        filterOnlyEffStat_.numPassPositiveEvents())
0688             : 0;
0689   }
0690 
0691   edm::LogPrint("GenXSecAnalyzer") << "\nAfter filter: final cross section = " << std::scientific
0692                                    << std::setprecision(3) << xsec_.value() << " +- " << xsec_.error() << " pb";
0693 
0694   edm::LogPrint("GenXSecAnalyzer") << "After filter: final fraction of events with negative weights = "
0695                                    << std::scientific << std::setprecision(3) << final_fract_neg_w << " +- "
0696                                    << final_fract_neg_w_unc;
0697 
0698   // L=[N*(1-2f)^2]/s
0699   double lumi_1M_evts =
0700       xsec_.value() > 0 ? 1e6 * (1 - 2 * final_fract_neg_w) * (1 - 2 * final_fract_neg_w) / xsec_.value() / 1e3 : 0;
0701   double lumi_1M_evts_unc =
0702       xsec_.value() > 0 ? (1 - 2 * final_fract_neg_w) * lumi_1M_evts *
0703                               sqrt(1e-6 + 16 * pow(final_fract_neg_w_unc, 2) / pow(1 - 2 * final_fract_neg_w, 2) +
0704                                    pow(xsec_.error() / xsec_.value(), 2))
0705                         : 0;
0706   edm::LogPrint("GenXSecAnalyzer") << "After filter: final equivalent lumi for 1M events (1/fb) = " << std::scientific
0707                                    << std::setprecision(3) << lumi_1M_evts << " +- " << lumi_1M_evts_unc;
0708 }
0709 
0710 DEFINE_FWK_MODULE(GenXSecAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!