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
0007
0008
0009
0010 #include <memory>
0011 #include <vector>
0012 #include <map>
0013
0014
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
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
0042 CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable double thisRunWeightPre_ = 0;
0043
0044
0045 CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable double thisRunWeight_ = 0;
0046
0047
0048
0049 CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable GenLumiInfoProduct product_;
0050
0051
0052
0053 CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable GenFilterInfo filterOnlyEffRun_;
0054
0055
0056 CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable GenFilterInfo hepMCFilterEffRun_;
0057
0058
0059
0060
0061 CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable std::map<int, GenLumiInfoProduct::XSec> previousLumiBlockLHEXSec_;
0062
0063
0064
0065
0066 CMS_THREAD_GUARD(GenXSecAnalyzer::mutex_) mutable std::map<int, GenLumiInfoProduct::XSec> currentLumiBlockLHEXSec_;
0067 };
0068 }
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
0086 GenLumiInfoProduct::XSec compute(const GenLumiInfoProduct &) const;
0087
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
0097
0098 mutable std::atomic<int> nMCs_;
0099
0100 mutable std::atomic<int> hepidwtup_;
0101
0102 mutable std::mutex mutex_;
0103
0104
0105 CMS_THREAD_GUARD(mutex_) mutable double totalWeightPre_;
0106
0107
0108 CMS_THREAD_GUARD(mutex_) mutable double totalWeight_;
0109
0110
0111 CMS_THREAD_GUARD(mutex_) mutable GenLumiInfoProduct::XSec xsecPreFilter_;
0112
0113
0114 CMS_THREAD_GUARD(mutex_) mutable GenLumiInfoProduct::XSec xsec_;
0115
0116
0117 CMS_THREAD_GUARD(mutex_) mutable GenFilterInfo filterOnlyEffStat_;
0118
0119
0120 CMS_THREAD_GUARD(mutex_) mutable GenFilterInfo hepMCFilterEffStat_;
0121
0122
0123
0124
0125
0126 CMS_THREAD_GUARD(mutex_) mutable std::vector<GenLumiInfoProduct::XSec> xsecBeforeMatching_;
0127
0128 CMS_THREAD_GUARD(mutex_) mutable std::vector<GenLumiInfoProduct::XSec> xsecAfterMatching_;
0129
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
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
0185
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
0227
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
0248 jetMatchEffStat_[10000].mergeProduct(tempInfo);
0249 double currentValue = theProcesses[ip].lheXSec().value();
0250 double currentError = theProcesses[ip].lheXSec().error();
0251
0252
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)
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
0265 thisRunWeightPre += totalw;
0266
0267 }
0268
0269 else {
0270 x = tempInfo;
0271 y = theProcesses[ip].lheXSec();
0272 thisRunWeightPre += totalw;
0273 }
0274
0275 runC->previousLumiBlockLHEXSec_[id] = theProcesses[ip].lheXSec();
0276 }
0277
0278 return;
0279 }
0280
0281 void GenXSecAnalyzer::globalEndRun(edm::Run const &iRun, edm::EventSetup const &) const {
0282
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
0300
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
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 ¤tValue,
0353 const double ¤tError,
0354 const double ¤tWeight) 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
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
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
0406 if (proc.killed().n() < 1) {
0407 tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
0408 continue;
0409 }
0410
0411
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
0431 double sigmaFin = hepxsec_value * fracAcc;
0432
0433
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
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
0479 sigSum += sigmaFin;
0480 err2Sum += deltaFin * deltaFin;
0481
0482 }
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
0521 double final_fract_neg_w = 0;
0522 double final_fract_neg_w_unc = 0;
0523
0524
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
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
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
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
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
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);