Page MenuHomeHEPForge
Files F7879311

No OneTemporary
Actions

View Options

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,353 +1,354 @@
#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>>;
-class Histo1DTuple : public YODA::Histo1D {
+template <class T>
+class TupleWrapper : public T {
public:
- Histo1DTuple(const YODA::Histo1D & h) : YODA::Histo1D(h) {}
+ 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(double x, double weight=1.0, double fraction=1.0) {
+ void fill(typename T::FillType x, double weight=1.0, double fraction=1.0) {
fills_.insert( {x, weight} );
}
void reset() {
fills_.clear();
}
- const Fills<YODA::Histo1D> & fills() const { return fills_; }
+ const Fills<T> & fills() const { return fills_; }
private:
// x / weight pairs
- Fills<YODA::Histo1D> fills_;
+ 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();
+ }
-#define RIVETAOPTR_COMMON(YODATYPE) \
- typedef shared_ptr<YODATYPE##Tuple> YODATYPE##TuplePtr; \
- \
- class YODATYPE##Ptr : 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. \
- */ \
- \
- YODATYPE##Ptr(size_t len_of_weightvec, const YODA::YODATYPE& p) { \
- for (size_t i = 0; i < len_of_weightvec; i++) \
- _persistent.push_back(make_shared<YODA::YODATYPE>(p)); \
- \
- return; \
- } \
- \
- YODA::YODATYPE##Ptr active() const { return _active; } \
- \
- /* @todo this probably need to loop over all? */ \
- bool operator!() const {return !active();} \
- operator bool() const {return bool(active());} \
- \
- YODA::YODATYPE* operator->() { \
- return active().get(); \
- } \
- \
- YODA::YODATYPE* operator->() const { \
- return active().get(); \
- } \
- \
- YODA::YODATYPE & operator*() { \
- return *active(); \
- } \
- \
- const YODA::YODATYPE & 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==(YODATYPE##Ptr a, YODATYPE##Ptr 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!=(YODATYPE##Ptr a, YODATYPE##Ptr b){ \
- return !(a == b); \
- } \
- \
- \
- friend bool operator<(YODATYPE##Ptr a, YODATYPE##Ptr 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); \
- return; \
- } \
- \
- /* this is for dev only---we shouldn't need this in real runs. */ \
- void unsetActiveWeight() { \
- _active.reset(); \
- return; \
- } \
- \
- void newSubEvent() { \
- YODATYPE##TuplePtr tmp \
- = make_shared<YODATYPE##Tuple>(_persistent[0]->clone()); \
- tmp->reset(); \
- _evgroup.push_back( tmp ); \
- _active = _evgroup.back(); \
- return; \
- } \
- \
- virtual YODA::AnalysisObjectPtr activeYODAPtr() const { \
- return _active; \
- } \
- \
- const vector<YODA::YODATYPE##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<YODA::YODATYPE##Ptr> _persistent; \
- \
- /* N of these, one for each event in evgroup */ \
- vector<YODATYPE##TuplePtr> _evgroup; \
- \
- YODA::YODATYPE##Ptr _active; \
- \
- friend class AnalysisHandler; \
+ 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.
- RIVETAOPTR_COMMON(Histo1D)
+ 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 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,266 +1,287 @@
#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,
- double x) {
+ 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);
+ 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 ) {
+ 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>
- T distance(T a, T b) {
+ 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 {
- using T = YODA::Histo1D;
-
- void Histo1DPtr::pushToPersistent(const vector<vector<double> >& weight) {
+ 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, {T::FillType(),
- 0.0});
+ 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() );
+
+ // 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;
// }
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:54 PM (1 d, 4 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805934
Default Alt Text
(30 KB)

Event Timeline