#include <GenXSecAnalyzer.h>

Inheritance diagram for GenXSecAnalyzer:

Public Member Functions
const double	final_xsec_error () const

const double	final_xsec_value () const

	GenXSecAnalyzer (const edm::ParameterSet &)

	~GenXSecAnalyzer () override

Public Member Functions inherited from edm::one::EDAnalyzer< edm::one::WatchRuns, edm::one::WatchLuminosityBlocks >
	EDAnalyzer ()=default

SerialTaskQueue *	globalLuminosityBlocksQueue () final

SerialTaskQueue *	globalRunsQueue () final

bool	wantsGlobalLuminosityBlocks () const final

bool	wantsGlobalRuns () const final

Public Member Functions inherited from edm::one::EDAnalyzerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

	EDAnalyzerBase ()

ModuleDescription const &	moduleDescription () const

bool	wantsStreamLuminosityBlocks () const

bool	wantsStreamRuns () const

	~EDAnalyzerBase () override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase const &)=delete

	EDConsumerBase (EDConsumerBase &&)=default

ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

EDConsumerBase &	operator= (EDConsumerBase &&)=default

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
void	analyze (const edm::Event &, const edm::EventSetup &) override

void	beginJob () override

void	beginLuminosityBlock (edm::LuminosityBlock const &, edm::EventSetup const &) override

void	beginRun (edm::Run const &, edm::EventSetup const &) override

void	combine (GenLumiInfoProduct::XSec &, double &, const GenLumiInfoProduct::XSec &, const double &)

void	combine (double &, double &, double &, const double &, const double &, const double &)

GenLumiInfoProduct::XSec	compute (const GenLumiInfoProduct &)

void	endJob () override

void	endLuminosityBlock (edm::LuminosityBlock const &, edm::EventSetup const &) override

void	endRun (edm::Run const &, edm::EventSetup const &) override

Private Attributes
std::map< int, GenLumiInfoProduct::XSec >	currentLumiBlockLHEXSec_

GenFilterInfo	filterOnlyEffRun_

GenFilterInfo	filterOnlyEffStat_

edm::EDGetTokenT< GenFilterInfo >	genFilterInfoToken_

edm::EDGetTokenT< GenLumiInfoProduct >	genLumiInfoToken_

int	hepidwtup_

GenFilterInfo	hepMCFilterEffRun_

GenFilterInfo	hepMCFilterEffStat_

edm::EDGetTokenT< GenFilterInfo >	hepMCFilterInfoToken_

std::map< int, GenFilterInfo >	jetMatchEffStat_

edm::EDGetTokenT< LHERunInfoProduct >	lheRunInfoToken_

int	nMCs_

std::map< int, GenLumiInfoProduct::XSec >	previousLumiBlockLHEXSec_

GenLumiInfoProduct	product_

double	thisRunWeight_

double	thisRunWeightPre_

double	totalWeight_

double	totalWeightPre_

GenLumiInfoProduct::XSec	xsec_

std::vector< GenLumiInfoProduct::XSec >	xsecAfterMatching_

std::vector< GenLumiInfoProduct::XSec >	xsecBeforeMatching_

GenLumiInfoProduct::XSec	xsecPreFilter_

Additional Inherited Members
Public Types inherited from edm::one::EDAnalyzerBase
typedef EDAnalyzerBase	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Static Public Member Functions inherited from edm::one::EDAnalyzerBase
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

Detailed Description

Definition at line 33 of file GenXSecAnalyzer.h.

Constructor & Destructor Documentation

GenXSecAnalyzer::GenXSecAnalyzer ( const edm::ParameterSet & iConfig )

explicit

Definition at line 7 of file GenXSecAnalyzer.cc.

References currentLumiBlockLHEXSec_, genFilterInfoToken_, genLumiInfoToken_, hepMCFilterInfoToken_, HLT_2018_cff::InputTag, jetMatchEffStat_, lheRunInfoToken_, previousLumiBlockLHEXSec_, xsecAfterMatching_, and xsecBeforeMatching_.

     : nMCs_(0),
       hepidwtup_(-9999),
       totalWeightPre_(0),
       thisRunWeightPre_(0),
       totalWeight_(0),
       thisRunWeight_(0),
       xsecPreFilter_(-1, -1),
       xsec_(-1, -1),
       product_(GenLumiInfoProduct(-9999)),
       filterOnlyEffRun_(0, 0, 0, 0, 0., 0., 0., 0.),
       hepMCFilterEffRun_(0, 0, 0, 0, 0., 0., 0., 0.),
       filterOnlyEffStat_(0, 0, 0, 0, 0., 0., 0., 0.),
       hepMCFilterEffStat_(0, 0, 0, 0, 0., 0., 0., 0.) {
   xsecBeforeMatching_.clear();
   xsecAfterMatching_.clear();
   jetMatchEffStat_.clear();
   previousLumiBlockLHEXSec_.clear();
   currentLumiBlockLHEXSec_.clear();
 
   genFilterInfoToken_ = consumes<GenFilterInfo, edm::InLumi>(edm::InputTag("genFilterEfficiencyProducer", ""));
   hepMCFilterInfoToken_ = consumes<GenFilterInfo, edm::InLumi>(edm::InputTag("generator", ""));
   genLumiInfoToken_ = consumes<GenLumiInfoProduct, edm::InLumi>(edm::InputTag("generator", ""));
   lheRunInfoToken_ = consumes<LHERunInfoProduct, edm::InRun>(edm::InputTag("externalLHEProducer", ""));
 }

