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
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 if (nullptr == runC)
0298 return;
0299
0300 std::lock_guard l{mutex_};
0301
0302
0303
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
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 ¤tValue,
0356 const double ¤tError,
0357 const double ¤tWeight) 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
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
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
0407 if (proc.killed().n() < 1) {
0408 tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
0409 continue;
0410 }
0411
0412
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
0432 double sigmaFin = hepxsec_value * fracAcc;
0433
0434
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
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 }
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
0518 double final_fract_neg_w = 0;
0519 double final_fract_neg_w_unc = 0;
0520
0521
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
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
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
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
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
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);