#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 ()

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

Public Member Functions inherited from edm::one::EDAnalyzerBase
	EDAnalyzerBase ()

ModuleDescription const &	moduleDescription () const

virtual	~EDAnalyzerBase ()

Public Member Functions inherited from edm::EDConsumerBase
	EDConsumerBase ()

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

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

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

std::vector < ProductHolderIndexAndSkipBit > const &	itemsToGetFromEvent () const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const

bool	registeredToConsume (ProductHolderIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

void	updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)

virtual	~EDConsumerBase ()

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

virtual void	beginJob () override

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

void	compute ()

virtual void	endJob () override

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

Private Attributes
edm::EDGetTokenT< GenFilterInfo >	genFilterInfoToken_

edm::EDGetTokenT < GenLumiInfoProduct >	genLumiInfoToken_

int	hepidwtup_

GenFilterInfo	jetMatchEffStat_

std::vector< GenLumiInfoProduct >	products_

GenFilterInfo	totalEffStat_

GenLumiInfoProduct::XSec	xsec_

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

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::one::EDAnalyzerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

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 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 30 of file GenXSecAnalyzer.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 6 of file GenXSecAnalyzer.cc.

References genFilterInfoToken_, genLumiInfoToken_, and products_.

                                                               :
   hepidwtup_(-1),
   xsec_(0,0),
   jetMatchEffStat_(0,0,0,0,0.,0.,0.,0.),
   totalEffStat_(0,0,0,0,0.,0.,0.,0.)
 {
   products_.clear();
   
   genFilterInfoToken_ = consumes<GenFilterInfo,edm::InLumi>(edm::InputTag("genFilterEfficiencyProducer",""));
   genLumiInfoToken_ = consumes<GenLumiInfoProduct,edm::InLumi>(edm::InputTag("generator",""));
 }

GenXSecAnalyzer::~GenXSecAnalyzer ( )

Definition at line 18 of file GenXSecAnalyzer.cc.

19 {

20 }

Member Function Documentation

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

overrideprivatevirtual

Implements edm::one::EDAnalyzerBase.

Definition at line 28 of file GenXSecAnalyzer.cc.

29 {

30 }

void GenXSecAnalyzer::beginJob ( void )

overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 23 of file GenXSecAnalyzer.cc.

References products_.

                           {
   products_.clear();  
 }

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

overrideprivatevirtual

Definition at line 34 of file GenXSecAnalyzer.cc.

34 {

35 }

void GenXSecAnalyzer::compute ( )

private

Definition at line 109 of file GenXSecAnalyzer.cc.