GenXSecAnalyzer::~GenXSecAnalyzer ( )

override

Definition at line 33 of file GenXSecAnalyzer.cc.

33 {}

Member Function Documentation

void GenXSecAnalyzer::analyze	(	const edm::Event &	,
		const edm::EventSetup &
	)

overrideprivate

Definition at line 67 of file GenXSecAnalyzer.cc.

Referenced by final_xsec_error().

67 {}

void GenXSecAnalyzer::beginJob ( void )

overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 35 of file GenXSecAnalyzer.cc.

References currentLumiBlockLHEXSec_, jetMatchEffStat_, previousLumiBlockLHEXSec_, xsecAfterMatching_, and xsecBeforeMatching_.

Referenced by final_xsec_error().

                                {
   xsecBeforeMatching_.clear();
   xsecAfterMatching_.clear();
   jetMatchEffStat_.clear();
   previousLumiBlockLHEXSec_.clear();
   currentLumiBlockLHEXSec_.clear();
 }

void GenXSecAnalyzer::beginLuminosityBlock	(	edm::LuminosityBlock const &	iLumi,
		edm::EventSetup const &
	)

overrideprivate

Definition at line 65 of file GenXSecAnalyzer.cc.

Referenced by final_xsec_error().

65 {}

void GenXSecAnalyzer::beginRun	(	edm::Run const &	iRun,
		edm::EventSetup const &
	)

overrideprivate

Definition at line 43 of file GenXSecAnalyzer.cc.

References currentLumiBlockLHEXSec_, filterOnlyEffRun_, hepMCFilterEffRun_, jetMatchEffStat_, nMCs_, previousLumiBlockLHEXSec_, product_, thisRunWeight_, thisRunWeightPre_, xsecAfterMatching_, and xsecBeforeMatching_.

Referenced by final_xsec_error().

                                                                        {
   // initialization for every different physics MC
 
   nMCs_++;
 
   thisRunWeightPre_ = 0;
   thisRunWeight_ = 0;
 
   product_ = GenLumiInfoProduct(-9999);
 
   filterOnlyEffRun_ = GenFilterInfo(0, 0, 0, 0, 0., 0., 0., 0.);
   hepMCFilterEffRun_ = GenFilterInfo(0, 0, 0, 0, 0., 0., 0., 0.);
 
   xsecBeforeMatching_.clear();
   xsecAfterMatching_.clear();
   jetMatchEffStat_.clear();
   previousLumiBlockLHEXSec_.clear();
   currentLumiBlockLHEXSec_.clear();
 
   return;
 }

void GenXSecAnalyzer::combine	(	GenLumiInfoProduct::XSec &	finalXSec,
		double &	totalw,
		const GenLumiInfoProduct::XSec &	thisRunXSec,
		const double &	thisw
	)

private

Definition at line 256 of file GenXSecAnalyzer.cc.

References GenLumiInfoProduct::XSec::error(), relativeConstraints::error, and GenLumiInfoProduct::XSec::value().

Referenced by endLuminosityBlock(), endRun(), and final_xsec_error().

                                                    {
   double value = finalXSec.value();
   double error = finalXSec.error();
   double thisValue = thisRunXSec.value();
   double thisError = thisRunXSec.error();
   combine(value, error, totalw, thisValue, thisError, thisw);
   finalXSec = GenLumiInfoProduct::XSec(value, error);
   return;
 }

void GenXSecAnalyzer::combine	(	double &	finalValue,
		double &	finalError,
		double &	finalWeight,
		const double &	currentValue,
		const double &	currentError,
		const double &	currentWeight
	)

private

Definition at line 234 of file GenXSecAnalyzer.cc.

References runTheMatrix::err, and mathSSE::sqrt().

                                                            {
   if (finalValue <= 0) {
     finalValue = currentValue;
     finalError = currentError;
     finalWeight += currentWeight;
   } else {
     double wgt1 = (finalError <= 0 || currentError <= 0) ? finalWeight : 1 / (finalError * finalError);
     double wgt2 = (finalError <= 0 || currentError <= 0) ? currentWeight : 1 / (currentError * currentError);
     double xsec = (wgt1 * finalValue + wgt2 * currentValue) / (wgt1 + wgt2);
     double err = (finalError <= 0 || currentError <= 0) ? 0 : 1.0 / std::sqrt(wgt1 + wgt2);
     finalValue = xsec;
     finalError = err;
     finalWeight += currentWeight;
   }
   return;
 }

