Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix setting cross section in Rivet analaysis #15462

Merged
merged 1 commit into from Aug 22, 2016
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
Expand Up @@ -50,6 +50,7 @@ class RivetAnalyzer : public edm::EDAnalyzer
std::string _outFileName;
bool _doFinalize;
bool _produceDQM;
double _xsection;

DQMStore *dbe;
std::vector<MonitorElement *> _mes;
Expand Down
52 changes: 30 additions & 22 deletions GeneratorInterface/RivetInterface/plugins/RivetAnalyzer.cc
Expand Up @@ -19,7 +19,8 @@ _outFileName(pset.getParameter<std::string>("OutputFile")),
//decide whether to finlaize tthe plots or not.
//deciding not to finalize them can be useful for further harvesting of many jobs
_doFinalize(pset.getParameter<bool>("DoFinalize")),
_produceDQM(pset.getParameter<bool>("ProduceDQMOutput"))
_produceDQM(pset.getParameter<bool>("ProduceDQMOutput")),
_xsection(-1.)
{
//retrive the analysis name from paarmeter set
std::vector<std::string> analysisNames = pset.getParameter<std::vector<std::string> >("AnalysisNames");
Expand Down Expand Up @@ -47,11 +48,10 @@ _produceDQM(pset.getParameter<bool>("ProduceDQMOutput"))
std::set< AnaHandle, CmpAnaHandle >::const_iterator ibeg = analyses.begin();
std::set< AnaHandle, CmpAnaHandle >::const_iterator iend = analyses.end();
std::set< AnaHandle, CmpAnaHandle >::const_iterator iana;
double xsection = -1.;
xsection = pset.getParameter<double>("CrossSection");
_xsection = pset.getParameter<double>("CrossSection");
for (iana = ibeg; iana != iend; ++iana){
if ((*iana)->needsCrossSection())
(*iana)->setCrossSection(xsection);
(*iana)->setCrossSection(_xsection);
}
if (_produceDQM){
// book stuff needed for DQM
Expand Down Expand Up @@ -89,26 +89,34 @@ void RivetAnalyzer::analyze(const edm::Event& iEvent,const edm::EventSetup& iSet
// get HepMC GenEvent
const HepMC::GenEvent *myGenEvent = evt->GetEvent();

//if you want to use an external weight we have to clone the GenEvent and change the weight
if ( _useExternalWeight ){

//if you want to use an external weight or set the cross section we have to clone the GenEvent and change the weight
if ( _useExternalWeight || _xsection > 0 ){
HepMC::GenEvent * tmpGenEvtPtr = new HepMC::GenEvent( *(evt->GetEvent()) );
if (tmpGenEvtPtr->weights().size() == 0) {
throw cms::Exception("RivetAnalyzer") << "Original weight container has 0 size ";
}
if (tmpGenEvtPtr->weights().size() > 1) {
edm::LogWarning("RivetAnalyzer") << "Original event weight size is " << tmpGenEvtPtr->weights().size() << ". Will change only the first one ";
}

if (_xsection > 0){
HepMC::GenCrossSection xsec;
xsec.set_cross_section(_xsection);
tmpGenEvtPtr->set_cross_section(xsec);
}

if ( _useExternalWeight ){
if (tmpGenEvtPtr->weights().size() == 0) {
throw cms::Exception("RivetAnalyzer") << "Original weight container has 0 size ";
}
if (tmpGenEvtPtr->weights().size() > 1) {
edm::LogWarning("RivetAnalyzer") << "Original event weight size is " << tmpGenEvtPtr->weights().size() << ". Will change only the first one ";
}

if(!_useLHEweights){
edm::Handle<GenEventInfoProduct> genEventInfoProduct;
iEvent.getByToken(_genEventInfoCollection, genEventInfoProduct);
tmpGenEvtPtr->weights()[0] = genEventInfoProduct->weight();
}else{
edm::Handle<LHEEventProduct> lheEventHandle;
iEvent.getByToken(_LHECollection,lheEventHandle);
const LHEEventProduct::WGT& wgt = lheEventHandle->weights().at(_LHEweightNumber);
tmpGenEvtPtr->weights()[0] = wgt.wgt;
if(!_useLHEweights){
edm::Handle<GenEventInfoProduct> genEventInfoProduct;
iEvent.getByToken(_genEventInfoCollection, genEventInfoProduct);
tmpGenEvtPtr->weights()[0] = genEventInfoProduct->weight();
}else{
edm::Handle<LHEEventProduct> lheEventHandle;
iEvent.getByToken(_LHECollection,lheEventHandle);
const LHEEventProduct::WGT& wgt = lheEventHandle->weights().at(_LHEweightNumber);
tmpGenEvtPtr->weights()[0] = wgt.wgt;
}
}
myGenEvent = tmpGenEvtPtr;

Expand Down