GenXSecAnalyzer Class Reference

#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::WatchRuns, edm::one::WatchLuminosityBlocks >
	EDAnalyzer ()=default
Public Member Functions inherited from edm::one::EDAnalyzerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
	EDAnalyzerBase ()
ModuleDescription const &	moduleDescription () const
virtual	~EDAnalyzerBase ()
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
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
void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)
virtual	~EDConsumerBase () noexcept(false)

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
virtual 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 &)
virtual void	endJob () override
virtual void	endLuminosityBlock (edm::LuminosityBlock const &, edm::EventSetup const &) override
virtual 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 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 34 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_, 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

(

)

Definition at line 34 of file GenXSecAnalyzer.cc.

{
}

Member Function Documentation

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

overrideprivatevirtual

Definition at line 84 of file GenXSecAnalyzer.cc.

Referenced by final_xsec_error().

{
}

void GenXSecAnalyzer::beginJob ( void  )

overrideprivatevirtual

Reimplemented from edm::one::EDAnalyzerBase.

Definition at line 39 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 &    )

overrideprivatevirtual

Definition at line 77 of file GenXSecAnalyzer.cc.

Referenced by final_xsec_error().

                                                                                           {

}

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

overrideprivatevirtual

Definition at line 50 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 319 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 290 of file GenXSecAnalyzer.cc.

References 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 332 of file GenXSecAnalyzer.cc.

References cuy::denominator, GenLumiInfoProduct::XSec::error(), GenLumiInfoProduct::getProcessInfos(), hepidwtup_, jetMatchEffStat_, GenLumiInfoProduct::ProcessInfo::killed(), GenLumiInfoProduct::ProcessInfo::lheXSec(), GenLumiInfoProduct::FinalStat::n(), GenLumiInfoProduct::ProcessInfo::nPassNeg(), GenLumiInfoProduct::ProcessInfo::nPassPos(), groupFilesInBlocks::ntotal, GenLumiInfoProduct::ProcessInfo::nTotalNeg(), GenLumiInfoProduct::ProcessInfo::nTotalPos(), cuy::numerator, funct::pow(), proc, mps_fire::result, GenLumiInfoProduct::ProcessInfo::selected(), mathSSE::sqrt(), GenLumiInfoProduct::FinalStat::sum(), GenLumiInfoProduct::FinalStat::sum2(), GenLumiInfoProduct::XSec::value(), xsecAfterMatching_, and xsecBeforeMatching_.

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 464 of file GenXSecAnalyzer.cc.

References GenLumiInfoProduct::XSec::error(), GenFilterInfo::filterEfficiency(), GenFilterInfo::filterEfficiencyError(), filterOnlyEffStat_, alignBH_cfg::fixed, hepidwtup_, hepMCFilterEffStat_, mps_fire::i, jetMatchEffStat_, plotBeamSpotDB::last, nMCs_, GenFilterInfo::numEventsPassed(), GenFilterInfo::numEventsTotal(), GenFilterInfo::numPassNegativeEvents(), GenFilterInfo::numPassPositiveEvents(), GenFilterInfo::numTotalNegativeEvents(), GenFilterInfo::numTotalPositiveEvents(), mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, GenFilterInfo::sumPassWeights(), GenFilterInfo::sumWeights(), GenFilterInfo::sumWeights2(), fftjetcommon_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 ? lumi_1M_evts*lumi_1M_evts*sqrt(16/(1-2*final_fract_neg_w)/(1-2*final_fract_neg_w)*final_fract_neg_w_unc*final_fract_neg_w_unc + xsec_.error()*xsec_.error()/xsec_.value()/xsec_.value()) : 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 &    )

overrideprivatevirtual

Definition at line 91 of file GenXSecAnalyzer.cc.

References combine(), currentLumiBlockLHEXSec_, GenLumiInfoProduct::XSec::error(), filterOnlyEffRun_, filterOnlyEffStat_, genFilterInfoToken_, genLumiInfoToken_, edm::LuminosityBlock::getByToken(), GenLumiInfoProduct::getHEPIDWTUP(), GenLumiInfoProduct::getProcessInfos(), hepidwtup_, hepMCFilterEffRun_, hepMCFilterEffStat_, hepMCFilterInfoToken_, hcalTTPDigis_cfi::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 &    )

overrideprivatevirtual

Definition at line 210 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_, findQualityFiles::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_);

}

const double GenXSecAnalyzer::final_xsec_error ( ) const

inline

Definition at line 41 of file GenXSecAnalyzer.h.

References analyze(), beginJob(), beginLuminosityBlock(), beginRun(), combine(), compute(), endJob(), endLuminosityBlock(), endRun(), GenLumiInfoProduct::XSec::error(), and xsec_.

{return xsec_.error();}

const double GenXSecAnalyzer::final_xsec_value ( ) const