GenLumiInfoProduct::XSec GenXSecAnalyzer::compute ( const GenLumiInfoProduct & iLumiInfo )

private

Definition at line 269 of file GenXSecAnalyzer.cc.

Referenced by endRun(), and final_xsec_error().

                                                                                    {
   // sum of cross sections and errors over different processes
   double sigSelSum = 0.0;
   double err2SelSum = 0.0;
   double sigSum = 0.0;
   double err2Sum = 0.0;
 
   std::vector<GenLumiInfoProduct::XSec> tempVector_before;
   std::vector<GenLumiInfoProduct::XSec> tempVector_after;
 
   // loop over different processes for each sample
   unsigned int vectorSize = iLumiInfo.getProcessInfos().size();
   for (unsigned int ip = 0; ip < vectorSize; ip++) {
     GenLumiInfoProduct::ProcessInfo proc = iLumiInfo.getProcessInfos()[ip];
     double hepxsec_value = proc.lheXSec().value();
     double hepxsec_error = proc.lheXSec().error() <= 0 ? 0 : proc.lheXSec().error();
     tempVector_before.push_back(GenLumiInfoProduct::XSec(hepxsec_value, hepxsec_error));
 
     sigSelSum += hepxsec_value;
     err2SelSum += hepxsec_error * hepxsec_error;
 
     // skips computation if jet matching efficiency=0
     if (proc.killed().n() < 1) {
       tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
       continue;
     }
 
     // computing jet matching efficiency for this process
     double fracAcc = 0;
     double ntotal = proc.nTotalPos() - proc.nTotalNeg();
     double npass = proc.nPassPos() - proc.nPassNeg();
     switch (hepidwtup_) {
       case 3:
       case -3:
         fracAcc = ntotal > 0 ? npass / ntotal : -1;
         break;
       default:
         fracAcc = proc.selected().sum() > 0 ? proc.killed().sum() / proc.selected().sum() : -1;
         break;
     }
 
     if (fracAcc <= 0) {
       tempVector_after.push_back(GenLumiInfoProduct::XSec(0.0, 0.0));
       continue;
     }
 
     // cross section after matching for this particular process
     double sigmaFin = hepxsec_value * fracAcc;
 
     // computing error on jet matching efficiency
     double relErr = 1.0;
     double efferr2 = 0;
     switch (hepidwtup_) {
       case 3:
       case -3: {
         double ntotal_pos = proc.nTotalPos();
         double effp = ntotal_pos > 0 ? (double)proc.nPassPos() / ntotal_pos : 0;
         double effp_err2 = ntotal_pos > 0 ? (1 - effp) * effp / ntotal_pos : 0;
 
         double ntotal_neg = proc.nTotalNeg();
         double effn = ntotal_neg > 0 ? (double)proc.nPassNeg() / ntotal_neg : 0;
         double effn_err2 = ntotal_neg > 0 ? (1 - effn) * effn / ntotal_neg : 0;
 
         efferr2 = ntotal > 0
                       ? (ntotal_pos * ntotal_pos * effp_err2 + ntotal_neg * ntotal_neg * effn_err2) / ntotal / ntotal
                       : 0;
         break;
       }
       default: {
         double denominator = pow(proc.selected().sum(), 4);
         double passw = proc.killed().sum();
         double passw2 = proc.killed().sum2();
         double failw = proc.selected().sum() - passw;
         double failw2 = proc.selected().sum2() - passw2;
         double numerator = (passw2 * failw * failw + failw2 * passw * passw);
 
         efferr2 = denominator > 0 ? numerator / denominator : 0;
         break;
       }
     }
     double delta2Veto = efferr2 / fracAcc / fracAcc;
 
     // computing total error on cross section after matching efficiency
 
     double sigma2Sum, sigma2Err;
     sigma2Sum = hepxsec_value * hepxsec_value;
     sigma2Err = hepxsec_error * hepxsec_error;
 
     double delta2Sum = delta2Veto + sigma2Err / sigma2Sum;
     relErr = (delta2Sum > 0.0 ? std::sqrt(delta2Sum) : 0.0);
     double deltaFin = sigmaFin * relErr;
 
     tempVector_after.push_back(GenLumiInfoProduct::XSec(sigmaFin, deltaFin));
 
     // sum of cross sections and errors over different processes
     sigSum += sigmaFin;
     err2Sum += deltaFin * deltaFin;
 
   }  // end of loop over different processes
   tempVector_before.push_back(GenLumiInfoProduct::XSec(sigSelSum, sqrt(err2SelSum)));
 
   double total_matcheff = jetMatchEffStat_[10000].filterEfficiency(hepidwtup_);
   double total_matcherr = jetMatchEffStat_[10000].filterEfficiencyError(hepidwtup_);
 
   double xsec_after = sigSelSum * total_matcheff;
   double xsecerr_after = (total_matcheff > 0 && sigSelSum > 0)
                              ? xsec_after * sqrt(err2SelSum / sigSelSum / sigSelSum +
                                                  total_matcherr * total_matcherr / total_matcheff / total_matcheff)
                              : 0;
 
   GenLumiInfoProduct::XSec result(xsec_after, xsecerr_after);
   tempVector_after.push_back(result);
 
   xsecBeforeMatching_ = tempVector_before;
   xsecAfterMatching_ = tempVector_after;
 
   return result;
 }

