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