inline

Definition at line 40 of file GenXSecAnalyzer.h.

References GenLumiInfoProduct::XSec::value(), and xsec_.

{return xsec_.value();}

Member Data Documentation

std::map<int, GenLumiInfoProduct::XSec> GenXSecAnalyzer::currentLumiBlockLHEXSec_

private

Definition at line 120 of file GenXSecAnalyzer.h.

Referenced by beginJob(), beginRun(), endLuminosityBlock(), endRun(), and GenXSecAnalyzer().

GenFilterInfo GenXSecAnalyzer::filterOnlyEffRun_

private

Definition at line 89 of file GenXSecAnalyzer.h.

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

GenFilterInfo GenXSecAnalyzer::filterOnlyEffStat_

private

Definition at line 95 of file GenXSecAnalyzer.h.

Referenced by endJob(), and endLuminosityBlock().

edm::EDGetTokenT<GenFilterInfo> GenXSecAnalyzer::genFilterInfoToken_

private

Definition at line 58 of file GenXSecAnalyzer.h.

Referenced by endLuminosityBlock(), and GenXSecAnalyzer().

edm::EDGetTokenT<GenLumiInfoProduct> GenXSecAnalyzer::genLumiInfoToken_

private

Definition at line 60 of file GenXSecAnalyzer.h.

Referenced by endLuminosityBlock(), and GenXSecAnalyzer().

int GenXSecAnalyzer::hepidwtup_

private

Definition at line 67 of file GenXSecAnalyzer.h.

Referenced by compute(), endJob(), endLuminosityBlock(), and endRun().

GenFilterInfo GenXSecAnalyzer::hepMCFilterEffRun_

private

Definition at line 92 of file GenXSecAnalyzer.h.

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

GenFilterInfo GenXSecAnalyzer::hepMCFilterEffStat_

private

Definition at line 98 of file GenXSecAnalyzer.h.

Referenced by endJob(), and endLuminosityBlock().

edm::EDGetTokenT<GenFilterInfo> GenXSecAnalyzer::hepMCFilterInfoToken_

private

Definition at line 59 of file GenXSecAnalyzer.h.

Referenced by endLuminosityBlock(), and GenXSecAnalyzer().

std::map<int, GenFilterInfo> GenXSecAnalyzer::jetMatchEffStat_

private

Definition at line 109 of file GenXSecAnalyzer.h.

Referenced by beginJob(), beginRun(), compute(), endJob(), endLuminosityBlock(), and GenXSecAnalyzer().

edm::EDGetTokenT<LHERunInfoProduct> GenXSecAnalyzer::lheRunInfoToken_

private

Definition at line 61 of file GenXSecAnalyzer.h.

Referenced by endRun(), and GenXSecAnalyzer().

int GenXSecAnalyzer::nMCs_

private

Definition at line 65 of file GenXSecAnalyzer.h.

Referenced by beginRun(), and endJob().

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

private

Definition at line 115 of file GenXSecAnalyzer.h.

Referenced by beginJob(), beginRun(), endLuminosityBlock(), and GenXSecAnalyzer().

GenLumiInfoProduct GenXSecAnalyzer::product_

private

Definition at line 85 of file GenXSecAnalyzer.h.

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

double GenXSecAnalyzer::thisRunWeight_

private

Definition at line 75 of file GenXSecAnalyzer.h.

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

double GenXSecAnalyzer::thisRunWeightPre_

private

Definition at line 71 of file GenXSecAnalyzer.h.

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

double GenXSecAnalyzer::totalWeight_

private

Definition at line 74 of file GenXSecAnalyzer.h.

Referenced by endRun().

double GenXSecAnalyzer::totalWeightPre_

private

Definition at line 70 of file GenXSecAnalyzer.h.

Referenced by endRun().

GenLumiInfoProduct::XSec GenXSecAnalyzer::xsec_

private

Definition at line 81 of file GenXSecAnalyzer.h.

Referenced by endJob(), endRun(), final_xsec_error(), and final_xsec_value().

std::vector<GenLumiInfoProduct::XSec> GenXSecAnalyzer::xsecAfterMatching_

private

Definition at line 107 of file GenXSecAnalyzer.h.

Referenced by beginJob(), beginRun(), compute(), endJob(), and GenXSecAnalyzer().

std::vector<GenLumiInfoProduct::XSec> GenXSecAnalyzer::xsecBeforeMatching_

private

Definition at line 105 of file GenXSecAnalyzer.h.

Referenced by beginJob(), beginRun(), compute(), endJob(), and GenXSecAnalyzer().

GenLumiInfoProduct::XSec GenXSecAnalyzer::xsecPreFilter_

private

Definition at line 78 of file GenXSecAnalyzer.h.

Referenced by endJob(), and endRun().