void GenXSecAnalyzer::endJob ( void )

overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 388 of file GenXSecAnalyzer.cc.

References MillePedeFileConverter_cfg::e, SiPixelPhase1Clusters_cfi::e3, GenLumiInfoProduct::XSec::error(), GenFilterInfo::filterEfficiency(), GenFilterInfo::filterEfficiencyError(), filterOnlyEffStat_, alignBH_cfg::fixed, hepidwtup_, hepMCFilterEffStat_, mps_fire::i, jetMatchEffStat_, dqmdumpme::last, nMCs_, GenFilterInfo::numEventsPassed(), GenFilterInfo::numEventsTotal(), GenFilterInfo::numPassNegativeEvents(), GenFilterInfo::numPassPositiveEvents(), GenFilterInfo::numTotalNegativeEvents(), GenFilterInfo::numTotalPositiveEvents(), funct::pow(), mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, GenFilterInfo::sumPassWeights(), GenFilterInfo::sumWeights(), GenFilterInfo::sumWeights2(), overlapproblemtsosanalyzer_cfi::title, GenLumiInfoProduct::XSec::value(), xsec_, xsecAfterMatching_, xsecBeforeMatching_, and xsecPreFilter_.

Referenced by final_xsec_error().

                              {
   edm::LogPrint("GenXSecAnalyzer") << "\n"
                                    << "------------------------------------"
                                    << "\n"
                                    << "GenXsecAnalyzer:"
                                    << "\n"
                                    << "------------------------------------";
 
   if (jetMatchEffStat_.empty()) {
     edm::LogPrint("GenXSecAnalyzer") << "------------------------------------"
                                      << "\n"
                                      << "Cross-section summary not available"
                                      << "\n"
                                      << "------------------------------------";
     return;
   }
 
   // fraction of negative weights
   double final_fract_neg_w = 0;
   double final_fract_neg_w_unc = 0;
 
   // below print out is only for combination of same physics MC samples and ME+Pythia MCs
 
   if (nMCs_ == 1 && hepidwtup_ != -1) {
     edm::LogPrint("GenXSecAnalyzer")
         << "-----------------------------------------------------------------------------------------------------------"
            "--------------------------------------------------------------- \n"
         << "Overall cross-section summary \n"
         << "-----------------------------------------------------------------------------------------------------------"
            "---------------------------------------------------------------";
     edm::LogPrint("GenXSecAnalyzer") << "Process\t\txsec_before [pb]\t\tpassed\tnposw\tnnegw\ttried\tnposw\tnnegw "
                                         "\txsec_match [pb]\t\t\taccepted [%]\t event_eff [%]";
 
     const unsigned sizeOfInfos = jetMatchEffStat_.size();
     const unsigned last = sizeOfInfos - 1;
     std::string* title = new std::string[sizeOfInfos];
     unsigned int i = 0;
     double jetmatch_eff = 0;
     double jetmatch_err = 0;
     double matching_eff = 1;
     double matching_efferr = 1;
 
     for (std::map<int, GenFilterInfo>::const_iterator iter = jetMatchEffStat_.begin(); iter != jetMatchEffStat_.end();
          ++iter, i++) {
       GenFilterInfo thisJetMatchStat = iter->second;
       GenFilterInfo thisEventEffStat =
           GenFilterInfo(thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
                         0,
                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents(),
                         0,
                         thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
                         thisJetMatchStat.numPassPositiveEvents() + thisJetMatchStat.numPassNegativeEvents(),
                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents(),
                         thisJetMatchStat.numTotalPositiveEvents() + thisJetMatchStat.numTotalNegativeEvents());
 
       jetmatch_eff = thisJetMatchStat.filterEfficiency(hepidwtup_);
       jetmatch_err = thisJetMatchStat.filterEfficiencyError(hepidwtup_);
 
       if (i == last) {
         title[i] = "Total";
 
         edm::LogPrint("GenXSecAnalyzer")
             << "-------------------------------------------------------------------------------------------------------"
                "------------------------------------------------------------------- ";
 
         // fill negative fraction of negative weights and uncertainty after matching
         final_fract_neg_w = thisEventEffStat.numEventsPassed() > 0
                                 ? thisJetMatchStat.numPassNegativeEvents() / thisEventEffStat.numEventsPassed()
                                 : 0;
         final_fract_neg_w_unc =
             thisJetMatchStat.numPassNegativeEvents() > 0
                 ? final_fract_neg_w * final_fract_neg_w / thisEventEffStat.numEventsPassed() *
                       sqrt(thisJetMatchStat.numPassPositiveEvents() * thisJetMatchStat.numPassPositiveEvents() /
                                thisJetMatchStat.numPassNegativeEvents() +
                            thisJetMatchStat.numPassPositiveEvents())
                 : 0;
       } else {
         title[i] = Form("%d", i);
       }
 
       edm::LogPrint("GenXSecAnalyzer") << title[i] << "\t\t" << std::scientific << std::setprecision(3)
                                        << xsecBeforeMatching_[i].value() << " +/- " << xsecBeforeMatching_[i].error()
                                        << "\t\t" << thisEventEffStat.numEventsPassed() << "\t"
                                        << thisJetMatchStat.numPassPositiveEvents() << "\t"
                                        << thisJetMatchStat.numPassNegativeEvents() << "\t"
                                        << thisEventEffStat.numEventsTotal() << "\t"
                                        << thisJetMatchStat.numTotalPositiveEvents() << "\t"
                                        << thisJetMatchStat.numTotalNegativeEvents() << "\t" << std::scientific
                                        << std::setprecision(3) << xsecAfterMatching_[i].value() << " +/- "
                                        << xsecAfterMatching_[i].error() << "\t\t" << std::fixed << std::setprecision(1)
                                        << (jetmatch_eff * 100) << " +/- " << (jetmatch_err * 100) << "\t" << std::fixed
                                        << std::setprecision(1) << (thisEventEffStat.filterEfficiency(+3) * 100)
                                        << " +/- " << (thisEventEffStat.filterEfficiencyError(+3) * 100);
 
       matching_eff = thisEventEffStat.filterEfficiency(+3);
       matching_efferr = thisEventEffStat.filterEfficiencyError(+3);
     }
     delete[] title;
 
     edm::LogPrint("GenXSecAnalyzer")
         << "-----------------------------------------------------------------------------------------------------------"
            "---------------------------------------------------------------";
 
     edm::LogPrint("GenXSecAnalyzer") << "Before matching: total cross section = " << std::scientific
                                      << std::setprecision(3) << xsecBeforeMatching_[last].value() << " +- "
                                      << xsecBeforeMatching_[last].error() << " pb";
 
     edm::LogPrint("GenXSecAnalyzer") << "After matching: total cross section = " << std::scientific
                                      << std::setprecision(3) << xsecAfterMatching_[last].value() << " +- "
                                      << xsecAfterMatching_[last].error() << " pb";
 
     edm::LogPrint("GenXSecAnalyzer") << "Matching efficiency = " << std::fixed << std::setprecision(1) << matching_eff
                                      << " +/- " << matching_efferr << "   [TO BE USED IN MCM]";
 
   } else if (hepidwtup_ == -1)
     edm::LogPrint("GenXSecAnalyzer") << "Before Filter: total cross section = " << std::scientific
                                      << std::setprecision(3) << xsecPreFilter_.value() << " +- "
                                      << xsecPreFilter_.error() << " pb";
 
   // hepMC filter efficiency
   double hepMCFilter_eff = 1.0;
   double hepMCFilter_err = 0.0;
   if (hepMCFilterEffStat_.sumWeights2() > 0) {
     hepMCFilter_eff = hepMCFilterEffStat_.filterEfficiency(-1);
     hepMCFilter_err = hepMCFilterEffStat_.filterEfficiencyError(-1);
     edm::LogPrint("GenXSecAnalyzer") << "HepMC filter efficiency (taking into account weights)= "
                                      << "(" << hepMCFilterEffStat_.sumPassWeights() << ")"
                                      << " / "
                                      << "(" << hepMCFilterEffStat_.sumWeights() << ")"
                                      << " = " << std::scientific << std::setprecision(3) << hepMCFilter_eff << " +- "
                                      << hepMCFilter_err;
 
     double hepMCFilter_event_total =
         hepMCFilterEffStat_.numTotalPositiveEvents() + hepMCFilterEffStat_.numTotalNegativeEvents();
     double hepMCFilter_event_pass =
         hepMCFilterEffStat_.numPassPositiveEvents() + hepMCFilterEffStat_.numPassNegativeEvents();
     double hepMCFilter_event_eff = hepMCFilter_event_total > 0 ? hepMCFilter_event_pass / hepMCFilter_event_total : 0;
     double hepMCFilter_event_err =
         hepMCFilter_event_total > 0
             ? sqrt((1 - hepMCFilter_event_eff) * hepMCFilter_event_eff / hepMCFilter_event_total)
             : -1;
     edm::LogPrint("GenXSecAnalyzer") << "HepMC filter efficiency (event-level)= "
                                      << "(" << hepMCFilter_event_pass << ")"
                                      << " / "
                                      << "(" << hepMCFilter_event_total << ")"
                                      << " = " << std::scientific << std::setprecision(3) << hepMCFilter_event_eff
                                      << " +- " << hepMCFilter_event_err;
   }
 
   // gen-particle filter efficiency
   if (filterOnlyEffStat_.sumWeights2() > 0) {
     double filterOnly_eff = filterOnlyEffStat_.filterEfficiency(-1);
     double filterOnly_err = filterOnlyEffStat_.filterEfficiencyError(-1);
 
     edm::LogPrint("GenXSecAnalyzer") << "Filter efficiency (taking into account weights)= "
                                      << "(" << filterOnlyEffStat_.sumPassWeights() << ")"
                                      << " / "
                                      << "(" << filterOnlyEffStat_.sumWeights() << ")"
                                      << " = " << std::scientific << std::setprecision(3) << filterOnly_eff << " +- "
                                      << filterOnly_err;
 
     double filterOnly_event_total =
         filterOnlyEffStat_.numTotalPositiveEvents() + filterOnlyEffStat_.numTotalNegativeEvents();
     double filterOnly_event_pass =
         filterOnlyEffStat_.numPassPositiveEvents() + filterOnlyEffStat_.numPassNegativeEvents();
     double filterOnly_event_eff = filterOnly_event_total > 0 ? filterOnly_event_pass / filterOnly_event_total : 0;
     double filterOnly_event_err = filterOnly_event_total > 0
                                       ? sqrt((1 - filterOnly_event_eff) * filterOnly_event_eff / filterOnly_event_total)
                                       : -1;
     edm::LogPrint("GenXSecAnalyzer") << "Filter efficiency (event-level)= "
                                      << "(" << filterOnly_event_pass << ")"
                                      << " / "
                                      << "(" << filterOnly_event_total << ")"
                                      << " = " << std::scientific << std::setprecision(3) << filterOnly_event_eff
                                      << " +- " << filterOnly_event_err << "    [TO BE USED IN MCM]";
 
     // fill negative fraction of negative weights and uncertainty after filter
     final_fract_neg_w =
         filterOnly_event_pass > 0 ? filterOnlyEffStat_.numPassNegativeEvents() / (filterOnly_event_pass) : 0;
     final_fract_neg_w_unc =
         filterOnlyEffStat_.numPassNegativeEvents() > 0
             ? final_fract_neg_w * final_fract_neg_w / filterOnly_event_pass *
                   sqrt(filterOnlyEffStat_.numPassPositiveEvents() * filterOnlyEffStat_.numPassPositiveEvents() /
                            filterOnlyEffStat_.numPassNegativeEvents() +
                        filterOnlyEffStat_.numPassPositiveEvents())
             : 0;
   }
 
   edm::LogPrint("GenXSecAnalyzer") << "\nAfter filter: final cross section = " << std::scientific
                                    << std::setprecision(3) << xsec_.value() << " +- " << xsec_.error() << " pb";
 
   edm::LogPrint("GenXSecAnalyzer") << "After filter: final fraction of events with negative weights = "
                                    << std::scientific << std::setprecision(3) << final_fract_neg_w << " +- "
                                    << final_fract_neg_w_unc;
 
   // L=[N*(1-2f)^2]/s
   double lumi_1M_evts =
       xsec_.value() > 0 ? 1e6 * (1 - 2 * final_fract_neg_w) * (1 - 2 * final_fract_neg_w) / xsec_.value() / 1e3 : 0;
   double lumi_1M_evts_unc =
       xsec_.value() > 0 ? (1 - 2 * final_fract_neg_w) * lumi_1M_evts *
                               sqrt(1e-6 + 16 * pow(final_fract_neg_w_unc, 2) / pow(1 - 2 * final_fract_neg_w, 2) +
                                    pow(xsec_.error() / xsec_.value(), 2))
                         : 0;
   edm::LogPrint("GenXSecAnalyzer") << "After filter: final equivalent lumi for 1M events (1/fb) = " << std::scientific
                                    << std::setprecision(3) << lumi_1M_evts << " +- " << lumi_1M_evts_unc;
 }

