Page MenuHomeHEPForge
Files F7878516

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,406 +1,386 @@
#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
/// these need to be multi-weighted eventually.
class Scatter1DPtr : public AnalysisObjectPtr {
public:
Scatter1DPtr() :
_scatter(YODA::Scatter1DPtr()) { }
Scatter1DPtr(const YODA::Scatter1D& p) :
_scatter(make_shared<YODA::Scatter1D>(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::Scatter2D& p) :
_scatter(make_shared<YODA::Scatter2D>(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::Scatter3D& p) :
_scatter(make_shared<YODA::Scatter3D>(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<valarray<double> >& weight) = 0;
virtual YODA::AnalysisObjectPtr activeYODAPtr() const = 0;
};
using Weight = double;
template <class T>
using Fill = pair<typename T::FillType, Weight>;
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;
template<>
+class TupleWrapper<YODA::Counter> : public YODA::Counter {
+public:
+ typedef shared_ptr<TupleWrapper<YODA::Counter>> Ptr;
+ TupleWrapper(const YODA::Counter & h) : YODA::Counter(h) {}
+ // todo: do we need to deal with users using fractions directly?
+ void fill( double weight=1.0, double fraction=1.0 ) {
+ fills_.insert( {YODA::Counter::FillType(),weight} );
+ }
+ void reset() { fills_.clear(); }
+ const Fills<YODA::Counter> & fills() const { return fills_; }
+private:
+ // x / weight pairs
+ Fills<YODA::Counter> fills_;
+};
+
+template<>
class TupleWrapper<YODA::Histo1D> : public YODA::Histo1D {
public:
typedef shared_ptr<TupleWrapper<YODA::Histo1D>> Ptr;
TupleWrapper(const YODA::Histo1D & h) : YODA::Histo1D(h) {}
// todo: do we need to deal with users using fractions directly?
void fill( double x, double weight=1.0, double fraction=1.0 ) {
+ if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
fills_.insert( { x , weight } );
}
void reset() { fills_.clear(); }
const Fills<YODA::Histo1D> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<YODA::Histo1D> fills_;
};
template<>
class TupleWrapper<YODA::Profile1D> : public YODA::Profile1D {
public:
typedef shared_ptr<TupleWrapper<YODA::Profile1D>> Ptr;
TupleWrapper(const YODA::Profile1D & h) : YODA::Profile1D(h) {}
// todo: do we need to deal with users using fractions directly?
void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
+ if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
+ if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
fills_.insert( { {x,y}, weight } );
}
void reset() { fills_.clear(); }
const Fills<YODA::Profile1D> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<YODA::Profile1D> fills_;
};
-template<>
-class TupleWrapper<YODA::Counter> : public YODA::Counter {
-public:
- typedef shared_ptr<TupleWrapper<YODA::Counter>> Ptr;
- TupleWrapper(const YODA::Counter & h) : YODA::Counter(h) {}
- // todo: do we need to deal with users using fractions directly?
- void fill( double weight=1.0, double fraction=1.0 ) {
- fills_.insert( {YODA::Counter::FillType(),weight} );
- }
- void reset() { fills_.clear(); }
- const Fills<YODA::Counter> & fills() const { return fills_; }
-private:
- // x / weight pairs
- Fills<YODA::Counter> fills_;
-};
template<>
class TupleWrapper<YODA::Histo2D> : public YODA::Histo2D {
public:
typedef shared_ptr<TupleWrapper<YODA::Histo2D>> Ptr;
TupleWrapper(const YODA::Histo2D & h) : YODA::Histo2D(h) {}
// todo: do we need to deal with users using fractions directly?
void fill( double x, double y, double weight=1.0, double fraction=1.0 ) {
+ if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
+ if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
fills_.insert( {{x,y}, weight} );
}
void reset() { fills_.clear(); }
const Fills<YODA::Histo2D> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<YODA::Histo2D> fills_;
};
template<>
class TupleWrapper<YODA::Profile2D> : public YODA::Profile2D {
public:
typedef shared_ptr<TupleWrapper<YODA::Profile2D>> Ptr;
TupleWrapper(const YODA::Profile2D & h) : YODA::Profile2D(h) {}
// todo: do we need to deal with users using fractions directly?
void fill( double x, double y, double z, double weight=1.0, double fraction=1.0 ) {
+ if ( std::isnan(x) ) throw YODA::RangeError("X is NaN");
+ if ( std::isnan(y) ) throw YODA::RangeError("Y is NaN");
+ if ( std::isnan(z) ) throw YODA::RangeError("Z is NaN");
fills_.insert( {{x,y,z}, weight} );
}
void reset() { fills_.clear(); }
const Fills<YODA::Profile2D> & fills() const { return fills_; }
private:
// x / weight pairs
Fills<YODA::Profile2D> 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 m = 0; m < len_of_weightvec; ++m)
- _persistent.push_back(make_shared<T>(p));
- }
+ Wrapper() = default;
- Wrapper() : _persistent(), _evgroup(), _active() {}
+ Wrapper(size_t len_of_weightvec, const T & p);
- typename T::Ptr active() const { assert(_active); return _active; }
+ typename T::Ptr active() const;
/* @todo this probably need to loop over all? */
bool operator!() const { return !active(); }
- operator bool() const { return bool(active()); }
- T * operator->() {
- return active().get();
- }
+ operator bool() const { return static_cast<bool>(active()); }
- T * operator->() const {
- return active().get();
- }
+ T * operator->() { return active().get(); }
+
+ T * operator->() const { return active().get(); }
- T & operator*() {
- return *active();
- }
+ T & operator*() { return *active(); }
- const T & operator*() const {
- 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 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();
- assert(_active);
- }
+ void newSubEvent();
- virtual YODA::AnalysisObjectPtr activeYODAPtr() const {
- return _active;
- }
+ virtual YODA::AnalysisObjectPtr activeYODAPtr() const { return _active; }
- const vector<typename T::Ptr> & persistent() const {
- return _persistent;
- }
+ const vector<typename T::Ptr> & persistent() const { return _persistent; }
/* to be implemented for each type */
void pushToPersistent(const vector<valarray<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;
};
// 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>;
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);
}
#endif
diff --git a/src/Tools/RivetYODA.cc b/src/Tools/RivetYODA.cc
--- a/src/Tools/RivetYODA.cc
+++ b/src/Tools/RivetYODA.cc
@@ -1,301 +1,335 @@
#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 {
+template <class T>
+Wrapper<T>::Wrapper(size_t len_of_weightvec, const T & p) {
+ for (size_t m = 0; m < len_of_weightvec; ++m)
+ _persistent.push_back(make_shared<T>(p));
+}
+
+
+template <class T>
+typename T::Ptr Wrapper<T>::active() const {
+ assert(_active);
+ return _active;
+}
+
+
+template <class T>
+void Wrapper<T>::newSubEvent() {
+ typename TupleWrapper<T>::Ptr tmp
+ = make_shared<TupleWrapper<T>>(_persistent[0]->clone());
+ tmp->reset();
+ _evgroup.push_back( tmp );
+ _active = _evgroup.back();
+ assert(_active);
+}
+
+
+
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;
using Rivet::TupleWrapper;
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>
typename T::BinType
fillT2binT(typename T::FillType a) {
return a;
}
template <>
YODA::Profile1D::BinType
fillT2binT<YODA::Profile1D>(YODA::Profile1D::FillType a) {
return get<0>(a);
}
template <>
YODA::Profile2D::BinType
fillT2binT<YODA::Profile2D>(YODA::Profile2D::FillType a) {
return { get<0>(a), get<1>(a) };
}
template <class T>
void commit(vector<typename T::Ptr> & persistent,
const vector< vector<Fill<T>> > & tuple,
const vector<valarray<double>> & weights ) {
// TODO check if all the xs are in the same bin anyway!
// Then no windowing needed
assert(persistent.size() == weights[0].size());
for ( const auto & x : tuple ) {
double maxwindow = 0.0;
for ( const auto & xi : x ) {
// TODO check for NOFILL here
// persistent[0] has the same binning as all the other persistent objects
double window = get_window_size<T>(persistent[0], fillT2binT<T>(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(fillT2binT<T>(xi.first) - wsize);
edgeset.insert(fillT2binT<T>(xi.first) + wsize);
}
vector< std::tuple<double,valarray<double>,double> > hfill;
double sumf = 0.0;
auto edgit = edgeset.begin();
double ehi = *edgit;
while ( ++edgit != edgeset.end() ) {
double elo = ehi;
ehi = *edgit;
valarray<double> sumw(0.0, persistent.size()); // need m copies of this
bool gap = true; // Check for gaps between the sub-windows.
for ( int i = 0; i < x.size(); ++i ) {
// check equals comparisons here!
if ( fillT2binT<T>(x[i].first) + wsize >= ehi
&&
fillT2binT<T>(x[i].first) - wsize <= elo ) {
sumw += x[i].second * weights[i];
gap = false;
}
}
if ( gap ) continue;
hfill.push_back( { (ehi + elo)/2.0, sumw, (ehi - elo) } );
sumf += ehi - elo;
}
for ( auto f : hfill )
for ( size_t m = 0; m < persistent.size(); ++m )
persistent[m]->fill( get<0>(f), get<1>(f)[m], get<2>(f)/sumf );
// Note the scaling to one single fill
}
}
template<>
void commit<YODA::Histo2D>(vector<YODA::Histo2D::Ptr> & persistent,
const vector< vector<Fill<YODA::Histo2D>> > & tuple,
const vector<valarray<double>> & weights)
{}
template<>
void commit<YODA::Profile2D>(vector<YODA::Profile2D::Ptr> & persistent,
const vector< vector<Fill<YODA::Profile2D>> > & tuple,
const vector<valarray<double>> & weights)
{}
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<typename TupleWrapper<T>::Ptr> & evgroup, 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 & it : evgroup ) {
const auto & subev = it->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(fillT2binT<T>(subev[j].first),
fillT2binT<T>(full[j].first))
>
distance(fillT2binT<T>(subev[j].first),
fillT2binT<T>(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<valarray<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 ( size_t m = 0; m < _persistent.size(); ++m )
for ( const auto & f : _evgroup[0]->fills() )
_persistent[m]->fill( f.first, f.second * weight[0][m] );
+
+
} else {
+
+
assert( _evgroup.size() == weight.size() );
// outer index is subevent, inner index is jets in the event
vector<vector<Fill<T>>> linedUpXs
= match_fills<T>(_evgroup, {typename T::FillType(), 0.0});
commit<T>( _persistent, linedUpXs, weight );
+
+
}
_evgroup.clear();
_active.reset();
}
template <>
void Wrapper<YODA::Counter>::pushToPersistent(const vector<valarray<double> >& weight) {
for ( size_t m = 0; m < _persistent.size(); ++m ) {
for ( size_t n = 0; n < _evgroup.size(); ++n ) {
for ( const auto & f : _evgroup[n]->fills() ) {
_persistent[m]->fill( f.second * weight[n][m] );
}
}
}
}
// explicitly instantiate all wrappers
template class Wrapper<YODA::Histo1D>;
template class Wrapper<YODA::Histo2D>;
template class Wrapper<YODA::Profile1D>;
template class Wrapper<YODA::Profile2D>;
template class Wrapper<YODA::Counter>;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:11 PM (1 d, 19 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805532
Default Alt Text
(24 KB)

Event Timeline