Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:26 

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   if (nullptr == runC)
0298     return;
0299 
0300   std::lock_guard l{mutex_};
0301 
0302   // compute cross section for this run first
0303   // set the correct combined LHE+filter cross sections
0304   unsigned int i = 0;
0305   std::vector<GenLumiInfoProduct::ProcessInfo> newInfos;
0306   for (std::map<int, GenLumiInfoProduct::XSec>::const_iterator iter = runC->currentLumiBlockLHEXSec_.begin();
0307        iter != runC->currentLumiBlockLHEXSec_.end();
0308        ++iter, i++) {
0309     GenLumiInfoProduct::ProcessInfo temp = runC->product_.getProcessInfos()[i];
0310     temp.setLheXSec(iter->second.value(), iter->second.error());
0311     newInfos.push_back(temp);
0312   }
0313   runC->product_.setProcessInfo(newInfos);
0314 
0315   const GenLumiInfoProduct::XSec thisRunXSecPre = compute(runC->product_);
0316   // xsection after matching before filters
0317   combine(xsecPreFilter_, totalWeightPre_, thisRunXSecPre, runC->thisRunWeightPre_);
0318 
0319   double thisHepFilterEff = 1;
0320   double thisHepFilterErr = 0;
0321 
0322   if (runC->hepMCFilterEffRun_.sumWeights2() > 0) {
0323     thisHepFilterEff = runC->hepMCFilterEffRun_.filterEfficiency(hepidwtup_);
0324     thisHepFilterErr = runC->hepMCFilterEffRun_.filterEfficiencyError(hepidwtup_);
0325     if (thisHepFilterEff < 0) {
0326       thisHepFilterEff = 1;
0327       thisHepFilterErr = 0;
0328     }
0329   }
0330 
0331   double thisGenFilterEff = 1;
0332   double thisGenFilterErr = 0;
0333 
0334   if (runC->filterOnlyEffRun_.sumWeights2() > 0) {
0335     thisGenFilterEff = runC->filterOnlyEffRun_.filterEfficiency(hepidwtup_);
0336     thisGenFilterErr = runC->filterOnlyEffRun_.filterEfficiencyError(hepidwtup_);
0337     if (thisGenFilterEff < 0) {
0338       thisGenFilterEff = 1;
0339       thisGenFilterErr = 0;
0340     }
0341   }
0342   double thisXsec = thisRunXSecPre.value() > 0 ? thisHepFilterEff * thisGenFilterEff * thisRunXSecPre.value() : 0;
0343   double thisErr =
0344       thisRunXSecPre.value() > 0
0345           ? thisXsec * sqrt(pow(TMath::Max(thisRunXSecPre.error(), (double)0) / thisRunXSecPre.value(), 2) +
0346                             pow(thisHepFilterErr / thisHepFilterEff, 2) + pow(thisGenFilterErr / thisGenFilterEff, 2))
0347           : 0;
0348   const GenLumiInfoProduct::XSec thisRunXSec = GenLumiInfoProduct::XSec(thisXsec, thisErr);
0349   combine(xsec_, totalWeight_, thisRunXSec, runC->thisRunWeight_);
0350 }
0351 
0352 void GenXSecAnalyzer::combine(double &finalValue,
0353                               double &finalError,
0354                               double &finalWeight,
0355                               const double &currentValue,
0356                               const double &currentError,
0357                               const double &currentWeight) const {
0358   if (finalValue <= 0) {
0359     finalValue = currentValue;
0360     finalError = currentError;
0361     finalWeight += currentWeight;
0362   } else {
0363     double wgt1 = (finalError <= 0 || currentError <= 0) ? finalWeight : 1 / (finalError * finalError);
0364     double wgt2 = (finalError <= 0 || currentError <= 0) ? currentWeight : 1 / (currentError * currentError);
0365     double xsec = (wgt1 * finalValue + wgt2 * currentValue) / (wgt1 + wgt2);
0366     double err = (finalError <= 0 || currentError <= 0) ? 0 : 1.0 / std::sqrt(wgt1 + wgt2);
0367     finalValue = xsec;
0368     finalError = err;
0369     finalWeight += currentWeight;
0370   }
0371   return;
0372 }
0373 
0374 void GenXSecAnalyzer::combine(GenLumiInfoProduct::XSec &finalXSec,
0375                               double &totalw,
0376                               const GenLumiInfoProduct::XSec &thisRunXSec,
0377                               const double &thisw) const {
0378   double value = finalXSec.value();
0379   double error = finalXSec.error();
0380   double thisValue = thisRunXSec.value();
0381   double thisError = thisRunXSec.error();
0382   combine(value, error, totalw, thisValue, thisError, thisw);
0383   finalXSec = GenLumiInfoProduct::XSec(value, error);
0384   return;
0385 }
0386 
0387 GenLumiInfoProduct::XSec GenXSecAnalyzer::compute(const GenLumiInfoProduct &iLumiInfo) const {
0388   // sum of cross sections and errors over different processes
0389   double sigSelSum = 0.0;
0390   double err2SelSum = 0.0;
0391 
0392   std::vector<GenLumiInfoProduct::XSec> tempVector_before;
0393   std::vector<GenLumiInfoProduct::XSec> tempVector_after;
0394 
0395   // loop over different processes for each sample
0396   unsigned int vectorSize = iLumiInfo.getProcessInfos().size();
0397   for (unsigned int ip = 0; ip < vectorSize; ip++) {
0398     GenLumiInfoProduct::ProcessInfo proc = iLumiInfo.getProcessInfos()[ip];
0399     double hepxsec_value = proc.lheXSec().value();
0400     double hepxsec_error = proc.lheXSec().error() <= 0 ? 0 : proc.lheXSec().error();
0401     tempVector_before.push_back(GenLumiInfoProduct::XSec(hepxsec_value, hepxsec_error));
0402 
0403     sigSelSum += hepxsec_value;
0404     err2SelSum += hepxsec_error * hepxsec_error;
0405 
0406     // skips computation if jet matching efficiency=0
0407     if (proc.killed().n() < 1) {
0408       tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
0409       continue;
0410     }
0411 
0412     // computing jet matching efficiency for this process
0413     double fracAcc = 0;
0414     double ntotal = proc.nTotalPos() - proc.nTotalNeg();
0415     double npass = proc.nPassPos() - proc.nPassNeg();
0416     switch (hepidwtup_) {
0417       case 3:
0418       case -3:
0419         fracAcc = ntotal > 0 ? npass / ntotal : -1;
0420         break;
0421       default:
0422         fracAcc = proc.selected().sum() > 0 ? proc.killed().sum() / proc.selected().sum() : -1;
0423         break;
0424     }
0425 
0426     if (fracAcc <= 0) {
0427       tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
0428       continue;
0429     }
0430 
0431     // cross section after matching for this particular process
0432     double sigmaFin = hepxsec_value * fracAcc;
0433 
0434     // computing error on jet matching efficiency
0435     double relErr = 1.0;
0436     double efferr2 = 0;
0437     switch (hepidwtup_) {
0438       case 3:
0439       case -3: {
0440         double ntotal_pos = proc.nTotalPos();
0441         double effp = ntotal_pos > 0 ? (double)proc.nPassPos() / ntotal_pos : 0;
0442         double effp_err2 = ntotal_pos > 0 ? (1 - effp) * effp / ntotal_pos : 0;
0443 
0444         double ntotal_neg = proc.nTotalNeg();
0445         double effn = ntotal_neg > 0 ? (double)proc.nPassNeg() / ntotal_neg : 0;
0446         double effn_err2 = ntotal_neg > 0 ? (1 - effn) * effn / ntotal_neg : 0;
0447 
0448         efferr2 = ntotal > 0
0449                       ? (ntotal_pos * ntotal_pos * effp_err2 + ntotal_neg * ntotal_neg * effn_err2) / ntotal / ntotal
0450                       : 0;
0451         break;
0452       }
0453       default: {
0454         double denominator = pow(proc.selected().sum(), 4);
0455         double passw = proc.killed().sum();
0456         double passw2 = proc.killed().sum2();
0457         double failw = proc.selected().sum() - passw;
0458         double failw2 = proc.selected().sum2() - passw2;
0459         double numerator = (passw2 * failw * failw + failw2 * passw * passw);
0460 
0461         efferr2 = denominator > 0 ? numerator / denominator : 0;
0462         break;
0463       }
0464     }
0465     double delta2Veto = efferr2 / fracAcc / fracAcc;
0466 
0467     // computing total error on cross section after matching efficiency
0468 
0469     double sigma2Sum, sigma2Err;
0470     sigma2Sum = hepxsec_value * hepxsec_value;
0471     sigma2Err = hepxsec_error * hepxsec_error;
0472 
0473     double delta2Sum = delta2Veto + sigma2Err / sigma2Sum;
0474     relErr = (delta2Sum > 0.0 ? std::sqrt(delta2Sum) : 0.0);
0475     double deltaFin = sigmaFin * relErr;
0476 
0477     tempVector_after.push_back(GenLumiInfoProduct::XSec(sigmaFin, deltaFin));
0478 
0479   }  // end of loop over different processes
0480   tempVector_before.push_back(GenLumiInfoProduct::XSec(sigSelSum, sqrt(err2SelSum)));
0481 
0482   double total_matcheff = jetMatchEffStat_[10000].filterEfficiency(hepidwtup_);
0483   double total_matcherr = jetMatchEffStat_[10000].filterEfficiencyError(hepidwtup_);
0484 
0485   double xsec_after = sigSelSum * total_matcheff;
0486   double xsecerr_after = (total_matcheff > 0 && sigSelSum > 0)
0487                              ? xsec_after * sqrt(err2SelSum / sigSelSum / sigSelSum +
0488                                                  total_matcherr * total_matcherr / total_matcheff / total_matcheff)
0489                              : 0;
0490 
0491   GenLumiInfoProduct::XSec result(xsec_after, xsecerr_after);
0492   tempVector_after.push_back(result);
0493 
0494   xsecBeforeMatching_ = tempVector_before;
0495   xsecAfterMatching_ = tempVector_after;
0496 
0497   return result;
0498 }
0499 
0500 void GenXSecAnalyzer::endJob() {
0501   edm::LogPrint("GenXSecAnalyzer") << "\n"
0502                                    << "------------------------------------"
0503                                    << "\n"
0504                                    << "GenXsecAnalyzer:"
0505                                    << "\n"
0506                                    << "------------------------------------";
0507 
0508   if (jetMatchEffStat_.empty()) {
0509     edm::LogPrint("GenXSecAnalyzer") << "------------------------------------"
0510                                      << "\n"
0511                                      << "Cross-section summary not available"
0512                                      << "\n"
0513                                      << "------------------------------------";
0514     return;
0515   }
0516 
0517   // fraction of negative weights
0518   double final_fract_neg_w = 0;
0519   double final_fract_neg_w_unc = 0;
0520 
0521   // below print out is only for combination of same physics MC samples and ME+Pythia MCs
0522 
0523   if (nMCs_ == 1 && hepidwtup_ != -1) {
0524     edm::LogPrint("GenXSecAnalyzer")
0525         << "-----------------------------------------------------------------------------------------------------------"
0526            "--------------------------------------------------------------- \n"
0527         << "Overall cross-section summary \n"
0528         << "-----------------------------------------------------------------------------------------------------------"
0529            "---------------------------------------------------------------";
0530     edm::LogPrint("GenXSecAnalyzer") << "Process\t\txsec_before [pb]\t\tpassed\tnposw\tnnegw\ttried\tnposw\tnnegw "
0531                                         "\txsec_match [pb]\t\t\taccepted [%]\t event_eff [%]";
0532 
0533     const unsigned sizeOfInfos = jetMatchEffStat_.size();
0534     const unsigned last = sizeOfInfos - 1;
0535     std::string *title = new std::string[sizeOfInfos];
0536     unsigned int i = 0;
0537     double jetmatch_eff = 0;
0538     double jetmatch_err = 0;
0539     double matching_eff = 1;
0540     double matching_efferr = 1;
0541 
0542     for (std::map<int, GenFilterInfo>::const_iterator iter = jetMatchEffStat_.begin(); iter != jetMatchEffStat_.end();
0543          ++iter, i++) {
0544       GenFilterInfo thisJetMatchStat = iter->second;
0545       GenFilterInfo thisEventEffStat =
0546           GenFilterInfo(thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
0547                         0,
0548                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents(),
0549                         0,
0550                         thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
0551                         thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
0552                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents(),
0553                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents());
0554 
0555       jetmatch_eff = thisJetMatchStat.filterEfficiency(hepidwtup_);
0556       jetmatch_err = thisJetMatchStat.filterEfficiencyError(hepidwtup_);
0557 
0558       if (i == last) {
0559         title[i] = "Total";
0560 
0561         edm::LogPrint("GenXSecAnalyzer")
0562             << "-------------------------------------------------------------------------------------------------------"
0563                "------------------------------------------------------------------- ";
0564 
0565         // fill negative fraction of negative weights and uncertainty after matching
0566         final_fract_neg_w = thisEventEffStat.numEventsPassed() > 0
0567                                 ? thisJetMatchStat.numPassNegativeEvents() / thisEventEffStat.numEventsPassed()
0568                                 : 0;
0569         final_fract_neg_w_unc =
0570             thisJetMatchStat.numPassNegativeEvents() > 0
0571                 ? final_fract_neg_w * final_fract_neg_w / thisEventEffStat.numEventsPassed() *
0572                       sqrt(thisJetMatchStat.numPassPositiveEvents() * thisJetMatchStat.numPassPositiveEvents() /
0573                                thisJetMatchStat.numPassNegativeEvents() +
0574                            thisJetMatchStat.numPassPositiveEvents())
0575                 : 0;
0576       } else {
0577         title[i] = Form("%d", i);
0578       }
0579 
0580       edm::LogPrint("GenXSecAnalyzer") << title[i] << "\t\t" << std::scientific << std::setprecision(3)
0581                                        << xsecBeforeMatching_[i].value() << " +/- " << xsecBeforeMatching_[i].error()
0582                                        << "\t\t" << thisEventEffStat.numEventsPassed() << "\t"
0583                                        << thisJetMatchStat.numPassPositiveEvents() << "\t"
0584                                        << thisJetMatchStat.numPassNegativeEvents() << "\t"
0585                                        << thisEventEffStat.numEventsTotal() << "\t"
0586                                        << thisJetMatchStat.numTotalPositiveEvents() << "\t"
0587                                        << thisJetMatchStat.numTotalNegativeEvents() << "\t" << std::scientific
0588                                        << std::setprecision(3) << xsecAfterMatching_[i].value() << " +/- "
0589                                        << xsecAfterMatching_[i].error() << "\t\t" << std::fixed << std::setprecision(1)
0590                                        << (jetmatch_eff * 100) << " +/- " << (jetmatch_err * 100) << "\t" << std::fixed
0591                                        << std::setprecision(1) << (thisEventEffStat.filterEfficiency(+3) * 100)
0592                                        << " +/- " << (thisEventEffStat.filterEfficiencyError(+3) * 100);
0593 
0594       matching_eff = thisEventEffStat.filterEfficiency(+3);
0595       matching_efferr = thisEventEffStat.filterEfficiencyError(+3);
0596     }
0597     delete[] title;
0598 
0599     edm::LogPrint("GenXSecAnalyzer")
0600         << "-----------------------------------------------------------------------------------------------------------"
0601            "---------------------------------------------------------------";
0602 
0603     edm::LogPrint("GenXSecAnalyzer") << "Before matching: total cross section = " << std::scientific
0604                                      << std::setprecision(3) << xsecBeforeMatching_[last].value() << " +- "
0605                                      << xsecBeforeMatching_[last].error() << " pb";
0606 
0607     edm::LogPrint("GenXSecAnalyzer") << "After matching: total cross section = " << std::scientific
0608                                      << std::setprecision(3) << xsecAfterMatching_[last].value() << " +- "
0609                                      << xsecAfterMatching_[last].error() << " pb";
0610 
0611     edm::LogPrint("GenXSecAnalyzer") << "Matching efficiency = " << std::fixed << std::setprecision(1) << matching_eff
0612                                      << " +/- " << matching_efferr << "   [TO BE USED IN MCM]";
0613 
0614   } else if (hepidwtup_ == -1)
0615     edm::LogPrint("GenXSecAnalyzer") << "Before Filter: total cross section = " << std::scientific
0616                                      << std::setprecision(3) << xsecPreFilter_.value() << " +- "
0617                                      << xsecPreFilter_.error() << " pb";
0618 
0619   // hepMC filter efficiency
0620   double hepMCFilter_eff = 1.0;
0621   double hepMCFilter_err = 0.0;
0622   if (hepMCFilterEffStat_.sumWeights2() > 0) {
0623     hepMCFilter_eff = hepMCFilterEffStat_.filterEfficiency(-1);
0624     hepMCFilter_err = hepMCFilterEffStat_.filterEfficiencyError(-1);
0625     edm::LogPrint("GenXSecAnalyzer") << "HepMC filter efficiency (taking into account weights)= "
0626                                      << "(" << hepMCFilterEffStat_.sumPassWeights() << ")"
0627                                      << " / "
0628                                      << "(" << hepMCFilterEffStat_.sumWeights() << ")"
0629                                      << " = " << std::scientific << std::setprecision(3) << hepMCFilter_eff << " +- "
0630                                      << hepMCFilter_err;
0631 
0632     double hepMCFilter_event_total =
0633         hepMCFilterEffStat_.numTotalPositiveEvents() + hepMCFilterEffStat_.numTotalNegativeEvents();
0634     double hepMCFilter_event_pass =
0635         hepMCFilterEffStat_.numPassPositiveEvents() + hepMCFilterEffStat_.numPassNegativeEvents();
0636     double hepMCFilter_event_eff = hepMCFilter_event_total > 0 ? hepMCFilter_event_pass / hepMCFilter_event_total : 0;
0637     double hepMCFilter_event_err =
0638         hepMCFilter_event_total > 0
0639             ? sqrt((1 - hepMCFilter_event_eff) * hepMCFilter_event_eff / hepMCFilter_event_total)
0640             : -1;
0641     edm::LogPrint("GenXSecAnalyzer") << "HepMC filter efficiency (event-level)= "
0642                                      << "(" << hepMCFilter_event_pass << ")"
0643                                      << " / "
0644                                      << "(" << hepMCFilter_event_total << ")"
0645                                      << " = " << std::scientific << std::setprecision(3) << hepMCFilter_event_eff
0646                                      << " +- " << hepMCFilter_event_err;
0647   }
0648 
0649   // gen-particle filter efficiency
0650   if (filterOnlyEffStat_.sumWeights2() > 0) {
0651     double filterOnly_eff = filterOnlyEffStat_.filterEfficiency(-1);
0652     double filterOnly_err = filterOnlyEffStat_.filterEfficiencyError(-1);
0653 
0654     edm::LogPrint("GenXSecAnalyzer") << "Filter efficiency (taking into account weights)= "
0655                                      << "(" << filterOnlyEffStat_.sumPassWeights() << ")"
0656                                      << " / "
0657                                      << "(" << filterOnlyEffStat_.sumWeights() << ")"
0658                                      << " = " << std::scientific << std::setprecision(3) << filterOnly_eff << " +- "
0659                                      << filterOnly_err;
0660 
0661     double filterOnly_event_total =
0662         filterOnlyEffStat_.numTotalPositiveEvents() + filterOnlyEffStat_.numTotalNegativeEvents();
0663     double filterOnly_event_pass =
0664         filterOnlyEffStat_.numPassPositiveEvents() + filterOnlyEffStat_.numPassNegativeEvents();
0665     double filterOnly_event_eff = filterOnly_event_total > 0 ? filterOnly_event_pass / filterOnly_event_total : 0;
0666     double filterOnly_event_err = filterOnly_event_total > 0
0667                                       ? sqrt((1 - filterOnly_event_eff) * filterOnly_event_eff / filterOnly_event_total)
0668                                       : -1;
0669     edm::LogPrint("GenXSecAnalyzer") << "Filter efficiency (event-level)= "
0670                                      << "(" << filterOnly_event_pass << ")"
0671                                      << " / "
0672                                      << "(" << filterOnly_event_total << ")"
0673                                      << " = " << std::scientific << std::setprecision(3) << filterOnly_event_eff
0674                                      << " +- " << filterOnly_event_err << "    [TO BE USED IN MCM]";
0675 
0676     // fill negative fraction of negative weights and uncertainty after filter
0677     final_fract_neg_w =
0678         filterOnly_event_pass > 0 ? filterOnlyEffStat_.numPassNegativeEvents() / (filterOnly_event_pass) : 0;
0679     final_fract_neg_w_unc =
0680         filterOnlyEffStat_.numPassNegativeEvents() > 0
0681             ? final_fract_neg_w * final_fract_neg_w / filterOnly_event_pass *
0682                   sqrt(filterOnlyEffStat_.numPassPositiveEvents() * filterOnlyEffStat_.numPassPositiveEvents() /
0683                            filterOnlyEffStat_.numPassNegativeEvents() +
0684                        filterOnlyEffStat_.numPassPositiveEvents())
0685             : 0;
0686   }
0687 
0688   edm::LogPrint("GenXSecAnalyzer") << "\nAfter filter: final cross section = " << std::scientific
0689                                    << std::setprecision(3) << xsec_.value() << " +- " << xsec_.error() << " pb";
0690 
0691   edm::LogPrint("GenXSecAnalyzer") << "After filter: final fraction of events with negative weights = "
0692                                    << std::scientific << std::setprecision(3) << final_fract_neg_w << " +- "
0693                                    << final_fract_neg_w_unc;
0694 
0695   // L=[N*(1-2f)^2]/s
0696   double lumi_1M_evts =
0697       xsec_.value() > 0 ? 1e6 * (1 - 2 * final_fract_neg_w) * (1 - 2 * final_fract_neg_w) / xsec_.value() / 1e3 : 0;
0698   double lumi_1M_evts_unc =
0699       xsec_.value() > 0 ? (1 - 2 * final_fract_neg_w) * lumi_1M_evts *
0700                               sqrt(1e-6 + 16 * pow(final_fract_neg_w_unc, 2) / pow(1 - 2 * final_fract_neg_w, 2) +
0701                                    pow(xsec_.error() / xsec_.value(), 2))
0702                         : 0;
0703   edm::LogPrint("GenXSecAnalyzer") << "After filter: final equivalent lumi for 1M events (1/fb) = " << std::scientific
0704                                    << std::setprecision(3) << lumi_1M_evts << " +- " << lumi_1M_evts_unc;
0705 }
0706 
0707 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!