void GenXSecAnalyzer::endLuminosityBlock	(	edm::LuminosityBlock const &	iLumi,
		edm::EventSetup const &
	)

overrideprivate

Definition at line 69 of file GenXSecAnalyzer.cc.

References combine(), currentLumiBlockLHEXSec_, runTheMatrix::err, GenLumiInfoProduct::XSec::error(), filterOnlyEffRun_, filterOnlyEffStat_, genFilterInfoToken_, rivetAnalyzer_cfi::genLumiInfo, genLumiInfoToken_, edm::LuminosityBlock::getByToken(), GenLumiInfoProduct::getHEPIDWTUP(), GenLumiInfoProduct::getProcessInfos(), hepidwtup_, hepMCFilterEffRun_, hepMCFilterEffStat_, hepMCFilterInfoToken_, triggerObjects_cff::id, edm::HandleBase::isValid(), jetMatchEffStat_, GenFilterInfo::mergeProduct(), GenLumiInfoProduct::mergeProduct(), previousLumiBlockLHEXSec_, product_, GenLumiInfoProduct::FinalStat::sum(), GenLumiInfoProduct::FinalStat::sum2(), GenFilterInfo::sumPassWeights(), thisRunWeight_, thisRunWeightPre_, GenLumiInfoProduct::XSec::value(), relativeConstraints::value, x, and y.

