#include <RivetHarvesting.h>

Inheritance diagram for RivetHarvesting:

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

void	beginJob () override

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

void	endJob () override

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

	RivetHarvesting (const edm::ParameterSet &)

	~RivetHarvesting () override

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

	EDAnalyzer ()

SerialTaskQueue *	globalLuminosityBlocksQueue ()

SerialTaskQueue *	globalRunsQueue ()

ModuleDescription const &	moduleDescription () const

std::string	workerType () const

	~EDAnalyzer () 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

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)

virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
std::vector< YODA::Point2D >	getPoint2DValsErrs (std::string filename, std::string path, std::string name)

Private Attributes
Rivet::AnalysisHandler	_analysisHandler

std::vector< std::string >	_analysisNames

std::vector< double >	_crossSections

std::vector< std::string >	_fileNames

edm::InputTag	_hepmcCollection

bool	_isFirstEvent

std::vector< double >	_lumis

std::string	_outFileName

std::vector< double >	_sumOfWeights

Additional Inherited Members
Public Types inherited from edm::EDAnalyzer
typedef EDAnalyzer	ModuleType

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

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

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &)

static bool	wantsGlobalLuminosityBlocks ()

static bool	wantsGlobalRuns ()

static bool	wantsStreamLuminosityBlocks ()

static bool	wantsStreamRuns ()

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 19 of file RivetHarvesting.h.

Constructor & Destructor Documentation

RivetHarvesting::RivetHarvesting ( const edm::ParameterSet & pset )

Definition at line 26 of file RivetHarvesting.cc.

References _analysisHandler, _analysisNames, _crossSections, _fileNames, _lumis, _sumOfWeights, Exception, edm::ParameterSet::getParameter(), and mps_fire::i.

                                                             : 
 _analysisHandler(),
 _fileNames(pset.getParameter<std::vector<std::string> >("FilesToHarvest")),
 _sumOfWeights(pset.getParameter<std::vector<double> >("VSumOfWeights")),
 _crossSections(pset.getParameter<std::vector<double> >("VCrossSections")),
 _outFileName(pset.getParameter<std::string>("OutputFile")),
 _isFirstEvent(true),
 _hepmcCollection(pset.getParameter<edm::InputTag>("HepMCCollection")),
 _analysisNames(pset.getParameter<std::vector<std::string> >("AnalysisNames"))
 {
 
   if (_sumOfWeights.size() != _fileNames.size() ||
       _sumOfWeights.size() != _crossSections.size() ||
       _fileNames.size()    != _crossSections.size()){
     throw cms::Exception("RivetHarvesting") << "Mismatch in vector sizes: FilesToHarvest: " << _sumOfWeights.size() << ", VSumOfWeights: " << _sumOfWeights.size() << ", VCrossSections: " << _crossSections.size();  
   }    
         
   
   //get the analyses
   _analysisHandler.addAnalyses(_analysisNames);
 
   //go through the analyses and check those that need the cross section
   const std::set< AnaHandle, CmpAnaHandle > & analyses = _analysisHandler.analyses();
 
   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");
   for (iana = ibeg; iana != iend; ++iana){
     if ((*iana)->needsCrossSection())
       (*iana)->setCrossSection(xsection);
   }
 
   double totalSumOfWeights = _sumOfWeights[0];
   _lumis.push_back(_sumOfWeights[0]/_crossSections[0]); 
   for (unsigned int i = 1; i < _sumOfWeights.size(); ++i){
     _lumis.push_back(_sumOfWeights[i]/_crossSections[i]);
     totalSumOfWeights += _sumOfWeights[i]*_lumis[0]/_lumis[i];
   }
   _analysisHandler.setSumOfWeights(totalSumOfWeights);  
 
 
 }

RivetHarvesting::~RivetHarvesting ( )

override

Definition at line 71 of file RivetHarvesting.cc.

71 {

72 }

Member Function Documentation