Referenced by endJob().

 {
   // for pure parton shower generator
 
   if(hepidwtup_== -1)
     {
       double sigSum = 0.0;
       double totalN = 0.0;
       for(unsigned int i=0; i < products_.size(); i++){
 
     GenLumiInfoProduct::ProcessInfo proc = products_[i].getProcessInfos()[0];     
     double hepxsec_value = proc.lheXSec().value();
 
     sigSum += proc.tried().sum() * hepxsec_value;
     totalN += proc.tried().sum();
       }
       double sigAve = totalN>1e-6? sigSum/totalN: 0;
       xsec_ = GenLumiInfoProduct::XSec(sigAve,-1);      
     }
   // for ME+parton shower MC
   else{
 
     double sum_numerator = 0;
     double sum_denominator = 0;
   
     for(unsigned int i=0; i < products_.size(); i++){
 
       double sigSelSum = 0.0;
       double sigSum = 0.0;
       double sigBrSum = 0.0;
       double err2Sum = 0.0;
       double errBr2Sum = 0.0;
 
       unsigned int vectorSize = products_[i].getProcessInfos().size();
       for(unsigned int ip=0; ip < vectorSize; ip++){
     GenLumiInfoProduct::ProcessInfo proc = products_[i].getProcessInfos()[ip];    
     double hepxsec_value = proc.lheXSec().value();
     double hepxsec_error = proc.lheXSec().error();
 
     if (!proc.killed().n())
       continue;
 
     double sigma2Sum, sigma2Err;
     sigma2Sum = hepxsec_value * hepxsec_value;
     sigma2Err = hepxsec_error * hepxsec_error;
 
     double sigmaAvg = hepxsec_value;
 
     double fracAcc = 0;
     double ntotal = proc.nTotalPos()-proc.nTotalNeg();
     double npass  = proc.nPassPos() -proc.nPassNeg();
     switch(hepidwtup_){
     case 3: case -3:
       fracAcc = ntotal > 1e-6? npass/ntotal: -1;
         break;
     default:
       fracAcc = proc.selected().sum() > 1e-6? proc.killed().sum() / proc.selected().sum():-1;
       break;
     }
 
     if(fracAcc<1e-6)continue;
 
     double fracBr = proc.accepted().sum() > 0.0 ?
       proc.acceptedBr().sum() / proc.accepted().sum() : 1;
     double sigmaFin = sigmaAvg * fracAcc * fracBr;
     double sigmaFinBr = sigmaFin * fracBr;
 
     double relErr = 1.0;
     if (proc.killed().n() > 1) {
       double efferr2=0;
       switch(hepidwtup_) {
       case 3: case -3:
         {
           double ntotal_pos = proc.nTotalPos();
           double effp  = ntotal_pos > 1e-6?
         (double)proc.nPassPos()/ntotal_pos:0;
           double effp_err2 = ntotal_pos > 1e-6?
         (1-effp)*effp/ntotal_pos: 0;
 
           double ntotal_neg = proc.nTotalNeg();
           double effn  = ntotal_neg > 1e-6?
         (double)proc.nPassNeg()/ntotal_neg:0;
           double effn_err2 = ntotal_neg > 1e-6?
         (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>1e-6?
         numerator/denominator:0;
           break;
         }
       }
       double delta2Veto = efferr2/fracAcc/fracAcc;
       double delta2Sum = delta2Veto
         + sigma2Err / sigma2Sum;
       relErr = (delta2Sum > 0.0 ?
             std::sqrt(delta2Sum) : 0.0);
     }
     double deltaFin = sigmaFin * relErr;
     double deltaFinBr = sigmaFinBr * relErr;
 
     sigSelSum += sigmaAvg;
     sigSum += sigmaFin;
     sigBrSum += sigmaFinBr;
     err2Sum += deltaFin * deltaFin;
     errBr2Sum += deltaFinBr * deltaFinBr;
 
 
 
       } // end of loop over different processes
       
       double dN = std::sqrt(errBr2Sum);
       sum_denominator  +=  (dN> 1e-6)? 1/dN/dN: 0;
       sum_numerator    +=  (dN> 1e-6)? sigBrSum/dN/dN: 0;
 
 
     } // end of loop over different samples
     double final_value = sum_denominator > 0? sum_numerator/sum_denominator : 0;
     double final_error = sum_denominator > 0? 1/sqrt(sum_denominator) : -1;
     xsec_ = GenLumiInfoProduct::XSec(final_value, final_error);
   }
   return;
 }

void GenXSecAnalyzer::endJob ( void )

overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 247 of file GenXSecAnalyzer.cc.

References compute(), GenLumiInfoProduct::XSec::error(), GenFilterInfo::filterEfficiency(), GenFilterInfo::filterEfficiencyError(), hepidwtup_, jetMatchEffStat_, products_, mathSSE::sqrt(), GenFilterInfo::sumPassWeights(), GenFilterInfo::sumWeights(), totalEffStat_, GenLumiInfoProduct::XSec::value(), and xsec_.

                         {
 
   edm::LogPrint("GenXSecAnalyzer") << "\n"
   << "------------------------------------" << "\n"
   << "GenXsecAnalyzer:" << "\n"
   << "------------------------------------";
   
   if(!products_.size()) {
     edm::LogPrint("GenXSecAnalyzer") << "------------------------------------" << "\n"
     << "Cross-section summary not available" << "\n"
     << "------------------------------------";
     return;
   }
   
   
   compute();
 
   double filterOnly_eff = totalEffStat_.filterEfficiency(hepidwtup_);
   double filterOnly_err = totalEffStat_.filterEfficiencyError(hepidwtup_);
   
   double jetmatching_eff_total = jetMatchEffStat_.filterEfficiency(hepidwtup_);
   double jetmatching_err_total = jetMatchEffStat_.filterEfficiencyError(hepidwtup_);
 
   edm::LogPrint("GenXSecAnalyzer") << "------------------------------------" << "\n"
   << "Overall cross-section summary:" << "\n"
   << "------------------------------------";
   
   edm::LogPrint("GenXSecAnalyzer") << "jet matching efficiency = " << 
     jetMatchEffStat_.sumPassWeights() << "/" << 
     jetMatchEffStat_.sumWeights() << " = " 
        << std::setprecision(6) << jetmatching_eff_total << " +- " << jetmatching_err_total << "\n";
 
   edm::LogPrint("GenXSecAnalyzer") << "Before filter: cross section = " << std::setprecision(6)  << xsec_.value() << " +- " << std::setprecision(6) << xsec_.error() <<  " pb";
 
   edm::LogPrint("GenXSecAnalyzer") << "Filter efficiency = " << std::setprecision(6)  
         << filterOnly_eff << " +- " << filterOnly_err;
 
 
   double xsec_after  = xsec_.value()*filterOnly_eff ;
   double error_after = xsec_after*sqrt(xsec_.error()*xsec_.error()/xsec_.value()/xsec_.value()+
                   filterOnly_err*filterOnly_err/filterOnly_eff/filterOnly_eff);
 
   edm::LogPrint("GenXSecAnalyzer") << "After filter: cross section = " 
         << std::setprecision(6) << xsec_after
         << " +- " 
         << std::setprecision(6) << error_after
         << " pb";
 
 }

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

overrideprivatevirtual

Definition at line 38 of file GenXSecAnalyzer.cc.

References alignCSCRings::e, genFilterInfoToken_, genLumiInfoToken_, edm::LuminosityBlock::getByToken(), hepidwtup_, i, edm::HandleBase::isValid(), jetMatchEffStat_, GenFilterInfo::mergeProduct(), products_, GenLumiInfoProduct::FinalStat::sum(), GenLumiInfoProduct::FinalStat::sum2(), totalEffStat_, and relativeConstraints::value.

                                                                                          {
 
   edm::Handle<GenFilterInfo> genFilter;
   iLumi.getByToken(genFilterInfoToken_,genFilter);
   if(genFilter.isValid())
     totalEffStat_.mergeProduct(*genFilter);
 
 
   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;
       }
       if(theProcesses[0].lheXSec().value()<1e-6){
     edm::LogError("GenXSecAnalyzer::endLuminosityBlock") << "cross section value = "  << theProcesses[0].lheXSec().value();
     return;
       }
     }
 
   // doing generic summing 
   for(unsigned int ip=0; ip < theProcesses_size; ip++)
     {
       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();
       jetMatchEffStat_.mergeProduct(GenFilterInfo(
                           theProcesses[ip].nPassPos(),
                           theProcesses[ip].nPassNeg(),
                           theProcesses[ip].nTotalPos(),
                           theProcesses[ip].nTotalNeg(),
                           passw,
                           passw2,
                           totalw,
                           totalw2)
                     );
 
     }
   // now determine if this LuminosityBlock has the same lheXSec as existing products
   bool sameMC = false;
   for(unsigned int i=0; i < products_.size(); i++){
 
     if(products_[i].mergeProduct(*genLumiInfo))
       sameMC = true;
     else if(!products_[i].samePhysics(*genLumiInfo))
       {
     edm::LogError("GenXSecAnalyzer::endLuminosityBlock") << "Merging samples that come from different physics processes";
     return;
       }
 
   }
   
   if(!sameMC)
       products_.push_back(*genLumiInfo);
   return;
 }