Referenced by final_xsec_error().

                                                                                               {
   edm::Handle<GenLumiInfoProduct> genLumiInfo;
   iLumi.getByToken(genLumiInfoToken_, genLumiInfo);
   if (!genLumiInfo.isValid())
     return;
   hepidwtup_ = genLumiInfo->getHEPIDWTUP();
 
   std::vector<GenLumiInfoProduct::ProcessInfo> theProcesses = genLumiInfo->getProcessInfos();
 
   unsigned int theProcesses_size = theProcesses.size();
 
   // if it's a pure parton-shower generator, check there should be only one element in thisProcessInfos
   // the error of lheXSec is -1
   if (hepidwtup_ == -1) {
     if (theProcesses_size != 1) {
       edm::LogError("GenXSecAnalyzer::endLuminosityBlock") << "Pure parton shower has thisProcessInfos size!=1";
       return;
     }
   }
 
   for (unsigned int ip = 0; ip < theProcesses_size; ip++) {
     if (theProcesses[ip].lheXSec().value() < 0) {
       edm::LogError("GenXSecAnalyzer::endLuminosityBlock")
           << "cross section of process " << ip << " value = " << theProcesses[ip].lheXSec().value();
       return;
     }
   }
 
   product_.mergeProduct(*genLumiInfo);
 
   edm::Handle<GenFilterInfo> genFilter;
   iLumi.getByToken(genFilterInfoToken_, genFilter);
   if (genFilter.isValid()) {
     filterOnlyEffStat_.mergeProduct(*genFilter);
     filterOnlyEffRun_.mergeProduct(*genFilter);
     thisRunWeight_ += genFilter->sumPassWeights();
   }
 
   edm::Handle<GenFilterInfo> hepMCFilter;
   iLumi.getByToken(hepMCFilterInfoToken_, hepMCFilter);
 
   if (hepMCFilter.isValid()) {
     hepMCFilterEffStat_.mergeProduct(*hepMCFilter);
     hepMCFilterEffRun_.mergeProduct(*hepMCFilter);
   }
 
   // doing generic summing for jet matching statistics
   // and computation of combined LHE information
   for (unsigned int ip = 0; ip < theProcesses_size; ip++) {
     int id = theProcesses[ip].process();
     GenFilterInfo& x = jetMatchEffStat_[id];
     GenLumiInfoProduct::XSec& y = currentLumiBlockLHEXSec_[id];
     GenLumiInfoProduct::FinalStat temp_killed = theProcesses[ip].killed();
     GenLumiInfoProduct::FinalStat temp_selected = theProcesses[ip].selected();
     double passw = temp_killed.sum();
     double passw2 = temp_killed.sum2();
     double totalw = temp_selected.sum();
     double totalw2 = temp_selected.sum2();
     GenFilterInfo tempInfo(theProcesses[ip].nPassPos(),
                            theProcesses[ip].nPassNeg(),
                            theProcesses[ip].nTotalPos(),
                            theProcesses[ip].nTotalNeg(),
                            passw,
                            passw2,
                            totalw,
                            totalw2);
 
     // matching statistics for all processes
     jetMatchEffStat_[10000].mergeProduct(tempInfo);
     double currentValue = theProcesses[ip].lheXSec().value();
     double currentError = theProcesses[ip].lheXSec().error();
 
     // this process ID has occurred before
     if (y.value() > 0) {
       x.mergeProduct(tempInfo);
       double previousValue = previousLumiBlockLHEXSec_[id].value();
 
       if (currentValue != previousValue)  // transition of cross section
       {
         double xsec = y.value();
         double err = y.error();
         combine(xsec, err, thisRunWeightPre_, currentValue, currentError, totalw);
         y = GenLumiInfoProduct::XSec(xsec, err);
       } else  // LHE cross section is the same as previous lumiblock
         thisRunWeightPre_ += totalw;
 
     }
     // this process ID has never occurred before
     else {
       x = tempInfo;
       y = theProcesses[ip].lheXSec();
       thisRunWeightPre_ += totalw;
     }
 
     previousLumiBlockLHEXSec_[id] = theProcesses[ip].lheXSec();
   }  // end
 
   return;
 }

