|
|
[Rivet-svn] r2122 - in trunk: . pyext src src/Core
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Tue Dec 1 19:39:13 GMT 2009
Author: buckley
Date: Tue Dec 1 19:39:13 2009
New Revision: 2122
Log:
Moving core class implementation files into a subdir to allow for more selective rebuilds
Added:
trunk/src/Core/
trunk/src/Core/Analysis.cc
- copied unchanged from r2121, trunk/src/Analysis.cc
trunk/src/Core/AnalysisHandler.cc
- copied unchanged from r2120, trunk/src/AnalysisHandler.cc
trunk/src/Core/AnalysisInfo.cc
- copied unchanged from r2121, trunk/src/AnalysisInfo.cc
trunk/src/Core/AnalysisLoader.cc
- copied unchanged from r2120, trunk/src/AnalysisLoader.cc
trunk/src/Core/Cuts.cc
- copied unchanged from r2120, trunk/src/Cuts.cc
trunk/src/Core/Event.cc
- copied unchanged from r2120, trunk/src/Event.cc
trunk/src/Core/HistoHandler.cc
- copied unchanged from r2120, trunk/src/HistoHandler.cc
trunk/src/Core/Jet.cc
- copied unchanged from r2120, trunk/src/Jet.cc
trunk/src/Core/Makefile.am
- copied, changed from r2121, trunk/src/Makefile.am
trunk/src/Core/Particle.cc
- copied unchanged from r2120, trunk/src/Particle.cc
trunk/src/Core/Projection.cc
- copied unchanged from r2120, trunk/src/Projection.cc
trunk/src/Core/ProjectionApplier.cc
- copied unchanged from r2120, trunk/src/ProjectionApplier.cc
trunk/src/Core/ProjectionHandler.cc
- copied unchanged from r2120, trunk/src/ProjectionHandler.cc
trunk/src/Core/Run.cc
- copied unchanged from r2120, trunk/src/Run.cc
Deleted:
trunk/src/Analysis.cc
trunk/src/AnalysisHandler.cc
trunk/src/AnalysisInfo.cc
trunk/src/AnalysisLoader.cc
trunk/src/Cuts.cc
trunk/src/Event.cc
trunk/src/HistoHandler.cc
trunk/src/Jet.cc
trunk/src/Particle.cc
trunk/src/Projection.cc
trunk/src/ProjectionApplier.cc
trunk/src/ProjectionHandler.cc
trunk/src/Run.cc
Modified:
trunk/configure.ac
trunk/pyext/rivet.i
trunk/src/Makefile.am
Modified: trunk/configure.ac
==============================================================================
--- trunk/configure.ac Tue Dec 1 17:04:36 2009 (r2121)
+++ trunk/configure.ac Tue Dec 1 19:39:13 2009 (r2122)
@@ -2,7 +2,7 @@
AC_PREREQ(2.59)
AC_INIT([Rivet],[1.2.0a0],[rivet at projects.hepforge.org],[Rivet])
-AC_CONFIG_SRCDIR([src/Analysis.cc])
+AC_CONFIG_SRCDIR([src/Core/Analysis.cc])
AC_CONFIG_HEADERS([include/Rivet/Config/DummyConfig.hh include/Rivet/Config/RivetConfig.hh include/Rivet/Config/BuildOptions.hh])
AM_INIT_AUTOMAKE(dist-bzip2)
AC_CONFIG_MACRO_DIR([m4])
@@ -255,6 +255,7 @@
AC_EMPTY_SUBST
AC_CONFIG_FILES(include/Makefile include/Rivet/Makefile)
AC_CONFIG_FILES(src/Makefile)
+AC_CONFIG_FILES(src/Core/Makefile)
AC_CONFIG_FILES(src/Tools/Makefile src/Tools/yaml-cpp/Makefile)
AC_CONFIG_FILES(src/Projections/Makefile)
AC_CONFIG_FILES(src/Analyses/Makefile)
Modified: trunk/pyext/rivet.i
==============================================================================
--- trunk/pyext/rivet.i Tue Dec 1 17:04:36 2009 (r2121)
+++ trunk/pyext/rivet.i Tue Dec 1 19:39:13 2009 (r2122)
@@ -81,8 +81,7 @@
virtual std::string collider() const;
virtual std::string year() const;
virtual std::string status() const;
- // const std::pair<ParticleName,ParticleName>& beams() const { return _beams; }
- // const std::vector<std::pair<double,double> >& energies() const { return _energies; }
+ virtual const std::vector<std::pair<double,double> >& energies() const;
virtual std::vector<std::string> authors() const;
virtual std::vector<std::string> references() const;
virtual const BeamPair& beams() const;
Copied: trunk/src/Core/Analysis.cc (from r2121, trunk/src/Analysis.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/Analysis.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2121, trunk/src/Analysis.cc)
@@ -0,0 +1,514 @@
+// -*- C++ -*-
+#include "Rivet/Rivet.hh"
+#include "Rivet/RivetAIDA.hh"
+#include "Rivet/AnalysisHandler.hh"
+#include "Rivet/AnalysisInfo.hh"
+#include "Rivet/Analysis.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "LWH/AIManagedObject.h"
+using namespace AIDA;
+
+namespace Rivet {
+
+
+ Analysis::Analysis(const string& name)
+ : _crossSection(-1.0),
+ _gotCrossSection(false),
+ _needsCrossSection(false),
+ _analysishandler(NULL),
+ _madeHistoDir(false)
+ {
+ ProjectionApplier::_allowProjReg = false;
+ _defaultname = name;
+ AnalysisInfo* ai = AnalysisInfo::make(name);
+ assert(ai != 0);
+ _info.reset(ai);
+ assert(_info.get() != 0);
+ //setBeams(ANY, ANY);
+ }
+
+
+ Analysis::~Analysis()
+ { }
+
+
+ IAnalysisFactory& Analysis::analysisFactory() {
+ return handler().analysisFactory();
+ }
+
+
+ ITree& Analysis::tree() {
+ return handler().tree();
+ }
+
+
+ IHistogramFactory& Analysis::histogramFactory() {
+ return handler().histogramFactory();
+ }
+
+
+ IDataPointSetFactory& Analysis::datapointsetFactory() {
+ return handler().datapointsetFactory();
+ }
+
+
+ const string Analysis::histoDir() const {
+ string path = "/" + name();
+ if (handler().runName().length() > 0) {
+ path = "/" + handler().runName() + path;
+ }
+ while (find_first(path, "//")) {
+ replace_all(path, "//", "/");
+ }
+ return path;
+ }
+
+
+ const string Analysis::histoPath(const string& hname) const {
+ const string path = histoDir() + "/" + hname;
+ return path;
+ }
+
+
+ Log& Analysis::getLog() const {
+ string logname = "Rivet.Analysis." + name();
+ return Log::getLog(logname);
+ }
+
+
+ size_t Analysis::numEvents() const {
+ return handler().numEvents();
+ }
+
+
+ double Analysis::sumOfWeights() const {
+ return handler().sumOfWeights();
+ }
+
+
+ ////////////////////////////////////////////////////////////
+ // Metadata
+
+ const AnalysisInfo& Analysis::info() const {
+ assert(_info.get() != 0);
+ return *_info;
+ }
+
+ string Analysis::name() const {
+ if (_info && !_info->name().empty()) return _info->name();
+ return _defaultname;
+ }
+
+ string Analysis::spiresId() const {
+ if (!_info) return "NONE";
+ return _info->spiresId();
+ }
+
+ vector<string> Analysis::authors() const {
+ if (!_info) return std::vector<std::string>();
+ return _info->authors();
+ }
+
+ string Analysis::summary() const {
+ if (!_info) return "NONE";
+ return _info->summary();
+ }
+
+ string Analysis::description() const {
+ if (!_info) return "NONE";
+ return _info->description();
+ }
+
+ string Analysis::runInfo() const {
+ if (!_info) return "NONE";
+ return _info->runInfo();
+ }
+
+ const std::pair<ParticleName,ParticleName>& Analysis::beams() const {
+ return info().beams();
+ }
+
+ const std::vector<std::pair<double,double> >& Analysis::energies() const {
+ return info().energies();
+ }
+
+ string Analysis::experiment() const {
+ if (!_info) return "NONE";
+ return _info->experiment();
+ }
+
+ string Analysis::collider() const {
+ if (!_info) return "NONE";
+ return _info->collider();
+ }
+
+ string Analysis::year() const {
+ if (!_info) return "NONE";
+ return _info->year();
+ }
+
+ vector<string> Analysis::references() const {
+ if (!_info) return vector<string>();
+ return _info->references();
+ }
+
+ string Analysis::status() const {
+ if (!_info) return "UNVALIDATED";
+ return _info->status();
+ }
+
+ const BeamPair Analysis::requiredBeams() const {
+ return make_pdgid_pair(info().beams());
+ }
+
+ Analysis& Analysis::setBeams(const ParticleName& beam1, const ParticleName& beam2) {
+ assert(_info.get() != 0);
+ _info->_beams = make_pair(beam1, beam2);
+ return *this;
+ }
+
+ bool Analysis::isCompatible(const ParticleName& beam1, const ParticleName& beam2) const {
+ BeamPair beams(beam1, beam2);
+ return compatible(beams, requiredBeams());
+ /// @todo Need to also check internal consistency of the analysis'
+ /// beam requirements with those of the projections it uses.
+ }
+
+ bool Analysis::isCompatible(const BeamPair& beams) const {
+ return compatible(beams, requiredBeams());
+ /// @todo Need to also check internal consistency of the analysis'
+ /// beam requirements with those of the projections it uses.
+ }
+
+ Analysis& Analysis::setCrossSection(const double& xs) {
+ _crossSection = xs;
+ _gotCrossSection = true;
+ return *this;
+ }
+
+ bool Analysis::needsCrossSection() const {
+ return _needsCrossSection;
+ }
+
+ Analysis& Analysis::setNeedsCrossSection(bool needed) {
+ _needsCrossSection = needed;
+ return *this;
+ }
+
+ double Analysis::crossSection() const {
+ if (!_gotCrossSection || _crossSection < 0) {
+ string errMsg = "You did not set the cross section for the analysis " + name();
+ throw Error(errMsg);
+ }
+ return _crossSection;
+ }
+
+ double Analysis::crossSectionPerEvent() const {
+ const double sumW = sumOfWeights();
+ assert(sumW > 0);
+ return _crossSection / sumW;
+ }
+
+ AnalysisHandler& Analysis::handler() const {
+ return *_analysishandler;
+ }
+
+
+
+ ////////////////////////////////////////////////////////////
+ // Histogramming
+
+
+ void Analysis::_cacheBinEdges() {
+ _cacheXAxisData();
+ if (_histBinEdges.empty()) {
+ getLog() << Log::TRACE << "Getting histo bin edges from AIDA for paper " << name() << endl;
+ _histBinEdges = getBinEdges(_dpsData);
+ }
+ }
+
+
+ void Analysis::_cacheXAxisData() {
+ if (_dpsData.empty()) {
+ getLog() << Log::TRACE << "Getting DPS x-axis data from AIDA for paper " << name() << endl;
+ _dpsData = getDPSXValsErrs(name());
+ }
+ }
+
+
+ string Analysis::_makeAxisCode(const size_t datasetId, const size_t xAxisId, const size_t yAxisId) const {
+ stringstream axisCode;
+ axisCode << "d";
+ if (datasetId < 10) axisCode << 0;
+ axisCode << datasetId;
+ axisCode << "-x";
+ if (xAxisId < 10) axisCode << 0;
+ axisCode << xAxisId;
+ axisCode << "-y";
+ if (yAxisId < 10) axisCode << 0;
+ axisCode << yAxisId;
+ return axisCode.str();
+ }
+
+
+ IHistogram1D* Analysis::bookHistogram1D(const size_t datasetId, const size_t xAxisId,
+ const size_t yAxisId, const string& title,
+ const string& xtitle, const string& ytitle)
+ {
+ const string axisCode = _makeAxisCode(datasetId, xAxisId, yAxisId);
+ return bookHistogram1D(axisCode, title, xtitle, ytitle);
+ }
+
+
+ IHistogram1D* Analysis::bookHistogram1D(const string& hname, const string& title,
+ const string& xtitle, const string& ytitle)
+ {
+ // Get the bin edges (only read the AIDA file once)
+ _cacheBinEdges();
+ getLog() << Log::TRACE << "Using histo bin edges for " << name() << ":" << hname << endl;
+ const BinEdges edges = _histBinEdges.find(hname)->second;
+ _makeHistoDir();
+ const string path = histoPath(hname);
+ IHistogram1D* hist = histogramFactory().createHistogram1D(path, title, edges);
+ getLog() << Log::TRACE << "Made histogram " << hname << " for " << name() << endl;
+ hist->setXTitle(xtitle);
+ hist->setYTitle(ytitle);
+ return hist;
+ }
+
+
+ IHistogram1D* Analysis::bookHistogram1D(const string& hname,
+ const size_t nbins, const double lower, const double upper,
+ const string& title,
+ const string& xtitle, const string& ytitle) {
+ _makeHistoDir();
+ const string path = histoPath(hname);
+ IHistogram1D* hist = histogramFactory().createHistogram1D(path, title, nbins, lower, upper);
+ getLog() << Log::TRACE << "Made histogram " << hname << " for " << name() << endl;
+ hist->setXTitle(xtitle);
+ hist->setYTitle(ytitle);
+ return hist;
+ }
+
+
+ IHistogram1D* Analysis::bookHistogram1D(const string& hname,
+ const vector<double>& binedges,
+ const string& title,
+ const string& xtitle, const string& ytitle) {
+ _makeHistoDir();
+ const string path = histoPath(hname);
+ IHistogram1D* hist = histogramFactory().createHistogram1D(path, title, binedges);
+ getLog() << Log::TRACE << "Made histogram " << hname << " for " << name() << endl;
+ hist->setXTitle(xtitle);
+ hist->setYTitle(ytitle);
+ return hist;
+ }
+
+
+ /////////////////
+
+
+ IProfile1D* Analysis::bookProfile1D(const size_t datasetId, const size_t xAxisId,
+ const size_t yAxisId, const string& title,
+ const string& xtitle, const string& ytitle) {
+ const string axisCode = _makeAxisCode(datasetId, xAxisId, yAxisId);
+ return bookProfile1D(axisCode, title, xtitle, ytitle);
+ }
+
+
+ IProfile1D* Analysis::bookProfile1D(const string& hname, const string& title,
+ const string& xtitle, const string& ytitle)
+ {
+ // Get the bin edges (only read the AIDA file once)
+ _cacheBinEdges();
+ getLog() << Log::TRACE << "Using profile histo bin edges for " << name() << ":" << hname << endl;
+ const BinEdges edges = _histBinEdges.find(hname)->second;
+ if (getLog().isActive(Log::TRACE)) {
+ stringstream edges_ss;
+ foreach (const double be, edges) {
+ edges_ss << " " << be;
+ }
+ getLog() << Log::TRACE << "Edges:" << edges_ss.str() << endl;
+ }
+ _makeHistoDir();
+ const string path = histoPath(hname);
+ IProfile1D* prof = histogramFactory().createProfile1D(path, title, edges);
+ getLog() << Log::TRACE << "Made profile histogram " << hname << " for " << name() << endl;
+ prof->setXTitle(xtitle);
+ prof->setYTitle(ytitle);
+ return prof;
+ }
+
+
+ IProfile1D* Analysis::bookProfile1D(const string& hname,
+ const size_t nbins, const double lower, const double upper,
+ const string& title,
+ const string& xtitle, const string& ytitle) {
+ _makeHistoDir();
+ const string path = histoPath(hname);
+ IProfile1D* prof = histogramFactory().createProfile1D(path, title, nbins, lower, upper);
+ getLog() << Log::TRACE << "Made profile histogram " << hname << " for " << name() << endl;
+ prof->setXTitle(xtitle);
+ prof->setYTitle(ytitle);
+ return prof;
+ }
+
+
+ IProfile1D* Analysis::bookProfile1D(const string& hname,
+ const vector<double>& binedges,
+ const string& title,
+ const string& xtitle, const string& ytitle) {
+ _makeHistoDir();
+ const string path = histoPath(hname);
+ IProfile1D* prof = histogramFactory().createProfile1D(path, title, binedges);
+ getLog() << Log::TRACE << "Made profile histogram " << hname << " for " << name() << endl;
+ prof->setXTitle(xtitle);
+ prof->setYTitle(ytitle);
+ return prof;
+ }
+
+
+ ///////////////////
+
+
+
+ IDataPointSet* Analysis::bookDataPointSet(const string& hname, const string& title,
+ const string& xtitle, const string& ytitle) {
+ _makeHistoDir();
+ const string path = histoPath(hname);
+ IDataPointSet* dps = datapointsetFactory().create(path, title, 2);
+ getLog() << Log::TRACE << "Made data point set " << hname << " for " << name() << endl;
+ dps->setXTitle(xtitle);
+ dps->setYTitle(ytitle);
+ return dps;
+ }
+
+
+ IDataPointSet* Analysis::bookDataPointSet(const string& hname,
+ const size_t npts, const double lower, const double upper,
+ const string& title,
+ const string& xtitle, const string& ytitle) {
+ IDataPointSet* dps = bookDataPointSet(hname, title, xtitle, ytitle);
+ for (size_t pt = 0; pt < npts; ++pt) {
+ const double binwidth = (upper-lower)/npts;
+ const double bincentre = lower + (pt + 0.5) * binwidth;
+ dps->addPoint();
+ IMeasurement* meas = dps->point(pt)->coordinate(0);
+ meas->setValue(bincentre);
+ meas->setErrorPlus(binwidth/2.0);
+ meas->setErrorMinus(binwidth/2.0);
+ }
+ return dps;
+ }
+
+
+ IDataPointSet* Analysis::bookDataPointSet(const size_t datasetId, const size_t xAxisId,
+ const size_t yAxisId, const string& title,
+ const string& xtitle, const string& ytitle) {
+ // Get the bin edges (only read the AIDA file once)
+ _cacheXAxisData();
+ // Build the axis code
+ const string axisCode = _makeAxisCode(datasetId, xAxisId, yAxisId);
+ //const map<string, vector<DPSXPoint> > xpoints = getDPSXValsErrs(papername);
+ getLog() << Log::TRACE << "Using DPS x-positions for " << name() << ":" << axisCode << endl;
+ IDataPointSet* dps = bookDataPointSet(axisCode, title, xtitle, ytitle);
+ const vector<DPSXPoint> xpts = _dpsData.find(axisCode)->second;
+ for (size_t pt = 0; pt < xpts.size(); ++pt) {
+ dps->addPoint();
+ IMeasurement* meas = dps->point(pt)->coordinate(0);
+ meas->setValue(xpts[pt].val);
+ meas->setErrorPlus(xpts[pt].errplus);
+ meas->setErrorMinus(xpts[pt].errminus);
+ }
+ getLog() << Log::TRACE << "Made DPS " << axisCode << " for " << name() << endl;
+ return dps;
+ }
+
+
+ ////////////////////
+
+
+ void Analysis::_makeHistoDir() {
+ if (!_madeHistoDir) {
+ if (! name().empty()) {
+ // vector<string> dirs;
+ // split(dirs, histoDir(), "/");
+ // string pathpart;
+ // foreach (const string& d, dirs) {
+ //tree().mkdir();
+ //}
+ tree().mkdirs(histoDir());
+ }
+ _madeHistoDir = true;
+ }
+ }
+
+
+ void Analysis::normalize(AIDA::IHistogram1D*& histo, const double norm) {
+ if (!histo) {
+ getLog() << Log::ERROR << "Failed to normalise histo=NULL in analysis "
+ << name() << "(norm=" << norm << ")" << endl;
+ return;
+ }
+ const string hpath = tree().findPath(dynamic_cast<const AIDA::IManagedObject&>(*histo));
+ getLog() << Log::TRACE << "Normalizing histo " << hpath << " to " << norm << endl;
+
+ double oldintg = 0.0;
+ int nBins = histo->axis().bins();
+ for (int iBin = 0; iBin != nBins; ++iBin) {
+ // Leaving out factor of binWidth because AIDA's "height" already includes a width factor.
+ oldintg += histo->binHeight(iBin); // * histo->axis().binWidth(iBin);
+ }
+ if (oldintg == 0.0) {
+ getLog() << Log::WARN << "Histo " << hpath << " has null integral during normalisation" << endl;
+ return;
+ }
+
+ // Scale by the normalisation factor.
+ scale(histo, norm/oldintg);
+ }
+
+
+ void Analysis::scale(AIDA::IHistogram1D*& histo, const double scale) {
+ if (!histo) {
+ getLog() << Log::ERROR << "Failed to scale histo=NULL in analysis "
+ << name() << "(scale=" << scale << ")" << endl;
+ return;
+ }
+ const string hpath = tree().findPath(dynamic_cast<const AIDA::IManagedObject&>(*histo));
+ getLog() << Log::TRACE << "Scaling histo " << hpath << endl;
+
+ vector<double> x, y, ex, ey;
+ for (size_t i = 0, N = histo->axis().bins(); i < N; ++i) {
+ x.push_back(0.5 * (histo->axis().binLowerEdge(i) + histo->axis().binUpperEdge(i)));
+ ex.push_back(histo->axis().binWidth(i)*0.5);
+
+ // "Bin height" is a misnomer in the AIDA spec: width is neglected.
+ // We'd like to do this: y.push_back(histo->binHeight(i) * scale);
+ y.push_back(histo->binHeight(i)*scale/histo->axis().binWidth(i));
+
+ // "Bin error" is a misnomer in the AIDA spec: width is neglected.
+ // We'd like to do this: ey.push_back(histo->binError(i) * scale);
+ ey.push_back(histo->binError(i)*scale/(0.5*histo->axis().binWidth(i)));
+ }
+
+ string title = histo->title();
+ string xtitle = histo->xtitle();
+ string ytitle = histo->ytitle();
+
+ tree().mkdir("/tmpnormalize");
+ tree().mv(hpath, "/tmpnormalize");
+
+ AIDA::IDataPointSet* dps = datapointsetFactory().createXY(hpath, title, x, y, ex, ey);
+ dps->setXTitle(xtitle);
+ dps->setYTitle(ytitle);
+
+ tree().rm(tree().findPath(dynamic_cast<AIDA::IManagedObject&>(*histo)));
+ tree().rmdir("/tmpnormalize");
+
+ // Set histo pointer to null - it can no longer be used.
+ histo = 0;
+ }
+
+
+}
Copied: trunk/src/Core/AnalysisHandler.cc (from r2120, trunk/src/AnalysisHandler.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/AnalysisHandler.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/AnalysisHandler.cc)
@@ -0,0 +1,297 @@
+// -*- C++ -*-
+#include "Rivet/Rivet.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/AnalysisHandler.hh"
+#include "Rivet/RivetAIDA.hh"
+#include "LWH/AIManagedObject.h"
+
+using namespace AIDA;
+
+namespace Rivet {
+
+
+ AnalysisHandler::AnalysisHandler(string basefilename,
+ string runname, HistoFormat storetype)
+ : _runname(runname), _nRun(0), _iRun(0), _numEvents(0), _sumOfWeights(0.0) {
+ _theAnalysisFactory = createAnalysisFactory();
+ _setupFactories(basefilename, storetype);
+ }
+
+
+ AnalysisHandler::AnalysisHandler(IAnalysisFactory& afac, string basefilename,
+ string runname, HistoFormat storetype)
+ : _runname(runname), _nRun(0), _iRun(0), _numEvents(0), _sumOfWeights(0.0),
+ _theAnalysisFactory(&afac) {
+ _setupFactories(basefilename, storetype);
+ }
+
+
+ AnalysisHandler::~AnalysisHandler()
+ {
+ }
+
+
+ Log& AnalysisHandler::getLog() {
+ return Log::getLog("Rivet.Analysis.Handler");
+ }
+
+
+ void AnalysisHandler::init(int i, int N) {
+ getLog() << Log::DEBUG << "Initialising the analysis handler" << endl;
+ _nRun = N;
+ _iRun = i;
+ _numEvents = 0;
+ _sumOfWeights = 0.0;
+ foreach (Analysis* a, _analyses) {
+ getLog() << Log::DEBUG << "Initialising analysis: " << a->name() << endl;
+ // Allow projection registration in the init phase onwards
+ a->_allowProjReg = true;
+ a->init();
+ //getLog() << Log::DEBUG << "Checking consistency of analysis: " << a->name() << endl;
+ //a->checkConsistency();
+ getLog() << Log::DEBUG << "Done initialising analysis: " << a->name() << endl;
+ }
+ getLog() << Log::DEBUG << "Analysis handler initialised" << endl;
+ }
+
+
+ void AnalysisHandler::analyze(const GenEvent& ge) {
+ Event event(ge);
+ _numEvents++;
+ // Weights
+ const double weight = event.weight();
+ _sumOfWeights += weight;
+ getLog() << Log::DEBUG << "Event #" << _numEvents << " weight = " << weight << endl;
+ #ifdef HEPMC_HAS_CROSS_SECTION
+ if (ge.cross_section()) {
+ const double xs = ge.cross_section()->cross_section();
+ setCrossSection(xs);
+ }
+ #endif
+ foreach (Analysis* a, _analyses) {
+ getLog() << Log::DEBUG << "About to run analysis " << a->name() << endl;
+ a->analyze(event);
+ getLog() << Log::DEBUG << "Finished running analysis " << a->name() << endl;
+ }
+ }
+
+
+ void AnalysisHandler::finalize() {
+ getLog() << Log::INFO << "Finalising analyses" << endl;
+ foreach (Analysis* a, _analyses) {
+ a->finalize();
+ }
+
+ // Print out number of events processed
+ getLog() << Log::INFO << "Processed " << _numEvents << " event" << (_numEvents == 1 ? "" : "s") << endl;
+
+ // Change AIDA histos into data point sets
+ getLog() << Log::DEBUG << "Converting histograms to scatter plots" << endl;
+ assert(_theTree != 0);
+ _normalizeTree(tree());
+
+ // Delete analyses
+ getLog() << Log::DEBUG << "Deleting analyses" << endl;
+ foreach (Analysis* a, _analyses) {
+ delete a;
+ }
+ _analyses.clear();
+
+ // Print out MCnet boilerplate
+ cout << endl;
+ cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << endl;
+ cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:0809.4638" << endl;
+ }
+
+
+ AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) {
+ Analysis* analysis = AnalysisLoader::getAnalysis(analysisname);
+ if (analysis) { // < Check for null analysis.
+ getLog() << Log::DEBUG << "Adding analysis '" << analysisname << "'" << endl;
+ analysis->_analysishandler = this;
+ _analyses.insert(analysis);
+ }
+ return *this;
+ }
+
+
+ AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) {
+ Analysis* toremove = 0;
+ foreach (Analysis* a, _analyses) {
+ if (a->name() == analysisname) {
+ toremove = a;
+ break;
+ }
+ }
+ if (toremove) {
+ getLog() << Log::DEBUG << "Removing analysis '" << analysisname << "'" << endl;
+ _analyses.erase(toremove);
+ delete toremove;
+ }
+ return *this;
+ }
+
+
+ /// Remove beam-incompatible analyses from the run list.
+ AnalysisHandler& AnalysisHandler::removeIncompatibleAnalyses(const BeamPair& beams) {
+ vector<Analysis*> todelete;
+ foreach (Analysis* a, _analyses) {
+ if (! a->isCompatible(beams)) {
+ todelete.push_back(a);
+ }
+ }
+ foreach (Analysis* a, todelete) {
+ getLog() << Log::WARN << "Removing incompatible analysis '"
+ << a->name() << "'" << endl;
+ _analyses.erase(a);
+ delete a;
+ }
+ return *this;
+ }
+
+
+ void AnalysisHandler::_setupFactories(string basefilename, HistoFormat storetype) {
+ string filename(basefilename), storetypestr("");
+ if (storetype == AIDAML) {
+ if (!endsWith(filename, ".aida")) filename += ".aida";
+ storetypestr = "xml";
+ } else if (storetype == FLAT) {
+ if (!endsWith(filename, ".data")) filename += ".data";
+ storetypestr = "flat";
+ } else if (storetype == ROOT) {
+ if (!endsWith(filename, ".root")) filename += ".root";
+ storetypestr = "root";
+ }
+ _theTree = _theAnalysisFactory->createTreeFactory()->create(filename, storetypestr, false, true);
+ _theHistogramFactory = _theAnalysisFactory->createHistogramFactory(tree());
+ _theDataPointSetFactory = _theAnalysisFactory->createDataPointSetFactory(tree());
+ }
+
+
+ void AnalysisHandler::commitData() {
+ tree().commit();
+ }
+
+
+ void AnalysisHandler::_normalizeTree(ITree& tree) {
+ Log& log = getLog();
+ const vector<string> paths = tree.listObjectNames("/", true); // args set recursive listing
+ log << Log::TRACE << "Number of objects in AIDA tree = " << paths.size() << endl;
+ const string tmpdir = "/RivetNormalizeTmp";
+ tree.mkdir(tmpdir);
+ foreach (const string& path, paths) {
+
+ IManagedObject* hobj = tree.find(path);
+ if (hobj) {
+ IHistogram1D* histo = dynamic_cast<IHistogram1D*>(hobj);
+ IProfile1D* prof = dynamic_cast<IProfile1D*>(hobj);
+ // If it's a normal histo:
+ if (histo) {
+ log << Log::TRACE << "Converting histo " << path << " to DPS" << endl;
+ tree.mv(path, tmpdir);
+ const size_t lastslash = path.find_last_of("/");
+ const string basename = path.substr(lastslash+1, path.length() - (lastslash+1));
+ const string tmppath = tmpdir + "/" + basename;
+ IHistogram1D* tmphisto = dynamic_cast<IHistogram1D*>(tree.find(tmppath));
+ if (tmphisto) {
+ //getLog() << Log::TRACE << "Temp histo " << tmppath << " exists" << endl;
+ datapointsetFactory().create(path, *tmphisto);
+ }
+ tree.rm(tmppath);
+ }
+ // If it's a profile histo:
+ else if (prof) {
+ log << Log::TRACE << "Converting profile histo " << path << " to DPS" << endl;
+ tree.mv(path, tmpdir);
+ const size_t lastslash = path.find_last_of("/");
+ const string basename = path.substr(lastslash+1, path.length() - (lastslash+1));
+ const string tmppath = tmpdir + "/" + basename;
+ IProfile1D* tmpprof = dynamic_cast<IProfile1D*>(tree.find(tmppath));
+ if (tmpprof) {
+ //getLog() << Log::TRACE << "Temp profile histo " << tmppath << " exists" << endl;
+ datapointsetFactory().create(path, *tmpprof);
+ }
+ tree.rm(tmppath);
+ }
+
+ }
+
+ }
+ tree.rmdir(tmpdir);
+ }
+
+
+ string AnalysisHandler::runName() const { return _runname; }
+ size_t AnalysisHandler::numEvents() const { return _numEvents; }
+ double AnalysisHandler::sumOfWeights() const { return _sumOfWeights; }
+
+ void AnalysisHandler::setSumOfWeights(const double& sum) {
+ _sumOfWeights=sum;
+ }
+
+
+ std::vector<std::string> AnalysisHandler::analysisNames() {
+ std::vector<std::string> rtn;
+ foreach (Analysis* a, _analyses) {
+ rtn.push_back(a->name());
+ }
+ return rtn;
+ }
+
+
+ AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector<std::string>& analysisnames) {
+ foreach (const string& aname, analysisnames) {
+ //getLog() << Log::DEBUG << "Adding analysis '" << aname << "'" << endl;
+ addAnalysis(aname);
+ }
+ return *this;
+ }
+
+
+ AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector<std::string>& analysisnames) {
+ foreach (const string& aname, analysisnames) {
+ removeAnalysis(aname);
+ }
+ return *this;
+ }
+
+
+
+ AIDA::IAnalysisFactory& AnalysisHandler::analysisFactory() {
+ return *_theAnalysisFactory;
+ }
+
+
+ AIDA::ITree& AnalysisHandler::tree() {
+ return *_theTree;
+ }
+
+
+ AIDA::IHistogramFactory& AnalysisHandler::histogramFactory() {
+ return *_theHistogramFactory;
+ }
+
+
+ AIDA::IDataPointSetFactory& AnalysisHandler::datapointsetFactory() {
+ return *_theDataPointSetFactory;
+ }
+
+
+ bool AnalysisHandler::needCrossSection() const {
+ bool rtn = false;
+ foreach (const Analysis* a, _analyses) {
+ if (!rtn) rtn = a->needsCrossSection();
+ if (rtn) break;
+ }
+ return rtn;
+ }
+
+
+ AnalysisHandler& AnalysisHandler::setCrossSection(double xs) {
+ foreach (Analysis* a, _analyses) {
+ a->setCrossSection(xs);
+ }
+ return *this;
+ }
+
+}
Copied: trunk/src/Core/AnalysisInfo.cc (from r2121, trunk/src/AnalysisInfo.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/AnalysisInfo.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2121, trunk/src/AnalysisInfo.cc)
@@ -0,0 +1,173 @@
+#include "Rivet/Rivet.hh"
+#include "Rivet/AnalysisInfo.hh"
+#include "Rivet/Tools/Utils.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "yaml.h"
+#include <iostream>
+#include <fstream>
+#include <unistd.h>
+
+namespace Rivet {
+
+
+ /// Ideas:
+ /// * search RIVET_DATA_PATH etc. for <name>.info.yaml
+ /// * how to determine the name?
+ /// * only populate pointer on Analysis when requested
+ /// * use smart pointer: deletes automatically when Analysis
+ /// goes out of scope
+
+
+ /// Static factory method
+ AnalysisInfo* AnalysisInfo::make(const std::string& ananame) {
+ // Build the list of directories to search
+ vector<string> dirs;
+ char* env = 0;
+ // First try to use the Rivet data path variable
+ env = getenv("RIVET_DATA_PATH");
+ if (env) dirs += split(env);
+ // Then try to use the Rivet data install path
+ dirs += getRivetDataPath();
+ // And the current dir
+ dirs += ".";
+
+ bool found = false;
+ string datapath = "";
+ foreach (const string& d, dirs) {
+ if (d.empty()) continue;
+ /// @todo Use system-independent separator (e.g. Boost.File)
+ datapath = d + "/" + ananame + ".info";
+ Log::getLog("Rivet.AnalysisInfo")
+ << Log::TRACE << "Looking for analysis data file '" << datapath << "'" << endl;
+ if (access(datapath.c_str(), R_OK) == 0) {
+ found = true;
+ break;
+ }
+ }
+
+ // Returned AI, in semi-null state
+ AnalysisInfo* ai = new AnalysisInfo();
+ ai->_beams = make_pair(ANY,ANY);
+ ai->_name = ananame;
+
+ /// If no ana data file found, return null AI
+ if (!found) return ai;
+
+ // Read data from YAML document
+ Log::getLog("Rivet.AnalysisInfo")
+ << Log::DEBUG << "Reading analysis data from " << datapath << endl;
+ std::ifstream io(datapath.c_str());
+ YAML::Parser parser(io);
+ YAML::Node doc;
+ try {
+ parser.GetNextDocument(doc);
+ //cout << doc << endl;
+ } catch (const YAML::ParserException& ex) {
+ Log::getLog("Rivet.AnalysisInfo")
+ << Log::ERROR << "Parse error when reading analysis data from "
+ << datapath << endl;
+ return ai;
+ }
+
+ for (YAML::Iterator it = doc.begin(); it != doc.end(); ++it) {
+ string key;
+ it.first() >> key;
+ stringstream sec;
+ sec << it.second();
+ const string secstr = sec.str().substr(0, sec.str().length()-1);
+ Log::getLog("Rivet.AnalysisInfo")
+ << Log::TRACE << key << ": " << secstr << endl;
+ try {
+ if (key == "Name") {
+ it.second() >> ai->_name;
+ } else if (key == "Summary") {
+ it.second() >> ai->_summary;
+ } else if (key == "Experiment") {
+ it.second() >> ai->_experiment;
+ } else if (key == "Beams") {
+ const YAML::Node& beams = it.second();
+ vector<ParticleName> beampair;
+ for (YAML::Iterator b = beams.begin(); b != beams.end(); ++b) {
+ string bstr;
+ *b >> bstr;
+ ParticleName beamname = getParticleNameEnum(bstr);
+ beampair += beamname;
+ }
+ assert(beampair.size() == 2);
+ ai->_beams = make_pair<ParticleName,ParticleName>(beampair[0], beampair[1]);
+ // } else if (key == "NeedCrossSection") {
+ // // it.second() >> ai->_needsCrossSection;
+ } else if (key == "Energies") {
+ const YAML::Node& energies = it.second();
+ vector<pair<double,double> > beam_energy_pairs;
+ for (YAML::Iterator be = energies.begin(); be != energies.end(); ++be) {
+ if (be->GetType() == YAML::CT_SCALAR) {
+ // If beam energy is a scalar, then assume symmetric beams each with half that energy
+ double sqrts;
+ *be >> sqrts;
+ beam_energy_pairs += make_pair(sqrts/2.0, sqrts/2.0);
+ } else if (be->GetType() == YAML::CT_SEQUENCE) {
+ const YAML::Node& beamenergy_strs = be.second();
+ vector<double> beamenergies;
+ for (YAML::Iterator e = beamenergy_strs.begin(); e != beamenergy_strs.end(); ++e) {
+ double beamenergy;
+ *e >> beamenergy;
+ beamenergies += beamenergy;
+ }
+ assert(beamenergies.size() == 2);
+ beam_energy_pairs += make_pair(beamenergies[0], beamenergies[1]);
+ } else {
+ assert(0 && "Beam energies have to be a list of either numbers or pairs of numbers");
+ }
+ }
+ ai->_energies = beam_energy_pairs;
+ } else if (key == "Collider") {
+ it.second() >> ai->_collider;
+ } else if (key == "SpiresID") {
+ it.second() >> ai->_spiresId;
+ } else if (key == "Status") {
+ it.second() >> ai->_status;
+ } else if (key == "RunInfo") {
+ it.second() >> ai->_runInfo;
+ } else if (key == "Description") {
+ it.second() >> ai->_description;
+ } else if (key == "Year") {
+ it.second() >> ai->_year;
+ } else if (key == "Authors") {
+ const YAML::Node& authors = it.second();
+ for (YAML::Iterator a = authors.begin(); a != authors.end(); ++a) {
+ string astr;
+ *a >> astr;
+ ai->_authors += astr;
+ }
+ } else if (key == "References") {
+ const YAML::Node& refs = it.second();
+ for (YAML::Iterator r = refs.begin(); r != refs.end(); ++r) {
+ string rstr;
+ *r >> rstr;
+ ai->_references += rstr;
+ }
+ }
+ } catch (const YAML::RepresentationException& ex) {
+ Log::getLog("Rivet.Analysis")
+ << Log::WARN << "Type error when reading analysis data '"
+ << key << "' from " << datapath << endl;
+ }
+ }
+ Log::getLog("Rivet.AnalysisInfo") << Log::DEBUG << ai << endl;
+ return ai;
+ }
+
+
+ string toString(const AnalysisInfo& ai) {
+ stringstream ss;
+ ss << ai.name();
+ ss << " - " << ai.summary();
+ // ss << " - " << ai.beams();
+ // ss << " - " << ai.energies();
+ ss << " (" << ai.status() << ")";
+ return ss.str();
+ }
+
+
+}
Copied: trunk/src/Core/AnalysisLoader.cc (from r2120, trunk/src/AnalysisLoader.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/AnalysisLoader.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/AnalysisLoader.cc)
@@ -0,0 +1,110 @@
+// -*- C++ -*-
+#include "Rivet/AnalysisLoader.hh"
+#include "Rivet/Tools/Utils.hh"
+#include "Rivet/Tools/osdir.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Analysis.hh"
+#include <dlfcn.h>
+
+namespace Rivet {
+
+
+ // Initialise static ptr collection
+ AnalysisLoader::AnalysisBuilderMap AnalysisLoader::_ptrs;
+
+
+ vector<string> AnalysisLoader::analysisNames() {
+ _loadAnalysisPlugins();
+ vector<string> names;
+ foreach (const AnalysisBuilderMap::value_type& p, _ptrs) names += p.first;
+ return names;
+ }
+
+
+ Analysis* AnalysisLoader::getAnalysis(const string& analysisname) {
+ _loadAnalysisPlugins();
+ AnalysisBuilderMap::const_iterator ai = _ptrs.find(analysisname);
+ if (ai == _ptrs.end()) return 0;
+ return ai->second->mkAnalysis();
+ }
+
+
+ vector<Analysis*> AnalysisLoader::getAllAnalyses() {
+ _loadAnalysisPlugins();
+ vector<Analysis*> analyses;
+ foreach (const AnalysisBuilderMap::value_type& p, _ptrs) {
+ analyses += p.second->mkAnalysis();
+ }
+ return analyses;
+ }
+
+
+ void AnalysisLoader::_registerBuilder(const AnalysisBuilderBase* ab) {
+ if (!ab) return;
+ const string name = ab->name();
+ if (_ptrs.find(name) != _ptrs.end()) {
+ // Protect against overwriting analyses by throwing an error
+ /// @todo Tidy this up!
+ Log::getLog("Rivet.AnalysisLoader") << Log::ERROR << "Tried to register a second plugin analysis called '" << name << "'" << endl;
+ exit(1);
+ }
+ Log::getLog("Rivet.AnalysisLoader") << Log::TRACE << "Registering a plugin analysis called '" << name << "'" << endl;
+ _ptrs[name] = ab;
+ }
+
+
+ void AnalysisLoader::_loadAnalysisPlugins() {
+ // Only run once
+ if (!_ptrs.empty()) return;
+
+ // Build the list of directories to search
+ vector<string> dirs;
+ char* env = 0;
+ // First use the Rivet analysis path variable
+ env = getenv("RIVET_ANALYSIS_PATH");
+ if (env) dirs += split(env);
+ // Then the Rivet library install path
+ dirs += getLibPath();
+ // And then the user's (non-system) library path
+ env = getenv("LD_LIBRARY_PATH");
+ if (env) dirs += split(env);
+ // Use the current dir, and its libtool subdir, too.
+ dirs += ".";
+ dirs += ".libs";
+
+ // Find plugin module library files
+ const string libsuffix = ".so";
+ vector<string> pluginfiles;
+ foreach (const string& d, dirs) {
+ if (d.empty()) continue;
+ oslink::directory dir(d);
+ while (dir) {
+ string filename = dir.next();
+ // Require that name *starts* with 'Rivet' with new loader
+ if (filename.find("Rivet") != 0) continue;
+ size_t posn = filename.find(libsuffix);
+ if (posn == string::npos || posn != filename.length()-libsuffix.length()) continue;
+ /// @todo Make sure this is an abs path
+ /// @todo Sys-dependent path separator instead of "/"
+ const string path = d + "/" + filename;
+ // Ensure no duplicate paths
+ if (find(pluginfiles.begin(), pluginfiles.end(), path) == pluginfiles.end()) {
+ pluginfiles += path;
+ }
+ }
+ }
+
+ // Load the plugin files
+ Log::getLog("Rivet.AnalysisLoader") << Log::TRACE << "Candidate analysis plugin libs: " << pluginfiles << endl;
+ foreach (const string& pf, pluginfiles) {
+ Log::getLog("Rivet.AnalysisLoader") << Log::TRACE << "Trying to load plugin analyses from file " << pf << endl;
+ void* handle = dlopen(pf.c_str(), RTLD_LAZY);
+ if (!handle) {
+ Log::getLog("Rivet.AnalysisLoader") << Log::WARN << "Cannot open " << pf << ": " << dlerror() << endl;
+ continue;
+ }
+ }
+ }
+
+
+}
Copied: trunk/src/Core/Cuts.cc (from r2120, trunk/src/Cuts.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/Cuts.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/Cuts.cc)
@@ -0,0 +1,77 @@
+// -*- C++ -*-
+
+#include "Rivet/Rivet.hh"
+#include "Rivet/Cuts.hh"
+
+using std::ostringstream;
+
+namespace Rivet {
+
+
+ Cuts& Cuts::addCut(const string& quantity, const Comparison& comparison, const double value) {
+ // If this quantity doesn't yet have any associated cuts, make a
+ // default cut with effectively infinitely loose cuts.
+ if (_cuts.find(quantity) == _cuts.end()) {
+ _cuts[quantity] = BinaryCut();
+ }
+ // Combine cuts in the most restrictive way.
+ switch (comparison) {
+ case LESS_EQ:
+ if (value < _cuts[quantity].getLowerThan()) {
+ _cuts[quantity].setLowerThan(value);
+ }
+ break;
+ case MORE_EQ:
+ if (value > _cuts[quantity].getHigherThan()) {
+ _cuts[quantity].setHigherThan(value);
+ }
+ break;
+ case EQUAL:
+ _cuts[quantity].setLowerThan(value);
+ _cuts[quantity].setHigherThan(value);
+ break;
+ }
+
+ // Allow method chaining.
+ return *this;
+ }
+
+
+
+ bool Cuts::checkConsistency() const {
+ for (Cuts::const_iterator c = begin(); c != end(); ++c) {
+ if (c->second.getLowerThan() < c->second.getLowerThan()) {
+ ostringstream msg;
+ msg << "Constraints on " << c->first << " are incompatible: "
+ << ">=" << c->second.getHigherThan() << " AND "
+ << "<=" << c->second.getLowerThan();
+ throw Error(msg.str());
+ }
+ }
+ return true;
+ }
+
+
+
+ ostream& Cuts::print(ostream & os) const {
+ for (Cuts::const_iterator cut = begin(); cut != end(); ++cut) {
+ os << endl;
+ os << setw(12) << std::left << cut->first;
+ if (cut->second.getHigherThan() > -numeric_limits<double>::max()) {
+ os << setw(3) << ">=";
+ os << setw(10) << std::right << cut->second.getHigherThan();
+ } else {
+ os << setw(13) << "";
+ }
+ if (cut->second.getLowerThan() < numeric_limits<double>::max()) {
+ os << setw(3) << "<=";
+ os << setw(10) << std::right << cut->second.getLowerThan();
+ } else {
+ os << setw(13) << "";
+ }
+ }
+ return os;
+ }
+
+
+}
Copied: trunk/src/Core/Event.cc (from r2120, trunk/src/Event.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/Event.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/Event.cc)
@@ -0,0 +1,108 @@
+#include "Rivet/Event.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Projections/Beam.hh"
+#include "Rivet/BeamConstraint.hh"
+#include "HepMC/GenEvent.h"
+
+namespace Rivet {
+
+
+ // Convert the GenEvent to use conventional units (GeV, mm)
+ void _geNormUnits(GenEvent& ge) {
+ // Specify units if supported
+ #ifdef HEPMC_HAS_UNITS
+ ge.use_units(HepMC::Units::GEV, HepMC::Units::MM);
+ #endif
+ }
+
+
+ void _geRot180x(GenEvent& ge) {
+ for (HepMC::GenEvent::particle_iterator ip = ge.particles_begin(); ip != ge.particles_end(); ++ip) {
+ const HepMC::FourVector& mom = (*ip)->momentum();
+ (*ip)->set_momentum(HepMC::FourVector(mom.px(), -mom.py(), -mom.pz(), mom.e()));
+ }
+ for (HepMC::GenEvent::vertex_iterator iv = ge.vertices_begin(); iv != ge.vertices_end(); ++iv) {
+ const HepMC::FourVector& pos = (*iv)->position();
+ (*iv)->set_position(HepMC::FourVector(pos.x(), -pos.y(), -pos.z(), pos.t()));
+ }
+ }
+
+
+ // Convert the GenEvent to use conventional alignment
+ // (proton or electron on +ve z-axis?)
+ // For example, FHerwig only produces DIS events in the
+ // unconventional orientation and has to be corrected
+ void _geNormAlignment(GenEvent& ge) {
+ if (!ge.valid_beam_particles()) return;
+ typedef pair<HepMC::GenParticle*, HepMC::GenParticle*> GPPair;
+ GPPair bps = ge.beam_particles();
+ const BeamPair beamids = make_pdgid_pair(bps.first->pdg_id(), bps.second->pdg_id());
+ Log::getLog("Rivet.Event") << Log::TRACE << "Beam IDs: " << beamids << endl;
+ const HepMC::GenParticle* plusgp = 0;
+ bool rot = false;
+
+ // Rotate e+- p and ppbar to put p along +z
+ /// @todo e+ e- convention? B-factories different from LEP?
+ // if (compatible(beamids, make_pdgid_pair(ELECTRON, PROTON)) ||
+ // compatible(beamids, make_pdgid_pair(POSITRON, PROTON)) ||
+ // compatible(beamids, make_pdgid_pair(ANTIPROTON, PROTON)) ) {
+ // Log::getLog("Rivet.Event") << Log::TRACE << "May need to rotate event..." << endl;
+ if (bps.first->pdg_id() != PROTON || bps.second->pdg_id() != PROTON) {
+ if (bps.first->pdg_id() == PROTON) {
+ plusgp = bps.first;
+ } else if (bps.second->pdg_id() == PROTON) {
+ plusgp = bps.second;
+ }
+ if (plusgp && plusgp->momentum().pz() < 0) {
+ rot = true;
+ }
+ }
+
+ // Do the rotation
+ if (rot) {
+ if (Log::getLog("Rivet.Event").isActive(Log::TRACE)) {
+ Log::getLog("Rivet.Event") << Log::TRACE << "Rotating event" << endl;
+ Log::getLog("Rivet.Event") << Log::TRACE << "Before rotation: "
+ << bps.first->pdg_id() << "@pz=" << bps.first->momentum().pz()/GeV << ", "
+ << bps.second->pdg_id() << "@pz=" << bps.second->momentum().pz()/GeV << endl;
+ }
+ _geRot180x(ge);
+ }
+ }
+
+
+ Event::Event(const GenEvent& ge)
+ : _genEvent(ge), _weight(1.0)
+ {
+ // Set the weight if there is one, otherwise default to 1.0
+ if (!_genEvent.weights().empty()) {
+ _weight = ge.weights()[0];
+ }
+
+ // Use Rivet's preferred units if possible
+ _geNormUnits(_genEvent);
+
+ // Use the conventional alignment
+ _geNormAlignment(_genEvent);
+
+ // Debug printout to check that copying/magling has worked
+ //_genEvent.print();
+ }
+
+
+ Event::Event(const Event& e)
+ : _genEvent(e._genEvent),
+ _weight(e._weight)
+ {
+ //
+ }
+
+
+ Event& Event::operator=(const Event& e) {
+ _genEvent = e._genEvent;
+ _weight = e._weight;
+ return *this;
+ }
+
+
+}
Copied: trunk/src/Core/HistoHandler.cc (from r2120, trunk/src/HistoHandler.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/HistoHandler.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/HistoHandler.cc)
@@ -0,0 +1,96 @@
+// -*- C++ -*-
+#include "Rivet/Rivet.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/HistoHandler.hh"
+#include "Rivet/Analysis.hh"
+#include <algorithm>
+
+namespace Rivet {
+
+
+ // Initialize instance pointer to null.
+ HistoHandler* HistoHandler::_instance = 0;
+
+
+ HistoHandler* HistoHandler::create() {
+ if (!_instance) {
+ _instance = new HistoHandler();
+ Log::getLog("Rivet.HistoHandler")
+ << Log::TRACE << "Created new HistoHandler at "
+ << _instance << endl;
+ }
+ return _instance;
+ }
+
+
+ // Get a logger.
+ Log& HistoHandler::getLog() const {
+ return Log::getLog("Rivet.HistoHandler");
+ }
+
+
+ void HistoHandler::clear() {
+ _namedhistos.clear();
+ }
+
+
+ // Delete contained pointers.
+ HistoHandler::~HistoHandler() {
+ clear();
+ }
+
+
+ const AnalysisObject* HistoHandler::registerAnalysisObject(const Analysis& parent,
+ const AnalysisObject& ao,
+ const string& name) {
+ getLog() << Log::TRACE << "Trying to register"
+ << " analysis object " << &ao
+ << " for parent " << &parent << "(" << parent.name() << ")"
+ << " with name '" << name << "'" << endl;
+
+ // If this name is already registered for this analysis, throw a complaint
+ NamedHistosMap::const_iterator nhs = _namedhistos.find(&parent);
+ if (nhs != _namedhistos.end()) {
+ NamedHistos::const_iterator nh = nhs->second.find(name);
+ if (nh != nhs->second.end()) {
+ stringstream ss;
+ ss << "Histogram \"" << name
+ << "\" already exists for parent analysis " << &parent;
+ throw Error(ss.str());
+ }
+ }
+
+ _namedhistos[&parent][name] = &ao;
+ //return *(_namedhistos[&parent][name]);
+ return const_cast<AnalysisObject*>(_namedhistos[&parent][name]);
+ }
+
+
+
+ AnalysisObject* HistoHandler::_getAnalysisObject(const Analysis& parent,
+ const string& name) const {
+ getLog() << Log::TRACE << "Searching for child histo '"
+ << name << "' of " << &parent << endl;
+
+ NamedHistosMap::const_iterator nhs = _namedhistos.find(&parent);
+ if (nhs == _namedhistos.end()) {
+ stringstream ss;
+ ss << "Couldn't find any histograms for parent analysis " << &parent;
+ throw Error(ss.str());
+ }
+
+ NamedHistos::const_iterator nh = nhs->second.find(name);
+ if (nh == nhs->second.end()) {
+ stringstream ss;
+ ss << "Couldn't find histogram \"" << name
+ << "\" for parent analysis " << &parent;
+ throw Error(ss.str());
+ }
+
+ //return *(nh->second);
+ AnalysisObject* rtn = const_cast<AnalysisObject*>(nh->second);
+ return rtn;
+ }
+
+
+}
Copied: trunk/src/Core/Jet.cc (from r2120, trunk/src/Jet.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/Jet.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/Jet.cc)
@@ -0,0 +1,215 @@
+#include "Rivet/Jet.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Tools/ParticleIdUtils.hh"
+#include "Rivet/ParticleName.hh"
+
+namespace Rivet {
+
+
+ Jet::Jet()
+ : ParticleBase()
+ {
+ clear();
+ }
+
+
+ Jet& Jet::setParticles(const vector<FourMomentum>& particles) {
+ _particles = particles;
+ _resetCaches();
+ return *this;
+ }
+
+
+ Jet& Jet::addParticle(const FourMomentum& particle) {
+ _particles.push_back(particle);
+ _resetCaches();
+ return *this;
+ }
+
+
+ Jet& Jet::addParticle(const Particle& particle) {
+ _fullParticles.push_back(particle);
+ _particles.push_back(particle.momentum());
+ _resetCaches();
+ return *this;
+ }
+
+
+ bool Jet::containsParticle(const Particle& particle) const {
+ const int barcode = particle.genParticle().barcode();
+ foreach (const Particle& p, particles()) {
+ if (p.genParticle().barcode() == barcode) return true;
+ }
+ return false;
+ }
+
+
+ bool Jet::containsParticleId(PdgId pid) const {
+ foreach (const Particle& p, particles()) {
+ if (p.pdgId() == pid) return true;
+ }
+ return false;
+ }
+
+
+ bool Jet::containsParticleId(const vector<PdgId>& pids) const {
+ foreach (const Particle& p, particles()) {
+ foreach (PdgId pid, pids) {
+ if (p.pdgId() == pid) return true;
+ }
+ }
+ return false;
+ }
+
+
+ /// @todo Jet::containsMatch(Matcher m) { ... if m(pid) return true; ... }
+
+
+ double Jet::totalEnergy() const {
+ return momentum().E();
+ }
+
+
+ double Jet::neutralEnergy() const {
+ double e_neutral = 0.0;
+ foreach (const Particle& p, particles()) {
+ const PdgId pid = p.pdgId();
+ if (PID::threeCharge(pid) == 0) {
+ e_neutral += p.momentum().E();
+ }
+ }
+ return e_neutral;
+ }
+
+
+ double Jet::hadronicEnergy() const {
+ double e_hadr = 0.0;
+ foreach (const Particle& p, particles()) {
+ const PdgId pid = p.pdgId();
+ if (PID::isHadron(pid)) {
+ e_hadr += p.momentum().E();
+ }
+ }
+ return e_hadr;
+ }
+
+
+ bool Jet::containsCharm() const {
+ foreach (const Particle& p, particles()) {
+ if (abs(p.pdgId()) == CQUARK) return true;
+ HepMC::GenVertex* gv = p.genParticle().production_vertex();
+ if (gv) {
+ foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) {
+ const PdgId pid = pi->pdg_id();
+ if (PID::isHadron(pid) && PID::hasCharm(pid)) return true;
+ }
+ }
+ }
+ return false;
+ }
+
+
+ bool Jet::containsBottom() const {
+ foreach (const Particle& p, particles()) {
+ if (abs(p.pdgId()) == BQUARK) return true;
+ HepMC::GenVertex* gv = p.genParticle().production_vertex();
+ if (gv) {
+ foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) {
+ const PdgId pid = pi->pdg_id();
+ if (PID::isHadron(pid) && PID::hasBottom(pid)) return true;
+ }
+ }
+ }
+ return false;
+ }
+
+
+ Jet& Jet::clear() {
+ _particles.clear();
+ _fullParticles.clear();
+ _resetCaches();
+ return *this;
+ }
+
+
+ double Jet::ptWeightedEta() const {
+ _calcPtAvgs();
+ assert(_okPtWeightedEta);
+ return _ptWeightedEta;
+ }
+
+
+ double Jet::ptWeightedPhi() const {
+ _calcPtAvgs();
+ assert(_okPtWeightedPhi);
+ return _ptWeightedPhi;
+ }
+
+
+ double Jet::eta() const {
+ return momentum().eta();
+
+ }
+
+
+ double Jet::phi() const {
+ return momentum().phi();
+ }
+
+
+ const FourMomentum& Jet::momentum() const {
+ _calcMomVector();
+ return _momentum;
+ }
+
+
+ // FourMomentum& Jet::momentum() {
+ // _calcMomVector();
+ // return _momentum;
+ // }
+
+
+ double Jet::ptSum() const {
+ return momentum().pT();
+ }
+
+
+ double Jet::EtSum() const {
+ return momentum().Et();
+ }
+
+
+ void Jet::_resetCaches() const {
+ _okPtWeightedPhi = false;
+ _okPtWeightedEta = false;
+ _okMomentum = false;
+ }
+
+
+ void Jet::_calcMomVector() const {
+ if (!_okMomentum) {
+ _momentum = accumulate(begin(), end(), FourMomentum());
+ _okMomentum = true;
+ }
+ }
+
+
+ void Jet::_calcPtAvgs() const {
+ if (!_okPtWeightedEta || !_okPtWeightedPhi) {
+ double ptwetasum(0.0), ptwdphisum(0.0), ptsum(0.0);
+ double phi0 = phi();
+ foreach (const FourMomentum& p, momenta()) {
+ double pt = p.pT();
+ ptsum += pt;
+ ptwetasum += pt * p.pseudorapidity();
+ ptwdphisum += pt * mapAngleMPiToPi(phi0 - p.azimuthalAngle());
+ }
+ _ptWeightedEta = ptwetasum/ptsum;
+ _okPtWeightedEta = true;
+ _ptWeightedPhi = mapAngleMPiToPi(phi0 + ptwdphisum/ptsum);
+ _okPtWeightedPhi = true;
+ }
+ }
+
+
+}
Copied and modified: trunk/src/Core/Makefile.am (from r2121, trunk/src/Makefile.am)
==============================================================================
--- trunk/src/Makefile.am Tue Dec 1 17:04:36 2009 (r2121, copy source)
+++ trunk/src/Core/Makefile.am Tue Dec 1 19:39:13 2009 (r2122)
@@ -1,24 +1,10 @@
-SUBDIRS = . Projections Tools Test
-if ENABLE_ANALYSES
-SUBDIRS += Analyses
-endif
+noinst_LTLIBRARIES = libRivetCore.la
-lib_LTLIBRARIES = libRivet.la
-
-libRivet_la_SOURCES = \
+libRivetCore_la_SOURCES = \
Event.cc Jet.cc Particle.cc \
ProjectionApplier.cc Projection.cc \
Analysis.cc AnalysisLoader.cc AnalysisInfo.cc \
AnalysisHandler.cc Run.cc ProjectionHandler.cc HistoHandler.cc
-libRivet_la_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/src/Tools/yaml-cpp $(CPPFLAGS)
-
-libRivet_la_LDFLAGS = -L$(prefix)/lib \
- -L$(GSLLIBPATH) -L$(GSLCBLASLIBPATH) -export-dynamic $(VERSIONINFOFLAGS)
-
-FASTJETLIBADD = $(FASTJETCONFIGLIBADD)
-
-libRivet_la_LIBADD = \
- Projections/libRivetProjections.la \
- Tools/libRivetTools.la -ldl -lm \
- $(GSL_LDFLAGS) $(FASTJETLIBADD)
+libRivetCore_la_CPPFLAGS = \
+ $(AM_CPPFLAGS) -I$(top_srcdir)/src/Tools/yaml-cpp $(CPPFLAGS)
Copied: trunk/src/Core/Particle.cc (from r2120, trunk/src/Particle.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/Particle.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/Particle.cc)
@@ -0,0 +1,16 @@
+#include "Rivet/Particle.hh"
+
+namespace Rivet {
+
+
+ bool Particle::hasAncestor(PdgId pdg_id) const {
+ GenVertex* prodVtx = genParticle().production_vertex();
+ if (prodVtx == 0) return false;
+ foreach (const GenParticle* ancestor, particles(prodVtx, HepMC::ancestors)) {
+ if (ancestor->pdg_id() == pdg_id) return true;
+ }
+ return false;
+ }
+
+
+}
Copied: trunk/src/Core/Projection.cc (from r2120, trunk/src/Projection.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/Projection.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/Projection.cc)
@@ -0,0 +1,63 @@
+// -*- C++ -*-
+#include "Rivet/Projection.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Event.hh"
+#include "Rivet/Cmp.hh"
+
+namespace Rivet {
+
+
+ Projection::Projection()
+ : _name("BaseProjection")
+ {
+ addBeamPair(ANY, ANY);
+ getLog() << Log::TRACE << "Creating " << name() << " at " << this << endl;
+ }
+
+
+ Projection:: ~Projection() {
+ getLog() << Log::TRACE << "Destroying " << name() << " at " << this << endl;
+ }
+
+
+ // int Projection::compare(const Projection& p) const {
+ // return mkNamedPCmp(p, "FS");
+ // }
+
+
+ bool Projection::before(const Projection& p) const {
+ const std::type_info& thisid = typeid(*this);
+ const std::type_info& otherid = typeid(p);
+ if (thisid == otherid) {
+ return compare(p) < 0;
+ } else {
+ return thisid.before(otherid);
+ }
+ }
+
+
+ const set<BeamPair> Projection::beamPairs() const {
+ set<BeamPair> ret = _beamPairs;
+ set<ConstProjectionPtr> projs = getProjections();
+ for (set<ConstProjectionPtr>::const_iterator ip = projs.begin(); ip != projs.end(); ++ip) {
+ ConstProjectionPtr p = *ip;
+ getLog() << Log::TRACE << "Proj addr = " << p << endl;
+ if (p) ret = intersection(ret, p->beamPairs());
+ }
+ return ret;
+ }
+
+
+ Cmp<Projection> Projection::mkNamedPCmp(const Projection& otherparent,
+ const string& pname) const {
+ return pcmp(*this, otherparent, pname);
+ }
+
+
+ Cmp<Projection> Projection::mkPCmp(const Projection& otherparent,
+ const string& pname) const {
+ return pcmp(*this, otherparent, pname);
+ }
+
+
+}
Copied: trunk/src/Core/ProjectionApplier.cc (from r2120, trunk/src/ProjectionApplier.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/ProjectionApplier.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/ProjectionApplier.cc)
@@ -0,0 +1,45 @@
+// -*- C++ -*-
+#include "Rivet/ProjectionApplier.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/Event.hh"
+
+namespace Rivet {
+
+
+ // NB. Allow proj registration in constructor by default -- explicitly disable for Analysis
+ ProjectionApplier::ProjectionApplier()
+ : _allowProjReg(true),
+ _projhandler(ProjectionHandler::create())
+ { }
+
+
+ ProjectionApplier::~ProjectionApplier() {
+ getProjHandler().removeProjectionApplier(*this);
+ }
+
+
+ const Projection& ProjectionApplier::_applyProjection(const Event& evt,
+ const string& name) const {
+ return evt.applyProjection(getProjection(name));
+ }
+
+
+ const Projection& ProjectionApplier::_applyProjection(const Event& evt,
+ const Projection& proj) const {
+ return evt.applyProjection(proj);
+ }
+
+
+ const Projection& ProjectionApplier::_addProjection(const Projection& proj,
+ const std::string& name) {
+ if (!_allowProjReg) {
+ getLog() << Log::ERROR << "Trying to register projection '"
+ << proj.name() << "' before init phase in '" << this->name() << "'." << endl;
+ exit(2);
+ }
+ const Projection& reg = getProjHandler().registerProjection(*this, proj, name);
+ return reg;
+ }
+
+
+}
Copied: trunk/src/Core/ProjectionHandler.cc (from r2120, trunk/src/ProjectionHandler.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/ProjectionHandler.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/ProjectionHandler.cc)
@@ -0,0 +1,284 @@
+// -*- C++ -*-
+#include "Rivet/Rivet.hh"
+#include "Rivet/Tools/Logging.hh"
+#include "Rivet/ProjectionHandler.hh"
+#include "Rivet/Cmp.hh"
+#include <algorithm>
+
+namespace Rivet {
+
+
+ // Initialize instance pointer to null.
+ ProjectionHandler* ProjectionHandler::_instance = 0;
+
+
+ //int Log_TRACE = Log::INFO;
+
+
+ ProjectionHandler* ProjectionHandler::create() {
+ if (!_instance) {
+ _instance = new ProjectionHandler();
+ Log::getLog("Rivet.ProjectionHandler")
+ << Log::TRACE << "Created new ProjectionHandler at " << _instance << endl;
+ }
+ return _instance;
+ }
+
+
+
+ // Get a logger.
+ Log& ProjectionHandler::getLog() const {
+ return Log::getLog("Rivet.ProjectionHandler");
+ }
+
+
+
+ void ProjectionHandler::clear() {
+ _projs.clear();
+ _namedprojs.clear();
+ }
+
+
+
+ // Delete contained pointers.
+ ProjectionHandler::~ProjectionHandler() {
+ clear();
+ }
+
+
+
+ // Take a Projection, compare it to the others on record, and return (by
+ // reference) an equivalent Projection which is guaranteed to be the
+ // (persistent) version that will be applied to an event.
+ const Projection& ProjectionHandler::registerProjection(const ProjectionApplier& parent,
+ const Projection& proj,
+ const string& name)
+ {
+ getLog() << Log::TRACE << "Trying to register"
+ << " projection " << &proj << " (" << proj.name() << ")"
+ << " for parent " << &parent << " (" << parent.name() << ")"
+ << " with name '" << name << "'" << endl;
+
+ // Check for duplicate use of "name" on "parent"
+ const bool dupOk = _checkDuplicate(parent, proj, name);
+ if (!dupOk) exit(1);
+
+ // Choose which version of the projection to register with this parent and name
+ const Projection* p = _getEquiv(proj);
+ if (p == 0) { // a null pointer is a non-match
+ // If there is no equivalent projection, clone proj and use the clone for registering
+ p = _clone(parent, proj);
+ }
+
+ // Do the registering
+ p = _register(parent, *p, name);
+
+ // Return registered proj
+ return *p;
+ }
+
+
+
+
+ // Attach and retrieve a projection as a pointer.
+ const Projection* ProjectionHandler::registerProjection(const ProjectionApplier& parent,
+ const Projection* proj,
+ const string& name) {
+ if (proj == 0) return 0;
+ const Projection& p = registerProjection(parent, *proj, name);
+ return &p;
+ }
+
+
+
+ // Clone neatly
+ const Projection* ProjectionHandler::_clone(const ProjectionApplier& parent,
+ const Projection& proj)
+ {
+ // Clone a new copy of the passed projection on the heap
+ getLog() << Log::TRACE << "Cloning projection " << proj.name() << " from " << &proj << endl;
+ const Projection* newproj = proj.clone();
+ getLog() << Log::TRACE << "Cloned projection " << proj.name() << " at " << newproj << endl;
+
+ // Copy all the child ProjHandles when cloning, since otherwise links to "stack parents"
+ // will be generated by their children, without any connection to the cloned parent
+ if (&proj != newproj) {
+ NamedProjsMap::const_iterator nps = _namedprojs.find(&proj);
+ if (nps != _namedprojs.end()) {
+ getLog() << Log::TRACE << "Cloning registered projections list: "
+ << &proj << " -> " << newproj << endl;
+ _namedprojs[newproj] = nps->second;
+ }
+ }
+
+ return newproj;
+ }
+
+
+
+ // Take a Projection, compare it to the others on record, and
+ // return (by reference) an equivalent Projection which is guaranteed to be
+ // the version that will be applied to an event.
+ const Projection* ProjectionHandler::_register(const ProjectionApplier& parent,
+ const Projection& proj,
+ const string& name)
+ {
+ ProjHandle ph(&proj);
+ getLog() << Log::TRACE << "Registering new projection at " << ph << endl;
+
+ // Add the passed Projection to _projs
+ _projs.insert(ph);
+
+ // Add the ProjApplier* => name location to the associative container
+ _namedprojs[&parent][name] = ph;
+
+ return ph.get();
+ }
+
+
+
+
+ // Try to find a equivalent projection in the system
+ const Projection* ProjectionHandler::_getEquiv(const Projection& proj) const
+ {
+ // Get class type using RTTI
+ const std::type_info& newtype = typeid(proj);
+ getLog() << Log::TRACE << "RTTI type of " << &proj << " is " << newtype.name() << endl;
+
+ // Compare to ALL projections via _projs collection
+ getLog() << Log::TRACE << "Comparing " << &proj
+ << " with " << _projs.size()
+ << " registered projection" << (_projs.size() == 1 ? "" : "s") << endl;
+ foreach (const ProjHandle& ph, _projs) {
+ // Make sure the concrete types match, using RTTI.
+ const std::type_info& regtype = typeid(*ph);
+ getLog() << Log::TRACE << "RTTI type comparison with " << ph << ": "
+ << newtype.name() << " vs. " << regtype.name() << endl;
+ if (newtype != regtype) continue;
+ getLog() << Log::TRACE << "RTTI type matches with " << ph << endl;
+
+ // Test for semantic match
+ if (pcmp(*ph, proj) != EQUIVALENT) {
+ getLog() << Log::TRACE << "Projections at "
+ << &proj << " and " << ph << " are not equivalent" << endl;
+ } else {
+ getLog() << Log::TRACE << "MATCH! Projections at "
+ << &proj << " and " << ph << " are equivalent" << endl;
+ return ph.get();
+ }
+ }
+
+ // If no match, just return a null pointer
+ return 0;
+ }
+
+
+
+ string ProjectionHandler::_getStatus() const {
+ ostringstream msg;
+ msg << "Current projection hierarchy:" << endl;
+ foreach (const NamedProjsMap::value_type& nps, _namedprojs) {
+ //const string parentname = nps.first->name();
+ msg << nps.first << endl; //"(" << parentname << ")" << endl;
+ foreach (const NamedProjs::value_type& np, nps.second) {
+ msg << " " << np.second << " (" << np.second->name()
+ << ", locally called '" << np.first << "')" << endl;
+ }
+ msg << endl;
+ }
+ return msg.str();
+ }
+
+
+
+ // Check that the same parent hasn't already used this name for something else
+ bool ProjectionHandler::_checkDuplicate(const ProjectionApplier& parent,
+ const Projection& proj,
+ const string& name) const
+ {
+ NamedProjsMap::const_iterator ipnps = _namedprojs.find(&parent);
+ if (ipnps != _namedprojs.end()) {
+ const NamedProjs pnps = ipnps->second;
+ const NamedProjs::const_iterator ipph = pnps.find(name);
+ if (ipph != pnps.end()) {
+ const ProjHandle pph = ipph->second;
+ getLog() << Log::ERROR << "Projection clash! "
+ << parent.name() << " (" << &parent << ") "
+ << "is trying to overwrite its registered '" << name << "' "
+ << "projection (" << pph << "="
+ << pph->name() << ") with a non-equivalent projection "
+ << "(" << &proj << "=" << proj.name() << ")" << endl;
+ getLog() << Log::ERROR << _getStatus();
+ return false;
+ }
+ }
+ return true;
+ }
+
+
+
+
+ void ProjectionHandler::removeProjectionApplier(ProjectionApplier& parent) {
+ NamedProjsMap::iterator npi = _namedprojs.find(&parent);
+ if (npi != _namedprojs.end()) _namedprojs.erase(npi);
+ //
+ const Projection* parentprojptr = dynamic_cast<Projection*>(&parent);
+ if (parentprojptr) {
+ ProjHandle parentph(parentprojptr);
+ ProjHandles::iterator pi = find(_projs.begin(), _projs.end(), parentph);
+ if (pi != _projs.end()) _projs.erase(pi);
+ }
+ }
+
+
+
+
+ set<const Projection*> ProjectionHandler::getChildProjections(const ProjectionApplier& parent,
+ ProjDepth depth) const
+ {
+ set<const Projection*> toplevel;
+ NamedProjs nps = _namedprojs.find(&parent)->second;
+ foreach (NamedProjs::value_type& np, nps) {
+ toplevel.insert(np.second.get());
+ }
+ if (depth == SHALLOW) {
+ // Only return the projections directly contained within the top level
+ return toplevel;
+ } else {
+ // Return recursively built projection list
+ set<const Projection*> alllevels = toplevel;
+ foreach (const Projection* p, toplevel) {
+ set<const Projection*> allsublevels = getChildProjections(*p, DEEP);
+ alllevels.insert(allsublevels.begin(), allsublevels.end());
+ }
+ return alllevels;
+ }
+ }
+
+
+
+
+ const Projection& ProjectionHandler::getProjection(const ProjectionApplier& parent,
+ const string& name) const {
+ getLog() << Log::TRACE << "Searching for child projection '"
+ << name << "' of " << &parent << endl;
+ NamedProjsMap::const_iterator nps = _namedprojs.find(&parent);
+ if (nps == _namedprojs.end()) {
+ ostringstream msg;
+ msg << "No projections registered for parent " << &parent;
+ throw Error(msg.str());
+ }
+ NamedProjs::const_iterator np = nps->second.find(name);
+ if (np == nps->second.end()) {
+ ostringstream msg;
+ msg << "No projection '" << name << "' found for parent " << &parent;
+ throw Error(msg.str());
+ }
+ // If it's registered with the projection handler, we must be able to safely
+ // dereference the Projection pointer to a reference...
+ return *(np->second);
+ }
+
+
+
+}
Copied: trunk/src/Core/Run.cc (from r2120, trunk/src/Run.cc)
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/src/Core/Run.cc Tue Dec 1 19:39:13 2009 (r2122, copy of r2120, trunk/src/Run.cc)
@@ -0,0 +1,155 @@
+// -*- C++ -*-
+#include "Rivet/Run.hh"
+#include "Rivet/AnalysisHandler.hh"
+#include "HepMC/IO_GenEvent.h"
+#include "Rivet/Projections/Beam.hh"
+#include <limits>
+
+namespace Rivet {
+
+
+ Run::Run(AnalysisHandler& ah)
+ : _ah(ah), _xs(-1.0)
+ { }
+
+
+ Run::~Run() { }
+
+
+ Run& Run::setCrossSection(const double xs) {
+ _xs = xs;
+ return *this;
+ }
+
+
+ Run& Run::setListAnalyses(const bool dolist) {
+ _listAnalyses = dolist;
+ return *this;
+ }
+
+
+ bool Run::prepareFile(const std::string& evtfile) {
+ if (evtfile == "-") {
+ _io.reset(new HepMC::IO_GenEvent(std::cin));
+ } else {
+ // Ignore the HepMC::IO_GenEvent(filename, ios) constructor, since it's only available from HepMC 2.4
+ _istr.reset(new std::fstream(evtfile.c_str(), std::ios::in));
+ _io.reset(new HepMC::IO_GenEvent(*_istr));
+ }
+ if (_io->rdstate() != 0) {
+ Log::getLog("Rivet.Run") << Log::ERROR << "Read error on file " << evtfile << endl;
+ return false;
+ }
+
+ return true;
+ }
+
+
+ bool Run::processEvent(bool firstEvent) {
+ // Fill event and check for a bad read state
+ shared_ptr<GenEvent> evt;
+ evt.reset(new GenEvent());
+ if (_io->rdstate() != 0 || !_io->fill_next_event(evt.get()) ) {
+ Log::getLog("Rivet.Run") << Log::DEBUG << "Read failed. End of file?" << endl;
+ return false;
+ }
+
+ // Get beam details from first event, and ensure they match for all following events
+ if (evt->particles_size() != 0) {
+ const BeamPair beams = beamIds(*evt);
+ const double sqrts = Rivet::sqrtS(*evt);
+ Log::getLog("Rivet.Run") << Log::DEBUG << "Beams: "
+ << beams << " @ " << sqrts/GeV << " GeV" << endl;
+ if (firstEvent) {
+ _beams = beams;
+ _sqrts = sqrts;
+ Log::getLog("Rivet.Run") << Log::INFO << "First event beams: "
+ << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
+ } else {
+ if (beams != _beams || !fuzzyEquals(sqrts, sqrtS())) {
+ Log::getLog("Rivet.Run") << Log::ERROR << "Event beams mismatch: "
+ << beams << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
+ << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl;
+ return false;
+ }
+ }
+ }
+
+ // Set up system based on properties of first event
+ if (firstEvent) {
+ // If empty
+ if (evt->particles_size() == 0) {
+ Log::getLog("Rivet.Run") << Log::ERROR << "Empty first event." << endl;
+ return false;
+ }
+
+ const size_t num_anas_requested = _ah.analysisNames().size();
+ _ah.removeIncompatibleAnalyses(beams());
+ if (num_anas_requested > 0 && _ah.analysisNames().size() == 0) {
+ Log::getLog("Rivet.Run") << Log::ERROR
+ << "All analyses were incompatible with the first event's beams\n"
+ << "Exiting, since this probably isn't intentional!" << endl;
+ return false;
+ }
+
+ if (_listAnalyses) {
+ foreach (const std::string& ana, _ah.analysisNames()) {
+ cout << ana << endl;
+ }
+ }
+
+ }
+
+
+ // Set cross-section if specified from command line
+ if (_xs > 0.0) {
+ if (firstEvent) {
+ Log::getLog("Rivet.Run") << Log::DEBUG
+ << "Setting user cross-section = " << _xs << " pb" << endl;
+ _ah.setCrossSection(_xs);
+ }
+ }
+ // Set cross-section if found in event
+ #ifdef HEPMC_HAS_CROSS_SECTION
+ else if (evt->cross_section()) {
+ const double xs = evt->cross_section()->cross_section(); //< in pb
+ Log::getLog("Rivet.Run") << Log::DEBUG
+ << "Setting cross-section = " << xs << " pb" << endl;
+ _ah.setCrossSection(xs);
+ }
+ #endif
+ // Complain about absence of cross-section if required!
+ else {
+ if (_ah.needCrossSection()) {
+ Log::getLog("Rivet.Run") << Log::ERROR
+ << "Total cross-section needed for at least one of the analyses. "
+ << "Please set it (on the command line with '-x' if using the 'rivet' program)" << endl;
+ return false;
+ }
+ }
+
+ /// @todo If NOT first event, check that beams aren't changed
+
+ // Analyze event
+ _ah.analyze(*evt);
+
+ return true;
+ }
+
+
+ bool Run::finalizeFile() {
+ return true;
+ }
+
+
+ const BeamPair& Run::beams() const {
+ return _beams;
+ }
+
+
+ double Run::sqrtS() const {
+ return _sqrts;
+ }
+
+
+}
Modified: trunk/src/Makefile.am
==============================================================================
--- trunk/src/Makefile.am Tue Dec 1 17:04:36 2009 (r2121)
+++ trunk/src/Makefile.am Tue Dec 1 19:39:13 2009 (r2122)
@@ -1,24 +1,21 @@
-SUBDIRS = . Projections Tools Test
+SUBDIRS = Core Projections Tools .
if ENABLE_ANALYSES
SUBDIRS += Analyses
endif
+SUBDIRS += Test
-lib_LTLIBRARIES = libRivet.la
+lib_LTLIBRARIES = libRivet.la
-libRivet_la_SOURCES = \
- Event.cc Jet.cc Particle.cc \
- ProjectionApplier.cc Projection.cc \
- Analysis.cc AnalysisLoader.cc AnalysisInfo.cc \
- AnalysisHandler.cc Run.cc ProjectionHandler.cc HistoHandler.cc
+libRivet_la_SOURCES =
-libRivet_la_CPPFLAGS = $(AM_CPPFLAGS) -I$(top_srcdir)/src/Tools/yaml-cpp $(CPPFLAGS)
-
-libRivet_la_LDFLAGS = -L$(prefix)/lib \
- -L$(GSLLIBPATH) -L$(GSLCBLASLIBPATH) -export-dynamic $(VERSIONINFOFLAGS)
-
-FASTJETLIBADD = $(FASTJETCONFIGLIBADD)
+libRivet_la_LDFLAGS = \
+ -L$(prefix)/lib \
+ -L$(GSLLIBPATH) -L$(GSLCBLASLIBPATH) \
+ -export-dynamic $(VERSIONINFOFLAGS)
libRivet_la_LIBADD = \
+ Core/libRivetCore.la \
Projections/libRivetProjections.la \
- Tools/libRivetTools.la -ldl -lm \
- $(GSL_LDFLAGS) $(FASTJETLIBADD)
+ Tools/libRivetTools.la \
+ -ldl -lm \
+ $(GSL_LDFLAGS) $(FASTJETCONFIGLIBADD)
More information about the Rivet-svn
mailing list
|