void RivetHarvesting::analyze	(	const edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

override

Definition at line 89 of file RivetHarvesting.cc.

References _analysisHandler, _hepmcCollection, _isFirstEvent, edm::Event::getByLabel(), and edm::HepMCProduct::GetEvent().

                                                                                {
   if (!_isFirstEvent)
     return;
 
   //initialize the analysis handles, all histograms are booked
   //we need at least one event to get the handler initialized
   edm::Handle<HepMCProduct> evt;
   iEvent.getByLabel(_hepmcCollection, evt);
 
   // get HepMC GenEvent
   const HepMC::GenEvent *myGenEvent = evt->GetEvent();
   _analysisHandler.init(*myGenEvent);
   //gain access to the histogram factory and change the histograms
   /*
   AIDA::ITree & tree = _analysisHandler.tree();
   tree.ls(".", true);
 
   //from Analysis.hh (cls 18Feb2014)
   //const vector<AnalysisObjectPtr>& analysisObjects() const {
   //return _analysisobjects;
   //}
 
 
 
   for (std::vector<std::string>::const_iterator iAna = _analysisNames.begin(); iAna != _analysisNames.end(); ++iAna){
     std::vector<std::string> listOfNames = tree.listObjectNames("./"+(*iAna), true);
     std::vector<std::string>::const_iterator iNameBeg = listOfNames.begin();
     std::vector<std::string>::const_iterator iNameEnd = listOfNames.end();
     for (std::vector<std::string>::const_iterator iName = iNameBeg; iName != iNameEnd; ++iName ){
       AIDA::IManagedObject * iObj = tree.find(*iName);
       if (!iObj){
         std::cout << *iName << " not found; SKIPPING!" << std::endl;
         continue;
       } 
       
       std::cout << *iName << " FOUND!" << std::endl;
       vector<string>::const_iterator iFile;
       vector<string>::const_iterator iFileBeg = _fileNames.begin(); 
       vector<string>::const_iterator iFileEnd = _fileNames.end();
       AIDA::IHistogram1D* histo = dynamic_cast<AIDA::IHistogram1D*>(iObj); 
       AIDA::IProfile1D*   prof  = dynamic_cast<AIDA::IProfile1D*>(iObj);
       string tmpdir = "/tmpdir";
       tree.mkdir(tmpdir);
       unsigned int ifc = 0;
       for (iFile = iFileBeg; iFile != iFileEnd; ++iFile) {
         std::cout << "opening file " << *iFile << std::endl;
         string name = *iName;
         string tostrip = *iAna+'/';
         name.replace(name.find(tostrip),tostrip.length(),"");
         name.replace(name.find("/"),1,"");
         cout << name << endl;
         //vector<DPSXYPoint> original = getDPSXYValsErrs(*iFile, *iAna, name);
         vector<Point2D> original = getPoint2DValsErrs(*iFile, *iAna, name);
         if (histo){
           const string tmppath = tmpdir + "/" + name;
           cout << tmppath << endl;
           IHistogram1D* tmphisto = _analysisHandler.histogramFactory().createCopy(tmppath, *histo);
           tmphisto->reset();
           for (unsigned int i = 0; i < original.size(); ++i){
             tmphisto->fill(original[i].xval, original[i].yval);
           }
           tmphisto->scale(_lumis[ifc]);
           histo->add(*tmphisto);
           //iObj = new AIDA::IHistogram1D(*(_analysisHandler.histogramFactory().add(*iName, *histo, *tmphisto)));
           tree.rm(tmppath);
           //delete tmphisto;
         } else if (prof) {
           std::cout << *iName << "is a profile, doing nothing " << std::endl;
         } else {
           std::cout << *iName << " is neither a IHistogram1D neither a IProfile1D. Doing nothing with it." << std::endl;
         }
         ++ifc;
       }
       cout << iObj << endl;
     }
   }
 
   tree.ls(".", true);
   */
   _isFirstEvent = false;
 }

void RivetHarvesting::beginJob ( void )

overridevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 74 of file RivetHarvesting.cc.

References AlCaHLTBitMon_QueryRunRegistry::string.

                               {
   //set the environment, very ugly but rivet is monolithic when it comes to paths
   char * cmsswbase    = getenv("CMSSW_BASE");
   char * cmsswrelease = getenv("CMSSW_RELEASE_BASE");
   std::string rivetref, rivetinfo;
   rivetref = "RIVET_REF_PATH=" + string(cmsswbase) + "/src/GeneratorInterface/RivetInterface/data:" + string(cmsswrelease) + "/src/GeneratorInterface/RivetInterface/data";
   rivetinfo = "RIVET_INFO_PATH=" + string(cmsswbase) + "/src/GeneratorInterface/RivetInterface/data:" + string(cmsswrelease) + "/src/GeneratorInterface/RivetInterface/data";
   putenv(strdup(rivetref.c_str()));
   putenv(strdup(rivetinfo.c_str()));
 }

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

override

Definition at line 85 of file RivetHarvesting.cc.

                                                                             {
   return;
 }

void RivetHarvesting::endJob ( void )

overridevirtual

Reimplemented from edm::EDAnalyzer.

Definition at line 177 of file RivetHarvesting.cc.

References _analysisHandler, and _outFileName.

                             {
   _analysisHandler.finalize();   
   _analysisHandler.writeData(_outFileName);
 }

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

override

Definition at line 173 of file RivetHarvesting.cc.

                                                                           {
   return;
 }

vector< YODA::Point2D > RivetHarvesting::getPoint2DValsErrs	(	std::string	filename,
		std::string	path,
		std::string	name
	)

private

Check path to make sure that this is a reference histogram.

Definition at line 182 of file RivetHarvesting.cc.

References MessageLogger_cfi::cerr, DEFINE_FWK_MODULE, common_cff::doc, MillePedeFileConverter_cfg::e, cppFunctionSkipper::exception, Exception, corrVsCorr::filename, hiPixelPairStep_cff::points, EnergyCorrector::pt, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                                                           {
   
     // Open YODA XML file
     TiXmlDocument doc(filename);
     doc.LoadFile();
     if (doc.Error()) {
       string err = "Error in " + string(doc.Value());
       err += ": " + string(doc.ErrorDesc());
       cerr << err << endl;
       throw cms::Exception("RivetHarvesting") << "Cannot open " << filename;
     }
 
     // Return value, to be populated
     vector<Point2D> rtn;
     
     try {
       // Walk down tree to get to the <paper> element
       const TiXmlNode* yodaN = doc.FirstChild("yoda");
       if (!yodaN) throw cms::Exception("RivetHarvesting") << "Couldn't get <yoda> root element";
       for (const TiXmlNode* dpsN = yodaN->FirstChild("dataPointSet"); dpsN; dpsN = dpsN->NextSibling()) {
         const TiXmlElement* dpsE = dpsN->ToElement();
         const string plotname = dpsE->Attribute("name");
         const string plotpath = dpsE->Attribute("path");
         if (plotpath != path && plotname != name)
           continue;
         //if (plotpath.find("/REF") != 0) {
         //  cerr << "Skipping non-reference histogram " << plotname << endl;
         //  continue;
         //}
 
         vector<Point2D> points;
         for (const TiXmlNode* dpN = dpsN->FirstChild("dataPoint"); dpN; dpN = dpN->NextSibling()) {
           const TiXmlNode* xMeasN = dpN->FirstChild("measurement");
           const TiXmlNode* yMeasN = xMeasN->NextSibling();
           if (xMeasN && yMeasN)  {
             const TiXmlElement* xMeasE = xMeasN->ToElement();
             const TiXmlElement* yMeasE = yMeasN->ToElement();
             const string xcentreStr   = xMeasE->Attribute("value");
             const string xerrplusStr  = xMeasE->Attribute("errorPlus");
             const string xerrminusStr = xMeasE->Attribute("errorMinus");
             const string ycentreStr   = yMeasE->Attribute("value");
             const string yerrplusStr  = yMeasE->Attribute("errorPlus");
             const string yerrminusStr = yMeasE->Attribute("errorMinus");
             //if (!centreStr) throw Error("Couldn't get a valid bin centre");
             //if (!errplusStr) throw Error("Couldn't get a valid bin err+");
             //if (!errminusStr) throw Error("Couldn't get a valid bin err-");
             istringstream xssC(xcentreStr);
             istringstream xssP(xerrplusStr);
             istringstream xssM(xerrminusStr);
             istringstream yssC(ycentreStr);
             istringstream yssP(yerrplusStr);
             istringstream yssM(yerrminusStr);
             double xcentre, xerrplus, xerrminus, ycentre, yerrplus, yerrminus;
             xssC >> xcentre; xssP >> xerrplus; xssM >> xerrminus;
             yssC >> ycentre; yssP >> yerrplus; yssM >> yerrminus;
             //cout << "  " << centre << " + " << errplus << " - " << errminus << endl;
             Point2D pt(xcentre, xerrminus, xerrplus, ycentre, yerrminus, yerrplus);
             points.push_back(pt);
           } else {
             cerr << "Couldn't get <measurement> tag" << endl;
           }
         }
 
         return points;
       }
 
     }
     // Write out the error
     catch (std::exception& e) {
       cerr << e.what() << endl;
       throw;
     }
 
     throw cms::Exception("RivetHarvesting") << "could not find " << path << "/" << name << " in file " << filename;
   
     return rtn;
   
 }