void GenXSecAnalyzer::endRun	(	edm::Run const &	iRun,
		edm::EventSetup const &
	)

overrideprivate

Definition at line 169 of file GenXSecAnalyzer.cc.

References combine(), compute(), gather_cfg::cout, currentLumiBlockLHEXSec_, GenLumiInfoProduct::XSec::error(), GenFilterInfo::filterEfficiency(), GenFilterInfo::filterEfficiencyError(), filterOnlyEffRun_, alignBH_cfg::fixed, edm::Run::getByToken(), GenLumiInfoProduct::getProcessInfos(), hepidwtup_, hepMCFilterEffRun_, LHERunInfoProduct::heprup(), mps_fire::i, lheRunInfoToken_, lhef::HEPRUP::LPRUP, Max(), funct::pow(), product_, writedatasetfile::run, GenLumiInfoProduct::ProcessInfo::setLheXSec(), GenLumiInfoProduct::setProcessInfo(), mathSSE::sqrt(), GenFilterInfo::sumWeights2(), groupFilesInBlocks::temp, thisRunWeight_, thisRunWeightPre_, totalWeight_, totalWeightPre_, GenLumiInfoProduct::XSec::value(), lhef::HEPRUP::XERRUP, lhef::HEPRUP::XMAXUP, xsec_, xsecPreFilter_, and lhef::HEPRUP::XSECUP.

