diff --git a/include/Rivet/Tools/RivetYODA.hh b/include/Rivet/Tools/RivetYODA.hh
--- a/include/Rivet/Tools/RivetYODA.hh
+++ b/include/Rivet/Tools/RivetYODA.hh
@@ -1,354 +1,350 @@
#ifndef RIVET_RIVETYODA_HH
#define RIVET_RIVETYODA_HH
/// @author Andy Buckley
/// @date 2009-01-30
/// @author David Grellscheid
/// @date 2011-07-18
/// @author David Grellscheid
/// @date 2016-09-27
#include "Rivet/Config/RivetCommon.hh"
#include "YODA/AnalysisObject.h"
#include "YODA/Counter.h"
#include "YODA/Histo1D.h"
#include "YODA/Histo2D.h"
#include "YODA/Profile1D.h"
#include "YODA/Profile2D.h"
#include "YODA/Scatter1D.h"
#include "YODA/Scatter2D.h"
#include "YODA/Scatter3D.h"
#include <map>
namespace YODA {
typedef std::shared_ptr<YODA::AnalysisObject> AnalysisObjectPtr;
typedef std::shared_ptr<YODA::Histo1D> Histo1DPtr;
typedef std::shared_ptr<YODA::Histo2D> Histo2DPtr;
typedef std::shared_ptr<YODA::Profile1D> Profile1DPtr;
typedef std::shared_ptr<YODA::Profile2D> Profile2DPtr;
typedef std::shared_ptr<YODA::Scatter1D> Scatter1DPtr;
typedef std::shared_ptr<YODA::Scatter2D> Scatter2DPtr;
typedef std::shared_ptr<YODA::Scatter3D> Scatter3DPtr;
typedef std::shared_ptr<YODA::Counter> CounterPtr;
}
namespace Rivet {
class AnalysisObjectPtr {
public:
virtual ~AnalysisObjectPtr() {}
virtual YODA::AnalysisObject* operator->() = 0;
virtual YODA::AnalysisObject* operator->() const = 0;
virtual const YODA::AnalysisObject & operator*() const = 0;
bool operator ==(const AnalysisObjectPtr& p) { return (this == &p); }
protected:
/// @todo do we need this?
// virtual void reset() = 0;
};
/// @todo
/// implement scatter1dptr and scatter2dptr here
class Scatter1DPtr : public AnalysisObjectPtr {
public:
Scatter1DPtr() :
_scatter(YODA::Scatter1DPtr()) { }
Scatter1DPtr(const YODA::Scatter1DPtr& p) :
_scatter(p) { }
bool operator!() const { return !_scatter; }
operator bool() const { return bool(_scatter); }
YODA::Scatter1D* operator->() { return _scatter.get(); }
YODA::Scatter1D* operator->() const { return _scatter.get(); }
YODA::Scatter1D & operator*() { return *_scatter; }
const YODA::Scatter1D & operator*() const { return *_scatter; }
protected:
YODA::Scatter1DPtr _scatter;
};
class Scatter2DPtr : public AnalysisObjectPtr {
public:
Scatter2DPtr(const YODA::Scatter2DPtr& p) :
_scatter(p) { }
Scatter2DPtr() :
_scatter(YODA::Scatter2DPtr()) { }
bool operator!() { return !_scatter; }
operator bool() { return bool(_scatter); }
YODA::Scatter2D* operator->() { return _scatter.get(); }
YODA::Scatter2D* operator->() const { return _scatter.get(); }
YODA::Scatter2D & operator*() { return *_scatter; }
const YODA::Scatter2D & operator*() const { return *_scatter; }
protected:
YODA::Scatter2DPtr _scatter;
};
class Scatter3DPtr : public AnalysisObjectPtr {
public:
Scatter3DPtr(const YODA::Scatter3DPtr& p) :
_scatter(p) { }
Scatter3DPtr() :
_scatter(YODA::Scatter3DPtr()) { }
bool operator!() { return !_scatter; }
operator bool() { return bool(_scatter); }
YODA::Scatter3D* operator->() { return _scatter.get(); }
YODA::Scatter3D* operator->() const { return _scatter.get(); }
YODA::Scatter3D & operator*() { return *_scatter; }
const YODA::Scatter3D & operator*() const { return *_scatter; }
protected:
YODA::Scatter3DPtr _scatter;
};
class MultiweightAOPtr : public AnalysisObjectPtr {
public:
virtual void newSubEvent() = 0;
/// @todo
/// rename to setActive(Idx)?
virtual void setActiveWeightIdx(unsigned int iWeight) = 0;
virtual void pushToPersistent(const vector<vector<double> >& weight) = 0;
virtual YODA::AnalysisObjectPtr activeYODAPtr() const = 0;
};
template <class T>
using Fill = pair<typename T::FillType, double>;
template <class T>
using Fills = multiset<Fill<T>>;
+// TODO TODO TODO
+// need to override the old fill method too!
+// otherwise we can't intercept existing fill calls in analysis code
+//
+// need to
+// TODO TODO TODO
+
template <class T>
class TupleWrapper : public T {
public:
typedef shared_ptr<TupleWrapper> Ptr;
TupleWrapper(const T & h) : T(h) {}
- using T::fill;
-
// todo: do we need to deal with users using fractions directly?
void fill(typename T::FillType x, double weight=1.0, double fraction=1.0) {
fills_.insert( {x, weight} );
}
void reset() {
fills_.clear();
}
const Fills<T> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<T> fills_;
};
template <class T>
class Wrapper : public MultiweightAOPtr {
public:
/* @todo
* some things are not really well-defined here
* for instance: fill() in the finalize() method and integral() in
* the analyze() method.
*/
Wrapper(size_t len_of_weightvec, const T & p) {
for (size_t i = 0; i < len_of_weightvec; i++)
_persistent.push_back(make_shared<T>(p));
}
typename T::Ptr active() const { return _active; }
/* @todo this probably need to loop over all? */
bool operator!() const { return !active(); }
operator bool() const { return bool(active()); }
T * operator->() {
return active().get();
}
T * operator->() const {
return active().get();
}
T & operator*() {
return *active();
}
const T & operator*() const {
return *active();
}
/* @todo
* these need to be re-thought out.
void reset() { active()->reset(); }
*/
/* @todo
* these probably need to loop over all?
* do we even want to provide equality?
*/
/* @todo
* how about no.
friend bool operator==(Wrapper a, Wrapper b){
if (a._persistent.size() != b._persistent.size())
return false;
for (size_t i = 0; i < a._persistent.size(); i++) {
if (a._persistent.at(i) != b._persistent.at(i)) {
return false;
}
}
return true;
}
friend bool operator!=(Wrapper a, Wrapper b){
return !(a == b);
}
friend bool operator<(Wrapper a, Wrapper b){
if (a._persistent.size() >= b._persistent.size())
return false;
for (size_t i = 0; i < a._persistent.size(); i++) {
if (*(a._persistent.at(i)) >= *(b._persistent.at(i))) {
return false;
}
}
return true;
}
*/
private:
void setActiveWeightIdx(unsigned int iWeight) {
_active = _persistent.at(iWeight);
}
/* this is for dev only---we shouldn't need this in real runs. */
void unsetActiveWeight() {
_active.reset();
}
void newSubEvent() {
typename TupleWrapper<T>::Ptr tmp
= make_shared<TupleWrapper<T>>(_persistent[0]->clone());
tmp->reset();
_evgroup.push_back( tmp );
_active = _evgroup.back();
}
virtual YODA::AnalysisObjectPtr activeYODAPtr() const {
return _active;
}
const vector<typename T::Ptr> & persistent() const {
return _persistent;
}
/* to be implemented for each type */
void pushToPersistent(const vector<vector<double> >& weight);
/* M of these, one for each weight */
vector<typename T::Ptr> _persistent;
/* N of these, one for each event in evgroup */
vector<typename TupleWrapper<T>::Ptr> _evgroup;
typename T::Ptr _active;
friend class AnalysisHandler;
};
-
-using Histo1DTuple = TupleWrapper<YODA::Histo1D>;
-
-
-
// every object listed here needs a virtual fill method in YODA,
// otherwise the Tuple fakery won't work.
using Histo1DPtr = Wrapper<YODA::Histo1D>;
-
- using Histo2DPtr = Wrapper<YODA::Histo2D>;
-
- //RIVETAOPTR_COMMON(Histo1D)
- // RIVETAOPTR_COMMON(Histo2D)
- // RIVETAOPTR_COMMON(Profile1D)
- // RIVETAOPTR_COMMON(Profile2D)
- // RIVETAOPTR_COMMON(Counter)
+ //using Histo2DPtr = Wrapper<YODA::Histo2D>;
+ using Profile1DPtr = Wrapper<YODA::Profile1D>;
+ using Profile2DPtr = Wrapper<YODA::Profile2D>;
+ using CounterPtr = Wrapper<YODA::Counter>;
using YODA::Counter;
using YODA::Histo1D;
using YODA::HistoBin1D;
using YODA::Histo2D;
using YODA::HistoBin2D;
using YODA::Profile1D;
using YODA::ProfileBin1D;
using YODA::Profile2D;
using YODA::ProfileBin2D;
using YODA::Scatter1D;
using YODA::Point1D;
using YODA::Scatter2D;
using YODA::Point2D;
using YODA::Scatter3D;
using YODA::Point3D;
/// Function to get a map of all the refdata in a paper with the
/// given @a papername.
map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername);
/// @todo Also provide a Scatter3D getRefData() version?
/// Get the file system path to the reference file for this paper.
string getDatafilePath(const string& papername);
/// Return the integral over the histogram bins
/// @deprecated Prefer to directly use the histo's integral() method.
DEPRECATED("Prefer to directly use the histo's integral() method.")
inline double integral(Histo1DPtr histo) {
return histo->integral();
}
}
#endif
diff --git a/src/Tools/RivetYODA.cc b/src/Tools/RivetYODA.cc
--- a/src/Tools/RivetYODA.cc
+++ b/src/Tools/RivetYODA.cc
@@ -1,287 +1,300 @@
#include "Rivet/Config/RivetCommon.hh"
#include "Rivet/Tools/RivetYODA.hh"
#include "Rivet/Tools/RivetPaths.hh"
#include "YODA/ReaderYODA.h"
#include "YODA/ReaderAIDA.h"
using namespace std;
namespace Rivet {
string getDatafilePath(const string& papername) {
/// Try to find YODA otherwise fall back to try AIDA
const string path1 = findAnalysisRefFile(papername + ".yoda");
if (!path1.empty()) return path1;
const string path2 = findAnalysisRefFile(papername + ".aida");
if (!path2.empty()) return path2;
throw Rivet::Error("Couldn't find ref data file '" + papername + ".yoda/aida" +
" in $RIVET_REF_PATH, '" + getRivetDataPath() + "', or '.'");
}
map<string, YODA::AnalysisObjectPtr> getRefData(const string& papername) {
const string datafile = getDatafilePath(papername);
// Make an appropriate data file reader and read the data objects
/// @todo Remove AIDA support some day...
YODA::Reader& reader = (datafile.find(".yoda") != string::npos) ? \
YODA::ReaderYODA::create() : YODA::ReaderAIDA::create();
vector<YODA::AnalysisObject *> aovec;
reader.read(datafile, aovec);
// Return value, to be populated
map<string, YODA::AnalysisObjectPtr> rtn;
foreach ( YODA::AnalysisObject* ao, aovec ) {
YODA::AnalysisObjectPtr refdata(ao);
if (!refdata) continue;
const string plotpath = refdata->path();
// Split path at "/" and only return the last field, i.e. the histogram ID
const size_t slashpos = plotpath.rfind("/");
const string plotname = (slashpos+1 < plotpath.size()) ? plotpath.substr(slashpos+1) : "";
rtn[plotname] = refdata;
}
return rtn;
}
}
namespace {
using Rivet::Fill;
using Rivet::Fills;
template <class T>
void fillAllPersistent(const vector<typename T::Ptr> & persistent,
typename T::FillType x, double w,
const vector<double> & weight) {
for ( size_t m = 0; m < persistent.size(); ++m ) {
persistent[m]->fill( x, w * weight[m] );
}
}
template <class T>
double get_window_size(const typename T::Ptr & histo,
typename T::BinType x) {
// the bin index we fall in
const auto binidx = histo->binIndexAt(x);
// gaps, overflow, underflow don't contribute
if ( binidx == -1 )
return 0;
const auto & b = histo->bin(binidx);
// if we don't have a valid neighbouring bin,
// we use infinite width
typename T::Bin b1(-1.0/0.0, 1.0/0.0);
// points in the top half compare to the upper neighbour
if ( x > b.xMid() ) {
int nextidx = binidx + 1;
if ( nextidx < histo->bins().size() )
b1 = histo->bin(nextidx);
}
else { // compare to the lower neighbour
int nextidx = binidx - 1;
if ( nextidx >= 0 )
b1 = histo->bin(nextidx);
}
// the factor 2 is arbitrary, could poss. be smaller
return min( b.width(), b1.width() ) / 2.0;
}
// template <class T>
// double getX(T a) {
// return a;
// }
// template <>
// double getX<tuple<double,double> >(tuple<double,double> a) {
// return get<0>(a);
// }
template <class T>
void commit(typename T::Ptr & persistent,
const vector< vector<Fill<T> > > & tuple,
const vector<double> weights /* generator weight over subevents */) {
for ( const auto & x : tuple ) {
double maxwindow = 0.0;
for ( const auto & xi : x ) {
// check for NOFILL here
double window = get_window_size<T>(persistent, xi.first);
if ( window > maxwindow )
maxwindow = window;
}
const double wsize = maxwindow;
// all windows have same size
set<double> edgeset;
// bin edges need to be in here!
for ( const auto & xi : x ) {
edgeset.insert((xi.first) - wsize);
edgeset.insert((xi.first) + wsize);
}
vector< std::tuple<double,double,double> > hfill;
double sumf = 0.0;
auto edgit = edgeset.begin();
double ehi = *edgit;
while ( ++edgit != edgeset.end() ) {
double elo = ehi;
ehi = *edgit;
double sumw = 0.0;
bool gap = true; // Check for gaps between the sub-windows.
for ( int i = 0; i < x.size(); ++i ) {
// check equals comparisons here!
if ( (x[i].first) + wsize >= ehi && (x[i].first) - wsize <= elo ) {
sumw += x[i].second * weights[i];
gap = false;
}
}
if ( gap ) continue;
hfill.push_back(make_tuple((ehi + elo)/2.0, sumw, ehi - elo));
sumf += ehi - elo;
}
for ( auto f : hfill )
persistent->fill( get<0>(f), get<1>(f), get<2>(f)/sumf );
// Note the scaling to one single fill
}
}
template <class T>
double distance(T a, T b) {
return abs(a - b);
}
template <>
double distance<tuple<double,double> >(tuple<double,double> a, tuple<double,double> b) {
return Rivet::sqr(get<0>(a) - get<0>(b)) + Rivet::sqr(get<1>(a) - get<1>(b));
}
}
/// fills is a vector of sub-event with an ordered set of x-values of
/// the fills in each sub-event. NOFILL should be an "impossible"
/// value for this histogram. Returns a vector of sub-events with
/// an ordered vector of fills (including NOFILLs) for each sub-event.
template <class T>
vector< vector<Fill<T> > >
match_fills(const vector< Fills<T> > & fills, const Fill<T> & NOFILL)
{
vector< vector<Fill<T> > > matched;
// First just copy subevents into vectors and find the longest vector.
unsigned int maxfill = 0; // length of biggest vector
int imax = 0; // index position of biggest vector
for ( const auto & subev : fills ) {
if ( subev.size() > maxfill ) {
maxfill = subev.size();
imax = matched.size();
}
matched.push_back(vector<Fill<T> >(subev.begin(), subev.end()));
}
// Now, go through all subevents with missing fills.
const vector<Fill<T>> & full = matched[imax]; // the longest one
for ( auto & subev : matched ) {
if ( subev.size() == maxfill ) continue;
// Add NOFILLs to the end;
while ( subev.size() < maxfill ) subev.push_back(NOFILL);
// Iterate from the back and shift all fill values backwards by
// swapping with NOFILLs so that they better match the full
// subevent.
for ( int i = maxfill - 1; i >= 0; --i ) {
if ( subev[i] == NOFILL ) continue;
int j = i;
while ( j + 1 < maxfill && subev[j + 1] == NOFILL &&
distance(subev[j].first, full[j].first) > distance(subev[j].first, full[j + 1].first) )
{
swap(subev[j], subev[j + 1]);
++j;
}
}
}
// transpose
vector<vector<Fill<T>>> result(maxfill,vector<Fill<T>>(matched.size()));
for (size_t i = 0; i < matched.size(); ++i)
for (size_t j = 0; j < maxfill; ++j)
result.at(j).at(i) = matched.at(i).at(j);
return result;
}
namespace Rivet {
template <class T>
void Wrapper<T>::pushToPersistent(const vector<vector<double> >& weight) {
// have we had subevents at all?
const bool have_subevents = _evgroup.size() > 1;
if ( ! have_subevents ) {
assert( _evgroup.size() == 1 && weight.size() == 1 );
// simple replay of all tuple entries
// each recorded fill is inserted into all persistent weightname histos
for ( const auto & f : _evgroup[0]->fills() ) {
fillAllPersistent<T>( _persistent, f.first, f.second, weight[0] );
}
} else {
assert( _evgroup.size() == weight.size() );
// All the fills across subevents
// each item in allFills is a subevent
vector<Fills<T>> allFills;
for ( const auto & ev : _evgroup )
allFills.push_back( ev->fills() );
vector< vector<Fill<T> > >
linedUpXs = match_fills<T>(allFills, {typename T::FillType(), 0.0});
for ( size_t m = 0; m < _persistent.size(); ++m ) {
commit<T>( _persistent[m], linedUpXs, weight[m] );
}
}
-
_evgroup.clear();
_active.reset();
}
- static Histo1DPtr foobar1 = Histo1DPtr(13, YODA::Histo1D() );
+ template <>
+ void Wrapper<YODA::Counter>::pushToPersistent(
+ const vector<vector<double> >& weight) {
+ for ( size_t m = 0; m < _persistent.size(); ++m ) {
+ for ( const auto & f : _evgroup[m]->fills() ) {
+ fillAllPersistent<YODA::Counter>( _persistent, f.first, f.second, weight[m] );
+ }
+ }
+ }
+
+ static auto foobar1 = Histo1DPtr(13, YODA::Histo1D() );
+
+ static auto foobar2 = CounterPtr(13, YODA::Counter() );
+
+ //static auto foobar3 = Profile1DPtr(13, YODA::Profile1D() );
// static Histo2DPtr foobar2 = Histo2DPtr(13, YODA::Histo2D() );
// void Histo2DPtr::pushToPersistent(const vector<vector<double> >& weight) {
// /// @todo
// return;
// }
// void Profile1DPtr::pushToPersistent(const vector<vector<double> >& weight) {
// /// @todo
// return;
// }
// void Profile2DPtr::pushToPersistent(const vector<vector<double> >& weight) {
// /// @todo
// return;
// }
// void CounterPtr::pushToPersistent(const vector<vector<double> >& weight) {
// /// @todo
// return;
// }
}