Referenced by final_xsec_error().

                                                                      {
   //xsection before matching
   edm::Handle<LHERunInfoProduct> run;
 
   if (iRun.getByToken(lheRunInfoToken_, run)) {
     const lhef::HEPRUP thisHeprup_ = run->heprup();
 
     for (unsigned int iSize = 0; iSize < thisHeprup_.XSECUP.size(); iSize++) {
       std::cout << std::setw(14) << std::fixed << thisHeprup_.XSECUP[iSize] << std::setw(14) << std::fixed
                 << thisHeprup_.XERRUP[iSize] << std::setw(14) << std::fixed << thisHeprup_.XMAXUP[iSize]
                 << std::setw(14) << std::fixed << thisHeprup_.LPRUP[iSize] << std::endl;
     }
     std::cout << " " << std::endl;
   }
 
   // compute cross section for this run first
   // set the correct combined LHE+filter cross sections
   unsigned int i = 0;
   std::vector<GenLumiInfoProduct::ProcessInfo> newInfos;
   for (std::map<int, GenLumiInfoProduct::XSec>::const_iterator iter = currentLumiBlockLHEXSec_.begin();
        iter != currentLumiBlockLHEXSec_.end();
        ++iter, i++) {
     GenLumiInfoProduct::ProcessInfo temp = product_.getProcessInfos()[i];
     temp.setLheXSec(iter->second.value(), iter->second.error());
     newInfos.push_back(temp);
   }
   product_.setProcessInfo(newInfos);
 
   const GenLumiInfoProduct::XSec thisRunXSecPre = compute(product_);
   // xsection after matching before filters
   combine(xsecPreFilter_, totalWeightPre_, thisRunXSecPre, thisRunWeightPre_);
 
   double thisHepFilterEff = 1;
   double thisHepFilterErr = 0;
 
   if (hepMCFilterEffRun_.sumWeights2() > 0) {
     thisHepFilterEff = hepMCFilterEffRun_.filterEfficiency(hepidwtup_);
     thisHepFilterErr = hepMCFilterEffRun_.filterEfficiencyError(hepidwtup_);
     if (thisHepFilterEff < 0) {
       thisHepFilterEff = 1;
       thisHepFilterErr = 0;
     }
   }
 
   double thisGenFilterEff = 1;
   double thisGenFilterErr = 0;
 
   if (filterOnlyEffRun_.sumWeights2() > 0) {
     thisGenFilterEff = filterOnlyEffRun_.filterEfficiency(hepidwtup_);
     thisGenFilterErr = filterOnlyEffRun_.filterEfficiencyError(hepidwtup_);
     if (thisGenFilterEff < 0) {
       thisGenFilterEff = 1;
       thisGenFilterErr = 0;
     }
   }
   double thisXsec = thisRunXSecPre.value() > 0 ? thisHepFilterEff * thisGenFilterEff * thisRunXSecPre.value() : 0;
   double thisErr =
       thisRunXSecPre.value() > 0
           ? thisXsec * sqrt(pow(TMath::Max(thisRunXSecPre.error(), (double)0) / thisRunXSecPre.value(), 2) +
                             pow(thisHepFilterErr / thisHepFilterEff, 2) + pow(thisGenFilterErr / thisGenFilterEff, 2))
           : 0;
   const GenLumiInfoProduct::XSec thisRunXSec = GenLumiInfoProduct::XSec(thisXsec, thisErr);
   combine(xsec_, totalWeight_, thisRunXSec, thisRunWeight_);
 }