Page MenuHomeHEPForge

No OneTemporary

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/Exsample2/exsample/binary_tree.h b/Exsample2/exsample/binary_tree.h
--- a/Exsample2/exsample/binary_tree.h
+++ b/Exsample2/exsample/binary_tree.h
@@ -1,883 +1,883 @@
// -*- C++ -*-
//
// binary_tree.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
//
// ExSample is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
#ifndef EXSAMPLE_binary_tree_h_included
#define EXSAMPLE_binary_tree_h_included
#include "utility.h"
namespace exsample {
/// \brief binary_tree represents a binary tree with the ability to
/// `cascade' visitor objects down the tree
template<class Value>
class binary_tree {
public:
///@name type definitions
//@{
/// define the object type
typedef Value value_type;
//@}
public:
///@name constructors
//@{
/// default constructor
binary_tree()
: neighbours_(0,0),
parent_(0), value_(),
children_()
{ }
/// construct giving key/cell and parent
binary_tree(const value_type& thevalue,
binary_tree * theparent = 0)
: neighbours_(0,0), parent_(theparent),
value_(new value_type(thevalue)),
children_()
{ }
/// binary_tree has a strict ownership; on copying
/// binary trees ownership is transferred
binary_tree(const binary_tree& x)
: neighbours_(x.neighbours_),
parent_(x.parent_), value_(),
children_() {
assert(x.root());
binary_tree& nc_x = const_cast<binary_tree&>(x);
value_.swap(nc_x.value_);
children_.first.swap(nc_x.children_.first);
children_.second.swap(nc_x.children_.second);
nc_x.parent_ = 0;
nc_x.neighbours_.first = 0;
nc_x.neighbours_.second = 0;
}
/// binary_tree has a strict ownership; on copying
/// binary trees ownership is transferred
binary_tree& operator=(const binary_tree& x) {
if (this == &x)
return *this;
assert(x.root());
binary_tree& nc_x = const_cast<binary_tree&>(x);
value_.swap(nc_x.value_);
children_.first.swap(nc_x.children_.first);
children_.second.swap(nc_x.children_.second);
neighbours_ = x.neighbours_;
parent_ = x.parent_;
nc_x.parent_ = 0;
nc_x.neighbours_.first = 0;
nc_x.neighbours_.second = 0;
return *this;
}
//@}
///@name standard-conforming leaf iterators
//@{
public:
class const_iterator;
/// iterator
class iterator {
public:
///@name type definitions for iterator traits
//@{
/// define the iterator category
typedef std::bidirectional_iterator_tag iterator_category;
/// define the difference_type
typedef int difference_type;
/// define the value type
typedef Value value_type;
/// define the reference type
typedef value_type& reference;
/// define the pointer type
typedef value_type * pointer;
//@}
public:
///@name constructors
//@{
/// default constructor
iterator() : pointee(0), post_end(0), pre_begin(0) { }
/// constructor taking pointee
iterator(binary_tree * p, std::size_t end = 0)
: pointee(p), post_end(end), pre_begin(0) { }
//@
public:
///@name comparisons
//@{
/// comparison
bool operator==(const iterator& x) const {
return ((pointee == x.pointee) &&
(post_end == x.post_end) &&
(pre_begin == x.pre_begin));
}
/// comparison
bool operator!=(const iterator& x) const { return !(*this == x); }
//@}
public:
///@name derefrence and indirection
//@{
/// dereference
reference operator*() { return pointee->value(); }
/// indirection
pointer operator->() { return &**this; }
/// return reference to the node
binary_tree& node() { return *pointee; }
//@}
/// return raw pointer to the element pointed to
binary_tree * get() const { return pointee; }
///@name biderectional iterator increment/decrements
//@{
/// pre-increment
iterator& operator++() {
if (post_end) { ++post_end; return *this; }
if (pre_begin) { --pre_begin; return *this; }
if(!(pointee->right_neighbour())) { post_end = 1; return *this; }
pointee = pointee->right_neighbour();
return *this;
}
/// pre-decrement
iterator& operator--() {
if (post_end) { --post_end; return *this; }
if (pre_begin) { ++pre_begin; return *this; }
if(!(pointee->left_neighbour())) { pre_begin = 1; return *this; }
pointee = pointee->left_neighbour();
return *this;
}
/// post-increment
iterator operator++(int) {
iterator tmp = *this;
++(*this);
return tmp;
}
/// post-decrement
iterator operator--(int) {
iterator tmp = *this;
--(*this);
return tmp;
}
//@}
private:
/// friend for conversion
friend class const_iterator;
/// the node pointed to
binary_tree * pointee;
/// the distance from --end() (if above --end())
std::size_t post_end;
/// the distance from begin() (if below begin())
std::size_t pre_begin;
};
/// return begin iterator
iterator begin() { return iterator(left_most()); }
/// return end iterator
iterator end() { return iterator(right_most(),1); }
/// return global begin iterator
iterator global_begin() {
if (!root())
return parent().global_begin();
return iterator(left_most());
}
/// return global end iterator
iterator global_end() {
if (!root())
return parent().global_end();
return iterator(right_most(),1);
}
/// const_iterator
class const_iterator {
public:
///@name type definitions for iterator traits
//@{
/// define the iterator category
typedef std::bidirectional_iterator_tag iterator_category;
/// define the difference type
typedef int difference_type;
/// define the value type
typedef const Value value_type;
/// define the reference type
typedef const value_type& reference;
/// define the pointer type
typedef const value_type * pointer;
//@}
public:
///@name constructors
//@{
/// default constructor
const_iterator() : pointee(0), post_end(0), pre_begin(0) { }
/// constructor taking pointee
- const_iterator(binary_tree * p, std::size_t end = 0)
+ const_iterator(const binary_tree * p, std::size_t end = 0)
: pointee(p), post_end(end), pre_begin(0) { }
/// conversion from iterator
const_iterator(const iterator& x)
: pointee(x.pointee), post_end(x.post_end), pre_begin(x.pre_begin) { }
//@}
public:
///@name comparisons
//@{
/// comparison
bool operator==(const const_iterator& x) const {
return ((pointee == x.pointee) &&
(post_end == x.post_end) &&
(pre_begin == x.pre_begin));
}
/// comparison
bool operator!=(const const_iterator& x) const { return !(*this == x); }
//@}
public:
///@name dereference and indirection
//@{
/// dereference
reference operator*() const { return pointee->value(); }
/// indirection
pointer operator->() const { return &**this; }
/// return reference to the node
const binary_tree& node() const { return *pointee; }
//@}
///@name biderectional iterator increment/decrements
//@{
/// pre-increment
const_iterator& operator++() {
if (post_end) { ++post_end; return *this; }
if (pre_begin) { --pre_begin; return *this; }
if(!(pointee->right_neighbour())) { post_end = 1; return *this; }
pointee = pointee->right_neighbour();
return *this;
}
/// pre-decrement
const_iterator& operator--() {
if (post_end) { --post_end; return *this; }
if (pre_begin) { ++pre_begin; return *this; }
if(!(pointee->left_neighbour())) { pre_begin = 1; return *this; }
pointee = pointee->left_neighbour();
return *this;
}
/// post-increment
const_iterator operator++(int) {
const_iterator tmp = *this;
++(*this);
return tmp;
}
/// post-decrement
const_iterator operator--(int) {
const_iterator tmp = *this;
--(*this);
return tmp;
}
//@}
private:
/// the node pointed to
- binary_tree * pointee;
+ const binary_tree * pointee;
/// the distance from --end() (if above --end())
std::size_t post_end;
/// the distance from begin() (if below begin())
std::size_t pre_begin;
};
/// return begin const_iterator
const_iterator begin() const { return const_iterator(left_most()); }
/// return end const_iterator
const_iterator end() const { return const_iterator(right_most(),1); }
/// return global begin iterator
const_iterator global_begin() const {
if (!root())
return parent().global_begin();
return iterator(left_most());
}
/// return global end iterator
const_iterator global_end() const {
if (!root())
return parent().global_end();
return iterator(right_most(),1);
}
private:
/// set the left neighbour
void left_neighbour(binary_tree * n) { neighbours_.first = n; }
/// set the right neighbour
void right_neighbour(binary_tree * n) { neighbours_.second = n; }
/// get the left neighbour
binary_tree * left_neighbour() const { return neighbours_.first; }
/// get the right neighbour
binary_tree * right_neighbour() const { return neighbours_.second; }
/// return the left-most leaf
binary_tree * left_most() {
if(leaf()) return this;
return left_child().left_most();
}
/// return the right-most leaf
binary_tree * right_most() {
if(leaf()) return this;
return right_child().right_most();
}
/// return the left-most leaf
const binary_tree * left_most() const {
if(leaf()) return this;
return left_child().left_most();
}
/// return the right-most leaf
const binary_tree * right_most() const {
if(leaf()) return this;
return right_child().right_most();
}
/// the iterator is a good friend
friend class binary_tree<value_type>::iterator;
/// the iterator is a good friend
friend class binary_tree<value_type>::const_iterator;
/// the left and right neighbours of this node
std::pair<binary_tree*,binary_tree*> neighbours_;
//@}
public:
/// return true, if this node is empty
bool empty() const { return root() && leaf() && !value_; }
/// clear this node
void clear() {
neighbours_ = std::make_pair<binary_tree*,binary_tree*>(0,0);
parent_ = 0;
value_.reset(0);
if (!leaf()) {
left_child().clear();
right_child().clear();
}
children_.first.reset(0);
children_.second.reset(0);
}
public:
/// split this node
std::pair<iterator,iterator> split(std::pair<value_type,value_type> children) {
assert(leaf());
children_.first.reset(new binary_tree(children.first,this));
children_.second.reset(new binary_tree(children.second,this));
children_.first->left_neighbour(neighbours_.first);
children_.first->right_neighbour(children_.second.get());
children_.second->left_neighbour(children_.first.get());
children_.second->right_neighbour(neighbours_.second);
// adjust original neighbours
if(neighbours_.first) {
neighbours_.first->right_neighbour(children_.first.get());
}
if (neighbours_.second) {
neighbours_.second->left_neighbour(children_.second.get());
}
neighbours_.first = 0; neighbours_.second = 0;
return std::make_pair(iterator(children_.first.get()),iterator(children_.second.get()));
}
public:
/// select using a selector
template<class Selector>
iterator select(const Selector& selector) {
if(leaf()) {
bool use = selector.use(value());
if (use) return iterator(this);
return global_end();
}
std::pair<bool,bool> which(selector.use(value(),left_child().value(),right_child().value()));
assert(!which.first || !which.second);
if (!which.first && !which.second) {
return global_end();
}
if (which.first) {
return left_child().select(selector);
}
else {
return right_child().select(selector);
}
return global_end();
}
/// generate a hash value for the sub-tree
/// selected by the given selector object
template<class Selector, unsigned long bits>
void subtree_hash(const Selector& selector, bit_container<bits>& bhash) {
bhash = bit_container<bits>();
unsigned long pos = 0;
do_subtree_hash<Selector,bits>(selector,bhash,pos);
}
/// accumulate values using a binary function
/// and accessor object
template<class Accessor, class BinaryOp>
typename BinaryOp::result_type accumulate(const Accessor& acc,
BinaryOp binary_op) const {
if (!leaf()) {
return
binary_op(left_child().accumulate(acc,binary_op),
right_child().accumulate(acc,binary_op));
}
return acc.get(value(),true);
}
/// accumulate values only from branches
/// matching a Selector
template<class Selector, class Accessor, class BinaryOp>
typename BinaryOp::result_type accumulate(const Selector& selector,
const Accessor& acc,
BinaryOp binary_op) const {
if (!leaf()) {
std::pair<bool,bool> which(selector.use(value(),left_child().value(),right_child().value()));
assert(which.first || which.second);
if (which.first && which.second) {
return
binary_op(left_child().accumulate(selector,acc,binary_op),
right_child().accumulate(selector,acc,binary_op));
} else if (which.first) {
return left_child().accumulate(selector,acc,binary_op);
} else if (which.second) {
return right_child().accumulate(selector,acc,binary_op);
}
}
return acc.get(value(),true);
}
/// accumulate values using a binary function
/// and accessor object, storing intermediate
/// values in nodes
template<class Accessor, class BinaryOp>
typename BinaryOp::result_type tree_accumulate(const Accessor& acc,
BinaryOp binary_op) {
if (!leaf()) {
acc.set(value()) =
binary_op(left_child().tree_accumulate(acc,binary_op),
right_child().tree_accumulate(acc,binary_op));
return acc.get(value(),false);
}
acc.set(value()) = acc.get(value(),true);
return acc.get(value(),true);
}
/// accumulate values only from branches
/// matching a Selector
template<class Selector, class Accessor, class BinaryOp>
typename BinaryOp::result_type tree_accumulate(const Selector& selector,
const Accessor& acc,
BinaryOp binary_op) {
if (!leaf()) {
std::pair<bool,bool> which(selector.use(value(),left_child().value(),right_child().value()));
assert(which.first || which.second);
if (which.first && which.second) {
acc.set(value()) =
binary_op(left_child().tree_accumulate(selector,acc,binary_op),
right_child().tree_accumulate(selector,acc,binary_op));
} else if (which.first) {
acc.set(value()) = left_child().tree_accumulate(selector,acc,binary_op);
} else if (which.second) {
acc.set(value()) = right_child().tree_accumulate(selector,acc,binary_op);
}
return acc.get(value(),false);
}
acc.set(value()) = acc.get(value(),true);
return acc.get(value(),true);
}
/// forward propagate a visitor to all children nodes
template<class Visitor>
void cascade(Visitor visitor) const {
if (leaf()) {
visitor.visit(value());
return;
} else visitor.visit(value(),left_child().value(),right_child().value());
left_child().cascade(visitor);
right_child().cascade(visitor);
}
/// succesively split using a generator
template<class Generator>
void generate(Generator generator) {
if (root())
value_.reset(new value_type(generator.root()));
if (generator.split()) {
std::pair<iterator,iterator> ch = split(generator.generate(value()));
ch.first.node().generate(generator);
ch.second.node().generate(generator);
}
}
public:
///@name Public member access
//@{
/// return the value held by this node
value_type& value() { return *value_; }
/// return the value held by this node
const value_type& value() const { return *value_; }
/// return true, if this is the root node
bool root() const { return !parent_; }
/// return true, if this node has got children
bool leaf() const { return !(children_.first.get() && children_.second.get()); }
//@}
public:
///@name put and get from streams
//@{
/// forward visitor writing out the tree to given ostream
template<class OStream>
struct ostream_visitor {
/// construct from ostream reference
explicit ostream_visitor(OStream& os) : os_(&os), first_time_(true) {}
/// visit a leaf node
void visit(const value_type&) {
(*os_) << "end_branch";
ostream_traits<OStream>::separator(*os_);
}
/// visit a branching
void visit(const value_type& parent,
const value_type& left, const value_type& right) {
if (first_time_) {
(*os_) << "root_node";
ostream_traits<OStream>::separator(*os_);
parent.put(*os_);
first_time_ = false;
}
(*os_) << "left_child";
ostream_traits<OStream>::separator(*os_);
left.put(*os_);
(*os_) << "right_child";
ostream_traits<OStream>::separator(*os_);
right.put(*os_);
}
private:
/// pointer to the ostream to write to
OStream* os_;
/// whether we are at the or not
bool first_time_;
};
/// generator reading binary tree from istream
template<class IStream>
struct istream_generator {
/// construct from istream reference
explicit istream_generator(IStream& is)
: is_(&is), children_(), tag_("") {}
/// copy constructor
istream_generator(const istream_generator& x)
: is_(x.is_), children_(), tag_("") {}
/// read the root node
value_type root() {
*is_ >> tag_;
assert(tag_ == "root_node");
value_type rnode;
rnode.get(*is_);
return rnode;
}
/// read children nodes
bool split() {
*is_ >> tag_;
if (tag_ == "end_branch") {
return false;
}
assert (tag_ == "left_child");
children_.first.get(*is_);
*is_ >> tag_;
assert(tag_ == "right_child");
children_.second.get(*is_);
return true;
}
/// return the children generated
std::pair<value_type,value_type> generate(const value_type&) {
return children_;
}
/// initialize a leaf
void initialize_leaf(const value_type&) {}
private:
/// pointer to the istream used
IStream* is_;
/// the children currently handled
std::pair<value_type,value_type> children_;
/// temporary storage for tags
std::string tag_;
};
/// put to ostream
template<class OStream>
void put(OStream& os) const {
if (empty()) {
os << "empty";
ostream_traits<OStream>::separator(os);
return;
} else if (root() && leaf()) {
os << "root_only";
ostream_traits<OStream>::separator(os);
value().put(os);
return;
} else {
os << "non_empty";
ostream_traits<OStream>::separator(os);
}
assert(root());
cascade(ostream_visitor<OStream>(os));
}
/// get from istream
template<class IStream>
void get(IStream& is) {
std::string state;
is >> state;
if (state == "empty") {
return;
}
if (state == "root_only") {
value_.reset(new value_type());
value().get(is);
return;
}
assert(empty());
generate(istream_generator<IStream>(is));
}
//@}
private:
/// calculate hash value
template<class Selector, unsigned long bits>
void do_subtree_hash(const Selector& selector,
bit_container<bits>& current,
unsigned long& position,
bool selected = true) const {
if (!leaf()) {
std::pair<bool,bool> which(false,false);
if (selected)
which = selector.use(value(),left_child().value(),right_child().value());
current.bit(position,which.first);
current.bit(position+1,which.second);
position += 2;
left_child().do_subtree_hash(selector,current,position,which.first && selected);
right_child().do_subtree_hash(selector,current,position,which.second && selected);
}
}
private:
///@name private member access
//@{
/// return the parent of this node
binary_tree& parent() { assert(parent_); return *parent_; }
/// return the parent of this node
const binary_tree& parent() const { assert(parent_); return *parent_; }
/// return the left child of this node
binary_tree& left_child() { assert(children_.first.get()); return *children_.first; }
/// return the left child of this node
const binary_tree& left_child() const { assert(children_.first.get()); return *children_.first; }
/// return the right child of this node
binary_tree& right_child() { assert(children_.second.get()); return *children_.second; }
/// return the right child of this node
const binary_tree& right_child() const { assert(children_.second.get()); return *children_.second; }
//@}
private:
/// the parent of this node
binary_tree * parent_;
/// the cell held by this node
boost::scoped_ptr<value_type> value_;
/// the children of this node
std::pair<boost::scoped_ptr<binary_tree>,
boost::scoped_ptr<binary_tree> > children_;
};
}
#endif // EXSAMPLE_binary_tree_h_included
diff --git a/Exsample2/exsample/cell.icc b/Exsample2/exsample/cell.icc
--- a/Exsample2/exsample/cell.icc
+++ b/Exsample2/exsample/cell.icc
@@ -1,409 +1,430 @@
// -*- C++ -*-
//
// cell.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
//
// ExSample is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
namespace exsample {
template<class Random>
void cell_info::select (Random& rnd_gen,
std::vector<double>& p) {
std::transform(lower_left_.begin(),lower_left_.end(),
upper_right_.begin(),p.begin(),
rnd_gen);
++attempted_;
}
template<class Random>
void cell_info::select (Random& rnd_gen,
std::vector<double>& p,
const std::vector<bool>& sample) {
conditional_transform(lower_left_.begin(),lower_left_.end(),
upper_right_.begin(),sample.begin(),
p.begin(),rnd_gen);
++attempted_;
}
template<class Random, class Function, class SlaveStatistics>
void cell_info::explore(Random& rnd_gen,
const adaption_info& ainfo,
Function* function, statistics* stats,
SlaveStatistics& opt) {
function->start_presampling();
unsigned long n_sampled = 0;
std::vector<double> ll = lower_left_;
std::vector<double> ur = upper_right_;
double val = 0.;
std::vector<double> pos (ll.size());
while (n_sampled < ainfo.presampling_points) {
std::transform(ll.begin(),ll.end(),
ur.begin(),pos.begin(),
rnd_gen);
val = function->evaluate(pos);
stats->presampled(val);
opt.select(val);
if (std::abs(val) > std::abs(overestimate_)) {
overestimate_ = std::abs(val);
last_max_position_ = pos;
}
selected(pos,std::abs(val),ainfo);
++n_sampled;
}
function->stop_presampling();
}
template<class Random, class Function>
void cell_info::explore (Random& rnd_gen,
const adaption_info& ainfo, Function* function) {
function->start_presampling();
unsigned long n_sampled = 0;
std::vector<double> ll = lower_left_;
std::vector<double> ur = upper_right_;
double val = 0.;
std::vector<double> pos (ll.size());
std::vector<double> next_to_largest (ll.size());
while (n_sampled < ainfo.presampling_points) {
std::transform(ll.begin(),ll.end(),
ur.begin(),pos.begin(),
rnd_gen);
val = function->evaluate(pos);
if (std::abs(val) > std::abs(overestimate_)) {
next_to_largest = last_max_position_;
overestimate_ = std::abs(val);
last_max_position_ = pos;
}
++n_sampled;
}
function->stop_presampling();
}
template<class OStream>
void cell_info::put (OStream& os) const {
os << overestimate_;
ostream_traits<OStream>::separator(os);
os << volume_;
ostream_traits<OStream>::separator(os);
os << lower_left_.size();
ostream_traits<OStream>::separator(os);
for (std::size_t k = 0; k < lower_left_.size(); ++k) {
os << lower_left_[k];
ostream_traits<OStream>::separator(os);
}
for (std::size_t k = 0; k < upper_right_.size(); ++k) {
os << upper_right_[k];
ostream_traits<OStream>::separator(os);
}
for (std::size_t k = 0; k < mid_point_.size(); ++k) {
os << mid_point_[k];
ostream_traits<OStream>::separator(os);
}
for (std::size_t k = 0; k < last_max_position_.size(); ++k) {
os << last_max_position_[k];
ostream_traits<OStream>::separator(os);
}
for (std::size_t k = 0; k < avg_weight_.size(); ++k) {
os << avg_weight_[k].first;
ostream_traits<OStream>::separator(os);
os << avg_weight_[k].second;
ostream_traits<OStream>::separator(os);
}
+ os << attempted_;
+ ostream_traits<OStream>::separator(os);
+ os << accepted_;
+ ostream_traits<OStream>::separator(os);
+ os << parametric_missing_map_.size();
+ ostream_traits<OStream>::separator(os);
+ for ( std::map<bit_container<parameter_hash_bits>,int>::const_iterator p =
+ parametric_missing_map_.begin(); p != parametric_missing_map_.end(); ++p ) {
+ p->first.put(os);
+ os << p->second;
+ ostream_traits<OStream>::separator(os);
+ }
}
template<class IStream>
void cell_info::get (IStream& is) {
std::size_t dim;
is >> overestimate_ >> volume_ >> dim;
lower_left_.resize(dim);
for (std::size_t k = 0; k < lower_left_.size(); ++k) {
is >> lower_left_[k];
}
upper_right_.resize(dim);
for (std::size_t k = 0; k < upper_right_.size(); ++k) {
is >> upper_right_[k];
}
mid_point_.resize(dim);
for (std::size_t k = 0; k < mid_point_.size(); ++k) {
is >> mid_point_[k];
}
last_max_position_.resize(dim);
for (std::size_t k = 0; k < last_max_position_.size(); ++k) {
is >> last_max_position_[k];
}
avg_weight_.resize(dim);
for (std::size_t k = 0; k < avg_weight_.size(); ++k) {
is >> avg_weight_[k].first >> avg_weight_[k].second;
}
+ is >> attempted_ >> accepted_ >> dim;
+ for ( size_t k = 0; k < dim; ++k ) {
+ bit_container<parameter_hash_bits> in;
+ in.get(is);
+ is >> parametric_missing_map_[in];
+ }
}
template<class Random, class Function>
std::pair<cell,cell >
cell::split (std::pair<std::size_t,double> split_d,
Random& rnd_gen,
Function* function,
const adaption_info& ainfo,
const std::vector<bool>& sampled) {
assert(!missing_events() && !info().parametric_compensating());
split_dimension_ = split_d.first;
split_point_ = split_d.second;
std::vector<double> lower_left1 = info().lower_left();
std::vector<double> upper_right1 = info().upper_right();
std::vector<double> lower_left2 = info().lower_left();
std::vector<double> upper_right2 = info().upper_right();
upper_right1[split_dimension_] = split_point_;
lower_left2[split_dimension_] = split_point_;
std::pair<cell,cell> children;
if (sampled.empty())
children = std::pair<cell,cell>(cell(lower_left1,upper_right1,ainfo),
cell(lower_left2,upper_right2,ainfo));
else
children = std::pair<cell,cell> (cell(lower_left1,upper_right1,sampled,ainfo),
cell(lower_left2,upper_right2,sampled,ainfo));
if (info().last_max_position()[split_dimension_] <= split_point_) {
children.first.info().overestimate(info().overestimate(),info().last_max_position());
children.second.info().explore(rnd_gen,ainfo,function);
} else {
children.second.info().overestimate(info().overestimate(),info().last_max_position());
children.first.info().explore(rnd_gen,ainfo,function);
}
cell_info_.reset(0);
children.first.integral(children.first.info().volume() * children.first.info().overestimate());
children.second.integral(children.second.info().volume() * children.second.info().overestimate());
return children;
}
template<class OStream>
void cell::put (OStream& os) const {
os << split_dimension_;
ostream_traits<OStream>::separator(os);
os << split_point_;
ostream_traits<OStream>::separator(os);
os << integral_;
ostream_traits<OStream>::separator(os);
+ os << missing_events_;
+ ostream_traits<OStream>::separator(os);
if (cell_info_) {
os << "has_cell_info";
ostream_traits<OStream>::separator(os);
cell_info_->put(os);
} else {
os << "has_no_cell_info";
ostream_traits<OStream>::separator(os);
}
}
template<class IStream>
void cell::get (IStream& is) {
std::string info_tag;
is >> split_dimension_ >> split_point_
- >> integral_ >> info_tag;
+ >> integral_ >> missing_events_
+ >> info_tag;
if (info_tag == "has_cell_info") {
cell_info_.reset(new cell_info());
cell_info_->get(is);
}
}
inline cell_info::cell_info()
: overestimate_(0.), volume_(0.),
lower_left_(), upper_right_(), mid_point_(),
last_max_position_(), avg_weight_(),
attempted_(0), accepted_(0) {}
inline cell_info::cell_info(const std::vector<double>& ll,
const std::vector<double>& ur,
const adaption_info& ainfo)
: overestimate_(0.), volume_(),
lower_left_(ll), upper_right_(ur), mid_point_(),
last_max_position_(),
avg_weight_(std::vector<std::pair<double,double> >
(ainfo.dimension,std::make_pair(0.,0.))),
attempted_(0), accepted_(0) {
std::vector<double> delta;
std::transform(ur.begin(),ur.end(),
ll.begin(),std::back_inserter(delta),
std::minus<double>());
volume_ =
std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
std::transform(ur.begin(),ur.end(),
ll.begin(),std::back_inserter(mid_point_),
std::plus<double>());
for (std::size_t k = 0; k < ainfo.dimension; ++k)
mid_point_[k] /= 2.;
}
inline cell_info::cell_info(const std::vector<double>& ll,
const std::vector<double>& ur,
const std::vector<bool>& sampled_variables,
const adaption_info& ainfo)
: overestimate_(0.), volume_(),
lower_left_(ll), upper_right_(ur), mid_point_(),
last_max_position_(),
avg_weight_(std::vector<std::pair<double,double> >
(ainfo.dimension,std::make_pair(0.,0.))),
attempted_(0), accepted_(0) {
std::vector<double> delta;
conditional_transform(ur.begin(),ur.end(),
ll.begin(),sampled_variables.begin(),
std::back_inserter(delta),
std::minus<double>());
volume_ =
std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
std::transform(ur.begin(),ur.end(),
ll.begin(),std::back_inserter(mid_point_),
std::plus<double>());
for (std::size_t k = 0; k < ainfo.dimension; ++k)
mid_point_[k] /= 2.;
}
inline int cell_info::parametric_missing(const bit_container<parameter_hash_bits>& id) const {
std::map<bit_container<parameter_hash_bits>,int>::const_iterator mit
= parametric_missing_map_.find(id);
if (mit == parametric_missing_map_.end())
return 0;
return mit->second;
}
inline void cell_info::parametric_missing(const bit_container<parameter_hash_bits>& id, int n) {
if (n == 0) {
std::map<bit_container<parameter_hash_bits>,int>::iterator mit
= parametric_missing_map_.find(id);
if (mit != parametric_missing_map_.end())
parametric_missing_map_.erase(mit);
return;
}
parametric_missing_map_[id] = n;
}
inline void cell_info::increase_parametric_missing(const bit_container<parameter_hash_bits>& id) {
std::map<bit_container<parameter_hash_bits>,int>::iterator mit
= parametric_missing_map_.find(id);
if (mit != parametric_missing_map_.end()) {
mit->second += 1;
if (mit->second == 0) parametric_missing_map_.erase(mit);
} else parametric_missing_map_[id] = 1;
}
inline void cell_info::decrease_parametric_missing(const bit_container<parameter_hash_bits>& id) {
std::map<bit_container<parameter_hash_bits>,int>::iterator mit
= parametric_missing_map_.find(id);
if (mit != parametric_missing_map_.end()) {
mit->second -= 1;
if (mit->second == 0) parametric_missing_map_.erase(mit);
} else assert(false);
}
inline void cell_info::selected(const std::vector<double>& p,
double weight,
const adaption_info& ainfo) {
for (std::size_t k = 0; k < p.size(); ++k) {
if (ainfo.adapt[k]) {
if (p[k] < mid_point_[k])
avg_weight_[k].first += weight;
else
avg_weight_[k].second += weight;
}
}
}
inline std::pair<std::size_t,double> cell_info::get_split (const adaption_info& ainfo,
bool& worth) const {
std::size_t split_d = 0;
double gain = 0.;
for (std::size_t k = 0; k < ainfo.dimension; ++k) {
double xgain = 0.;
double left = avg_weight_[k].first;
double right = avg_weight_[k].second;
if (left+right > 0.) {
xgain = std::abs(left-right)/(left+right);
}
if (xgain > gain) {
gain = xgain;
split_d = k;
}
}
worth = (gain >= ainfo.gain_threshold);
return std::make_pair(split_d,mid_point_[split_d]);
}
inline bool cell_info::contains_parameter (const std::vector<double>& point,
const std::vector<bool>& sampled) const {
std::vector<double>::const_iterator p = point.begin();
std::vector<double>::const_iterator l = lower_left_.begin();
std::vector<double>::const_iterator u = upper_right_.begin();
std::vector<bool>::const_iterator f = sampled.begin();
for (; p < point.end(); ++p, ++f, ++l, ++u)
if (!(*f)) {
if (((*l) > (*p)) ||
((*u) < (*p)))
return false;
}
return true;
}
inline cell::cell()
: split_dimension_(0), split_point_(0.),
integral_(0.), missing_events_(0),
cell_info_(0) {}
inline cell::cell(const std::vector<double>& ll,
const std::vector<double>& ur,
const adaption_info& ainfo)
: split_dimension_(0), split_point_(0.),
integral_(0.), missing_events_(0),
cell_info_(new cell_info(ll,ur,ainfo)) {}
inline cell::cell(const std::vector<double>& ll,
const std::vector<double>& ur,
const std::vector<bool>& sampled_variables,
const adaption_info& ainfo)
: split_dimension_(0), split_point_(0.),
integral_(0.), missing_events_(0),
cell_info_(new cell_info(ll,ur,sampled_variables,ainfo)) {}
inline cell::cell(const cell& x)
: split_dimension_(x.split_dimension_),
split_point_(x.split_point_),
integral_(x.integral_),
missing_events_(x.missing_events_),
cell_info_(0) {
if (x.cell_info_)
cell_info_.reset(new cell_info(*x.cell_info_));
}
inline cell& cell::operator=(const cell& x) {
if (this == &x)
return *this;
split_dimension_ = x.split_dimension_;
split_point_ = x.split_point_;
integral_ = x.integral_;
missing_events_ = x.missing_events_;
if (x.cell_info_)
cell_info_.reset(new cell_info(*x.cell_info_));
return *this;
}
}
diff --git a/Exsample2/exsample/exponential_generator.h b/Exsample2/exsample/exponential_generator.h
--- a/Exsample2/exsample/exponential_generator.h
+++ b/Exsample2/exsample/exponential_generator.h
@@ -1,210 +1,220 @@
// -*- C++ -*-
//
// exponential_generator.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
//
// ExSample is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
#ifndef EXSAMPLE_exponential_generator_h_included
#define EXSAMPLE_exponential_generator_h_included
#include "cell.h"
#include "selectors.h"
#include "statistics.h"
#include "linear_interpolator.h"
#include "binary_tree.h"
namespace exsample {
/// \brief Exception thrown, if the exponential_generator has just changed its
/// state. The attempt of generating an event should be repeated.
struct exponential_regenerate{};
/// \brief The generator for sudakov-type distributions.
template<class Function, class Random>
class exponential_generator {
public:
/// default constructor
exponential_generator()
: function_(0), check_events_(0), adaption_info_(), root_cell_(),
rnd_gen_(), did_split_(false), initialized_(false),
evolution_variable_(0), evolution_cutoff_(0.),
sample_variables_(), sample_other_variables_(),
parameter_splits_(),
last_cell_(), last_point_(), last_value_(0.),
last_parameter_bin_(), exponents_(),
last_exponent_integrand_(),
last_exponent_(), compensating_(false),
integral_accessor_(), missing_accessor_(),
parametric_selector_(), exponent_selector_(),
parametric_sampler_(), attempts_(0), accepts_(0),
splits_(0) {}
public:
/// initialize this generator
void initialize();
/// finalize this generator
void finalize() {}
/// generate an event, returning
/// the sign of the weight or zero
/// for an event below the evolution cutoff
double generate();
/// return the last sampled phase space point
const std::vector<double>& last_point() const { return last_point_; }
/// return the last evaluated function
double last_value() const { return last_value_; }
/// indicate that the last generated point has been rejected
void reject() {
last_cell_->info().reject();
}
public:
/// return true, if this generator has been initialized
bool initialized() const { return initialized_; }
/// return true, if at least one split has been performed
bool did_split() const { return did_split_; }
/// access the function
Function& function() { return *function_; }
/// set the function
void function(Function * f) { function_ = f; }
/// access the adaption_info object
adaption_info& sampling_parameters() { return adaption_info_; }
+ public:
+
+ /// put to ostream
+ template<class OStream>
+ void put(OStream& os) const;
+
+ /// get from istream
+ template<class IStream>
+ void get(IStream& is);
+
private:
/// check for and possibly split
/// the last selected cell
bool split();
/// get the projection of the density integrating over every
/// variable to be sampled, except the evolution variable for the
/// indicated parameter point. the k'th entry in
/// last_exponent_integrand_ is the value in the evolution
/// variable bin from evolution_splits_[k] to
/// evolution_splits_[k+1]
void get_exponent();
/// compensate
void compensate();
/// get all parameter points to build
/// all possible sub tree hashes
std::set<std::vector<double> > parameter_points();
/// get all parameter points to build
/// all possible sub tree hashes
void recursive_parameter_points(std::set<std::vector<double> >&,
std::vector<double>&,
size_t);
/// function to be sampled
Function * function_;
/// the number of events after which
/// a cell is checked for splits
unsigned long check_events_;
/// the adaption info object
adaption_info adaption_info_;
/// the root cell
binary_tree<cell> root_cell_;
/// the random number generator to be used
rnd_generator<Random> rnd_gen_;
/// wether a split has already been performed
bool did_split_;
/// wether this generator has been initialized
bool initialized_;
/// the position of the evolution variable
std::size_t evolution_variable_;
/// the cutoff on the evolution variable
double evolution_cutoff_;
/// flags of variables to be sampled
/// including the evolution variable
std::vector<bool> sample_variables_;
/// flags of variables to be sampled
/// excluding the evolution variable
std::vector<bool> sample_other_variables_;
/// the splits in any parameter done so far
/// (including the evolution variable)
std::map<std::size_t,std::vector<double> > parameter_splits_;
/// the last selected cell
binary_tree<cell>::iterator last_cell_;
/// the last sampled phasespace point
std::vector<double> last_point_;
/// the last function value
double last_value_;
/// the last parameter bin id
bit_container<parameter_hash_bits> last_parameter_bin_;
/// map parameter bin ids to exponent interpolations
std::map<bit_container<parameter_hash_bits>,linear_interpolator > exponents_;
/// the last exponent integrand
std::vector<double> last_exponent_integrand_;
/// the last exponent
std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator last_exponent_;
/// wether or not we are compensating
bool compensating_;
/// the integral accessor to be used
integral_accessor integral_accessor_;
/// the missing events accessor to be used
parametric_missing_accessor missing_accessor_;
/// the parametric selector to be used
parametric_selector parametric_selector_;
/// the parametric selector to be used for parameter bins
parametric_selector exponent_selector_;
/// the parametric sampler to be used
parametric_sampling_selector<rnd_generator<Random> > parametric_sampler_;
/// the number of trials in the veto loo so far
unsigned long attempts_;
/// the number of accepted events so far
unsigned long accepts_;
/// number of splits done
unsigned long splits_;
};
}
#include "exponential_generator.icc"
#endif // EXSAMPLE_exponential_generator_h_included
diff --git a/Exsample2/exsample/exponential_generator.icc b/Exsample2/exsample/exponential_generator.icc
--- a/Exsample2/exsample/exponential_generator.icc
+++ b/Exsample2/exsample/exponential_generator.icc
@@ -1,324 +1,377 @@
// -*- C++ -*-
//
// exponential_generator.icc is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
//
// ExSample is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
namespace exsample {
template<class Function, class Random>
void exponential_generator<Function,Random>::initialize() {
adaption_info_.dimension = function_->dimension();
adaption_info_.lower_left = function_->support().first;
adaption_info_.upper_right = function_->support().second;
if (adaption_info_.adapt.empty())
adaption_info_.adapt = std::vector<bool>(adaption_info_.dimension,true);
evolution_variable_ = function_->evolution_variable();
evolution_cutoff_ = function_->evolution_cutoff();
sample_variables_ = function_->variable_flags();
sample_other_variables_ = sample_variables_;
sample_other_variables_[evolution_variable_] = false;
last_point_.resize(adaption_info_.dimension);
parametric_selector_ = parametric_selector(&last_point_,sample_other_variables_);
exponent_selector_ = parametric_selector(&last_point_,sample_variables_);
missing_accessor_ = parametric_missing_accessor(&last_parameter_bin_);
parametric_sampler_ = parametric_sampling_selector<rnd_generator<Random> >
(&last_point_,&last_parameter_bin_,sample_other_variables_,rnd_gen_);
if (initialized_) return;
splits_ = 0;
for ( std::size_t k = 0; k < adaption_info_.dimension; ++k ) {
if ( sample_other_variables_[k] )
continue;
parameter_splits_[k].push_back(adaption_info_.lower_left[k]);
parameter_splits_[k].push_back(adaption_info_.upper_right[k]);
}
root_cell_ =
binary_tree<cell>(cell(adaption_info_.lower_left,
adaption_info_.upper_right,
sample_other_variables_,
adaption_info_));
root_cell_.value().info().explore(rnd_gen_,adaption_info_,function_);
root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
last_exponent_integrand_.resize(1);
check_events_ = adaption_info_.presampling_points;
initialized_ = true;
}
template<class Function, class Random>
bool exponential_generator<Function,Random>::split () {
if (adaption_info_.freeze_grid == accepts_)
return false;
if (compensating_)
return false;
if (!(*last_cell_).info().bad(adaption_info_)) return false;
bool dosplit = false;
std::pair<std::size_t,double> sp =
(*last_cell_).info().get_split(adaption_info_,dosplit);
if (!dosplit) return false;
if (!adaption_info_.adapt[sp.first]) return false;
if (splits_ == parameter_hash_bits/2)
return false;
++splits_;
last_cell_.node().split((*last_cell_).split(sp,rnd_gen_,function_,adaption_info_,sample_other_variables_));
if ( !sample_other_variables_[sp.first] ) {
if ( std::find(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end(),sp.second)
== parameter_splits_[sp.first].end() ) {
parameter_splits_[sp.first].push_back(sp.second);
std::sort(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end());
if ( sp.first == evolution_variable_ ) {
last_exponent_integrand_.push_back(0.);
}
}
}
did_split_ = true;
last_point_ = function_->parameter_point();
root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
exponents_.clear();
get_exponent();
return true;
}
template<class Function, class Random>
void exponential_generator<Function,Random>::get_exponent () {
last_parameter_bin_.reset();
root_cell_.subtree_hash (exponent_selector_,last_parameter_bin_);
last_exponent_ = exponents_.find(last_parameter_bin_);
if (last_exponent_ != exponents_.end())
return;
exponents_[last_parameter_bin_] = linear_interpolator();
last_exponent_ = exponents_.find(last_parameter_bin_);
double old_evo = last_point_[evolution_variable_];
std::vector<double>::iterator exp_it = last_exponent_integrand_.begin();
for (std::vector<double>::iterator esp = parameter_splits_[evolution_variable_].begin();
esp < prior(parameter_splits_[evolution_variable_].end()); ++esp, ++exp_it) {
last_point_[evolution_variable_] = (*esp + *next(esp))/2.;
*exp_it = root_cell_.accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
}
exp_it = prior(last_exponent_integrand_.end());
double total = 0.;
for (std::vector<double>::iterator esp = prior(parameter_splits_[evolution_variable_].end());
esp > parameter_splits_[evolution_variable_].begin(); --esp, --exp_it) {
last_exponent_->second.set_interpolation(*esp,total);
total += (*exp_it) * ((*esp) - (*prior(esp)));
}
last_exponent_->second.set_interpolation(parameter_splits_[evolution_variable_].front(),total);
last_point_[evolution_variable_] = old_evo;
}
template<class Function, class Random>
std::set<std::vector<double> >
exponential_generator<Function,Random>::parameter_points() {
std::set<std::vector<double> > res;
std::vector<double> pt(adaption_info_.dimension,0.);
recursive_parameter_points(res,pt,0);
return res;
}
template<class Function, class Random>
void exponential_generator<Function,Random>::
recursive_parameter_points(std::set<std::vector<double> >& res,
std::vector<double>& pt,
size_t current) {
if ( current == adaption_info_.dimension ) {
res.insert(pt);
return;
}
if ( sample_variables_[current] ) {
recursive_parameter_points(res,pt,current+1);
return;
}
for ( std::vector<double>::const_iterator sp =
parameter_splits_[current].begin();
sp != prior(parameter_splits_[current].end()); ++sp ) {
pt[current] = (*sp + *next(sp))/2.;
recursive_parameter_points(res,pt,current+1);
}
}
template<class Function, class Random>
void exponential_generator<Function,Random>::compensate() {
if (!did_split_) {
exponents_.clear();
root_cell_.value().info().overestimate(last_value_,last_point_);
root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
last_point_ = function_->parameter_point();
get_exponent();
return;
}
std::vector<double> themaxpoint = last_point_;
std::set<std::vector<double> > id_points
= parameter_points();
for ( std::set<std::vector<double> >::const_iterator id =
id_points.begin(); id != id_points.end(); ++id ) {
last_point_ = *id;
get_exponent();
}
std::map<bit_container<parameter_hash_bits>,linear_interpolator >
old_exponents = exponents_;
double old_oe = last_cell_->info().overestimate();
last_cell_->info().overestimate(last_value_,themaxpoint);
last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
exponents_.clear();
for ( std::set<std::vector<double> >::const_iterator id =
id_points.begin(); id != id_points.end(); ++id ) {
last_point_ = *id;
get_exponent();
std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
old_exp = old_exponents.find(last_parameter_bin_);
std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
new_exp = exponents_.find(last_parameter_bin_);
assert(old_exp != old_exponents.end() && new_exp != exponents_.end());
double old_norm = 1. - std::exp(-(old_exp->second)(adaption_info_.lower_left[evolution_variable_]));
double new_norm = 1. - std::exp(-(new_exp->second)(adaption_info_.lower_left[evolution_variable_]));
for (binary_tree<cell>::iterator it = root_cell_.begin();
it != root_cell_.end(); ++it) {
if ( !it->info().contains_parameter(last_point_,sample_other_variables_) )
continue;
double old_int = 0.;
double new_int = 0.;
for ( std::vector<double>::const_iterator sp = parameter_splits_[evolution_variable_].begin();
sp != prior(parameter_splits_[evolution_variable_].end()); ++sp ) {
if ( *sp >= it->info().lower_left()[evolution_variable_] &&
*sp < it->info().upper_right()[evolution_variable_] ) {
double xl = *sp;
double xxl = *next(sp);
double old_al =
(old_exp->second.interpolation()[xxl] - old_exp->second.interpolation()[xl]) /
(xxl-xl);
double old_bl =
(xxl * old_exp->second.interpolation()[xl] -
xl * old_exp->second.interpolation()[xxl]) /
(xxl-xl);
double new_al =
(new_exp->second.interpolation()[xxl] - new_exp->second.interpolation()[xl]) /
(xxl-xl);
double new_bl =
(xxl * new_exp->second.interpolation()[xl] -
xl * new_exp->second.interpolation()[xxl]) /
(xxl-xl);
if ( std::abs(old_al) > std::numeric_limits<double>::epsilon() ) {
old_int += (exp(-(old_al*xl+old_bl)) - exp(-(old_al*xxl+old_bl)))/old_al;
} else {
old_int += (xxl-xl)*exp(-old_bl);
}
if ( std::abs(new_al) > std::numeric_limits<double>::epsilon() ) {
new_int += (exp(-(new_al*xl+new_bl)) - exp(-(new_al*xxl+new_bl)))/new_al;
} else {
new_int += (xxl-xl)*exp(-new_bl);
}
}
}
double scaling;
if (it != last_cell_) {
if (old_int > std::numeric_limits<double>::epsilon() &&
new_int > std::numeric_limits<double>::epsilon())
scaling = ((old_norm * new_int) /
(new_norm * old_int)) - 1.;
else
scaling = 0.;
} else {
if (old_int > std::numeric_limits<double>::epsilon() &&
new_int > std::numeric_limits<double>::epsilon())
scaling = ((last_value_ * old_norm * new_int) /
(old_oe * new_norm * old_int)) - 1.;
else
scaling = 0.;
}
it->info().parametric_missing(last_parameter_bin_,
it->info().parametric_missing(last_parameter_bin_) +
static_cast<int>(round(scaling * it->info().attempted())));
if (it->info().parametric_missing(last_parameter_bin_) != 0) {
compensating_ = true;
}
}
}
last_point_ = function_->parameter_point();
}
template<class Function, class Random>
double exponential_generator<Function,Random>::generate () {
if (compensating_) {
compensating_ = false;
for (binary_tree<cell>::iterator it = root_cell_.begin();
it != root_cell_.end(); ++it)
if (it->info().parametric_compensating()) {
compensating_ = true;
break;
}
parametric_sampler_.compensate(compensating_);
}
last_point_ = function_->parameter_point();
if (last_point_[evolution_variable_] < evolution_cutoff_) {
return 0.;
}
unsigned long n_hit_miss = 0;
unsigned long n_select = 0;
double minus_log_r;
root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
get_exponent();
while (true) {
n_select = 0;
minus_log_r = -std::log(rnd_gen_()) + last_exponent_->second(last_point_[evolution_variable_]);
if (!last_exponent_->second.invertible(minus_log_r)) {
return 0.;
}
try {
last_point_[evolution_variable_] = last_exponent_->second.unique_inverse(minus_log_r);
} catch (constant_interpolation& c) {
last_point_[evolution_variable_] = rnd_gen_(c.range.first,c.range.second);
}
assert(!std::isnan(last_point_[evolution_variable_]) &&
!std::isinf(last_point_[evolution_variable_]));
if (last_point_[evolution_variable_] < evolution_cutoff_) {
return 0.;
}
++attempts_;
if (compensating_) {
root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
}
if (parameter_splits_[evolution_variable_].size() > 2)
root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
if (did_split_)
while ((last_cell_ = root_cell_.select(parametric_sampler_)) == root_cell_.end()) {
root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
if(++n_select > adaption_info_.maxtry)
throw selection_maxtry();
}
else
last_cell_ = root_cell_.begin();
last_cell_->info().select(rnd_gen_,last_point_,sample_other_variables_);
last_value_ = function_->evaluate(last_point_);
assert(last_value_ >= 0.);
last_cell_->info().selected(last_point_,last_value_,adaption_info_);
if (last_value_ > last_cell_->info().overestimate()) {
if ( std::abs(last_value_)/last_cell_->info().overestimate() > 2. ) {
last_value_ =
last_value_*(1.+exp(2.*(2.-std::abs(last_value_)/last_cell_->info().overestimate())));
}
compensate();
throw exponential_regenerate();
}
if (last_cell_->info().attempted() % check_events_ == 0) {
if (split()) {
throw exponential_regenerate();
}
}
if (last_value_/last_cell_->info().overestimate() > rnd_gen_())
break;
if(++n_hit_miss > adaption_info_.maxtry)
throw hit_and_miss_maxtry();
}
if (last_value_ == 0.)
return 0.;
++accepts_;
++check_events_;
last_cell_->info().accept();
return 1.;
}
+ template<class Function, class Random>
+ template<class OStream>
+ void exponential_generator<Function,Random>::put (OStream& os) const {
+ os << check_events_; ostream_traits<OStream>::separator(os);
+ adaption_info_.put(os);
+ root_cell_.put(os);
+ os << did_split_; ostream_traits<OStream>::separator(os);
+ os << initialized_; ostream_traits<OStream>::separator(os);
+ os << evolution_variable_; ostream_traits<OStream>::separator(os);
+ os << evolution_cutoff_; ostream_traits<OStream>::separator(os);
+ os << sample_variables_; ostream_traits<OStream>::separator(os);
+ os << sample_other_variables_; ostream_traits<OStream>::separator(os);
+ os << parameter_splits_; ostream_traits<OStream>::separator(os);
+ // last_cell_ is selected new so we ignore it here
+ os << last_point_; ostream_traits<OStream>::separator(os);
+ os << last_value_; ostream_traits<OStream>::separator(os);
+ last_parameter_bin_.put(os);
+ os << exponents_.size(); ostream_traits<OStream>::separator(os);
+ for ( std::map<bit_container<parameter_hash_bits>,linear_interpolator >::const_iterator
+ ex = exponents_.begin(); ex != exponents_.end() ; ++ex ) {
+ ex->first.put(os);
+ ex->second.put(os);
+ }
+ os << last_exponent_integrand_; ostream_traits<OStream>::separator(os);
+ os << compensating_; ostream_traits<OStream>::separator(os);
+ os << attempts_; ostream_traits<OStream>::separator(os);
+ os << accepts_; ostream_traits<OStream>::separator(os);
+ os << splits_; ostream_traits<OStream>::separator(os);
+ }
+
+ template<class Function, class Random>
+ template<class IStream>
+ void exponential_generator<Function,Random>::get (IStream& is) {
+ is >> check_events_;
+ adaption_info_.get(is);
+ root_cell_.get(is);
+ is >> did_split_ >> initialized_ >> evolution_variable_
+ >> evolution_cutoff_ >> sample_variables_ >> sample_other_variables_
+ >> parameter_splits_;
+ // last_cell_ is selected new so we ignore it here
+ is >> last_point_ >> last_value_;
+ last_parameter_bin_.get(is);
+ size_t dim; is >> dim;
+ for ( size_t k = 0; k < dim ; ++k ) {
+ bit_container<parameter_hash_bits> key;
+ key.get(is);
+ exponents_[key].get(is);
+ }
+ is >> last_exponent_integrand_;
+ last_exponent_ = exponents_.find(last_parameter_bin_);
+ is >> compensating_ >> attempts_ >> accepts_ >> splits_;
+ }
+
}
diff --git a/Exsample2/exsample/utility.h b/Exsample2/exsample/utility.h
--- a/Exsample2/exsample/utility.h
+++ b/Exsample2/exsample/utility.h
@@ -1,322 +1,339 @@
// -*- C++ -*-
//
// utility.h is part of ExSample -- A Library for Sampling Sudakov-Type Distributions
//
// Copyright (C) 2008-2011 Simon Platzer -- simon.plaetzer@desy.de
//
// ExSample is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
#ifndef EXSAMPLE_utility_h_included
#define EXSAMPLE_utility_h_included
#include "config.h"
namespace exsample {
+ /// \brief separate quantities written to an ostream
+ template<class OStream>
+ struct ostream_traits {
+
+ /// put the separator to the ostream
+ static void separator(OStream& os) { os << "\n"; }
+
+ };
+
+#ifdef EXSAMPLE_has_ThePEG
+
+ /// \brief separate quantities written to a ThePEG::PersistentOStream
+ template<>
+ struct ostream_traits<ThePEG::PersistentOStream> {
+
+ /// put the separator to the ostream
+ static void separator(ThePEG::PersistentOStream&) { }
+
+ };
+
+#endif // EXSAMPLE_has_ThePEG
+
/// \brief Compile time conversion of unsigned long to bool
template<unsigned long>
struct static_binary {
enum { value = 1 };
};
/// \brief Compile time conversion of unsigned long to bool
template<>
struct static_binary<0> {
enum { value = 0 };
};
/// \brief Fixed-size, packed vector of bools
template<unsigned long bits>
struct bit_container {
enum {
/// the number of bits contained
n_bits = bits,
/// the number of bits fitting in a unsigned long
uint_bits = CHAR_BIT * sizeof(unsigned long),
/// the number of unsigned long segments needed
n_segments = bits / uint_bits + static_binary<bits % uint_bits>::value
};
/// the default constructor
bit_container() {
for (std::size_t i = 0; i < n_segments; ++i)
segments[i] = 0;
}
/// put all values to false
void reset() {
for (std::size_t i = 0; i < n_segments; ++i)
segments[i] = 0;
}
/// compare for equality
bool operator==(const bit_container& x) const {
for (std::size_t i = 0; i < n_segments; ++i)
if(segments[i] != x.segments[i])
return false;
return true;
}
/// compare for ordering
bool operator<(const bit_container& x) const {
for (std::size_t i = 0; i < n_segments; ++i)
if(segments[i] != x.segments[i])
return (segments[i] < x.segments[i]);
return false;
}
/// set the k'th bit
void bit(unsigned long k, bool value) {
assert(k<n_bits);
if (value)
segments[n_segments-k/uint_bits-1] |= (1ul << (k % uint_bits));
else
segments[n_segments-k/uint_bits-1] &= ~(1ul << (k % uint_bits));
}
/// get the k'th bit
- bool bit(unsigned long k) {
+ bool bit(unsigned long k) const {
assert(k<n_bits);
return (segments[n_segments-k/uint_bits-1] & (1ul << (k % uint_bits)));
}
/// print to ostream
template<class OStream>
void dump(OStream& os) const {
for ( unsigned int k = 0; k < n_segments; ++k )
os << segments[k] << " ";
}
+ /// put to ostream
+ template<class OStream>
+ void put(OStream& os) const {
+ for ( size_t k = 0; k < n_segments; ++k ) {
+ os << segments[k];
+ ostream_traits<OStream>::separator(os);
+ }
+ }
+
+ /// get from istream
+ template<class IStream>
+ void get(IStream& is) {
+ for ( size_t k = 0; k < n_segments; ++k ) {
+ is >> segments[k];
+ }
+ }
+
private:
/// segments needed to keep the hash value
unsigned long segments[n_segments];
};
/// \brief square a number
template<class T>
T sqr(T x) {
return x*x;
}
/// \brief cube a number
template<class T>
T cube(T x) {
return x*x*x;
}
/// \brief Round a floating point value to an integer value of the
/// same type.
template<class T>
T round(T x) {
T f = std::floor(x);
T c = std::ceil(x);
if (x < (f+c)/2.)
return f;
return c;
}
/// \brief Calculate fast powers of two.
inline std::size_t two_to(std::size_t n) {
assert(n <= sizeof(std::size_t)*CHAR_BIT);
return (1 << n);
}
- /// \brief separate quantities written to an ostream
- template<class OStream>
- struct ostream_traits {
-
- /// put the separator to the ostream
- static void separator(OStream& os) { os << "\n"; }
-
- };
-
-#ifdef EXSAMPLE_has_ThePEG
-
- /// \brief separate quantities written to a ThePEG::PersistentOStream
- template<>
- struct ostream_traits<ThePEG::PersistentOStream> {
-
- /// put the separator to the ostream
- static void separator(ThePEG::PersistentOStream&) { }
-
- };
-
-#endif // EXSAMPLE_has_ThePEG
-
/// \brief Fast, zero memory-overhead one-dimensional
/// histogram with 2^n equally spaced bins
template<class Statistics>
struct fast_small_histogram {
/// default constructor
fast_small_histogram()
: depth(0), bins(0) {}
/// copy constructor
fast_small_histogram(const fast_small_histogram& x)
: depth(x.depth), bins(0) {
if (x.bins) {
bins.reset(new Statistics[two_to(depth)]);
for(std::size_t k = 0; k < two_to(depth); ++k)
bins[k] = x.bins[k];
}
}
/// assignment
fast_small_histogram& operator=(const fast_small_histogram& x) {
if (&x == this)
return *this;
depth = x.depth;
bins.reset(0);
if (x.bins) {
bins.reset(new Statistics[two_to(depth)]);
for(std::size_t k = 0; k < two_to(depth); ++k)
bins[k] = x.bins[k];
}
return *this;
}
/// construct from depth d, creating 2^d bins
explicit fast_small_histogram(std::size_t d)
: depth(d), bins(0) {
bins.reset(new Statistics[two_to(d)]);
}
/// return the bin from event belongs to given outer boundaries
Statistics& bin(double lower,
double upper,
double event) {
double thelower = lower;
double theupper = upper;
std::size_t bindex = 0;
std::size_t current_depth = 0;
while (true) {
double cut
= (thelower+theupper)/2.;
if (event < cut) {
theupper = cut;
} else {
thelower = cut;
bindex += two_to(depth-current_depth-1);
}
if(++current_depth == depth)
break;
}
return bins[bindex];
}
/// the depth, defining a histogram of 2^depth bins
std::size_t depth;
/// the contained statistics objects
boost::scoped_array<Statistics> bins;
/// put histogram to an ostream
template<class OStream>
void put(OStream& os) const {
os << depth;
ostream_traits<OStream>::separator(os);
for (std::size_t k = 0; k < two_to(depth); ++k) {
bins[k].put(os);
}
}
/// get histogram from an istream
template<class IStream>
void get(IStream& is) {
is >> depth;
bins.reset(new Statistics[two_to(depth)]);
for(std::size_t k = 0; k < two_to(depth); ++k) {
bins[k].get(is);
}
}
};
/// \brief Generalize the transform algorithm to only apply
/// depending on a range of flags accompanying the input range
template<class FirstInputIterator,
class SecondInputIterator,
class FlagIterator,
class OutputIterator,
class BinaryOperation>
OutputIterator conditional_transform(FirstInputIterator first1,
FirstInputIterator last1,
SecondInputIterator first2,
FlagIterator firstf,
OutputIterator result,
BinaryOperation binary_op) {
for (; first1 != last1; ++first1, ++first2, ++firstf, ++result)
if (*firstf)
*result = binary_op(*first1, *first2);
return result;
}
/// \brief calculate a volume given lower left and upper right
/// corners
inline double volume(const std::vector<double>& lower_left,
const std::vector<double>& upper_right) {
std::vector<double> delta;
std::transform(upper_right.begin(),upper_right.end(),
lower_left.begin(),std::back_inserter(delta),
std::minus<double>());
return
std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
}
/// \brief calculate a volume given lower left and upper right
/// corners, taking into account only part of the dimensions, which
/// are flagged with true in the correspponding random access
/// container
inline double volume(const std::vector<double>& lower_left,
const std::vector<double>& upper_right,
const std::vector<bool>& flags) {
std::vector<double> delta;
conditional_transform(upper_right.begin(),upper_right.end(),
lower_left.begin(),flags.begin(),
std::back_inserter(delta),
std::minus<double>());
return
std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
}
/// \brief Exception thrown if the maximum number of attempts to
/// select a cell has been reached.
struct selection_maxtry{};
/// \brief Exception thrown, if the maximum number of misses has
/// been reached.
struct hit_and_miss_maxtry{};
/// \brief Random generator traits.
template<class Random>
struct rnd_generator {
///Generate uniform random number on [0,1]
double operator()() const {
return Random::rnd();
}
///Generate uniform random number on [0,a]
double operator()(double a) const {
return a*Random::rnd();
}
///Generate uniform random number on [a,b]
double operator()(double a, double b) const {
return (a + (b-a)*Random::rnd());
}
};
}
#endif // EXSAMPLE_utility_h_included
diff --git a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
--- a/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
+++ b/MatrixElement/Matchbox/Base/MatchboxMEBase.cc
@@ -1,1187 +1,1187 @@
// -*- C++ -*-
//
// MatchboxMEBase.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MatchboxMEBase class.
//
#include "MatchboxMEBase.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDF/PDF.h"
#include "ThePEG/PDT/PDT.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/DiagramDrawer.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
MatchboxMEBase::MatchboxMEBase()
: MEBase(),
theOneLoop(false),
theOneLoopNoBorn(false) {}
MatchboxMEBase::~MatchboxMEBase() {}
Ptr<MatchboxFactory>::tcptr MatchboxMEBase::factory() const { return theFactory; }
void MatchboxMEBase::factory(Ptr<MatchboxFactory>::tcptr f) { theFactory = f; }
Ptr<Tree2toNGenerator>::tptr MatchboxMEBase::diagramGenerator() const { return factory()->diagramGenerator(); }
Ptr<ProcessData>::tptr MatchboxMEBase::processData() const { return factory()->processData(); }
unsigned int MatchboxMEBase::getNLight() const { return factory()->nLight(); }
double MatchboxMEBase::factorizationScaleFactor() const { return factory()->factorizationScaleFactor(); }
double MatchboxMEBase::renormalizationScaleFactor() const { return factory()->renormalizationScaleFactor(); }
bool MatchboxMEBase::fixedCouplings() const { return factory()->fixedCouplings(); }
bool MatchboxMEBase::fixedQEDCouplings() const { return factory()->fixedQEDCouplings(); }
bool MatchboxMEBase::checkPoles() const { return factory()->checkPoles(); }
bool MatchboxMEBase::verbose() const { return factory()->verbose(); }
void MatchboxMEBase::getDiagrams() const {
if ( diagramGenerator() && processData() ) {
vector<Ptr<Tree2toNDiagram>::ptr> diags;
for ( vector<PDVector>::const_iterator p = subProcesses().begin();
p != subProcesses().end(); ++p ) {
vector<Ptr<Tree2toNDiagram>::ptr>& res =
processData()->diagramMap()[*p];
if ( res.empty() ) {
res = diagramGenerator()->generate(*p,orderInAlphaS(),orderInAlphaEW());
}
copy(res.begin(),res.end(),back_inserter(diags));
}
if ( diags.empty() )
return;
for ( vector<Ptr<Tree2toNDiagram>::ptr>::iterator d = diags.begin();
d != diags.end(); ++d ) {
add(*d);
}
return;
}
throw Exception()
<< "MatchboxMEBase::getDiagrams() expects a Tree2toNGenerator and ProcessData object.\n"
<< "Please check your setup." << Exception::abortnow;
}
Selector<MEBase::DiagramIndex>
MatchboxMEBase::diagrams(const DiagramVector & diags) const {
if ( phasespace() ) {
return phasespace()->selectDiagrams(diags);
}
throw Exception()
<< "MatchboxMEBase::diagrams() expects a MatchboxPhasespace object.\n"
<< "Please check your setup." << Exception::abortnow;
return Selector<MEBase::DiagramIndex>();
}
Selector<const ColourLines *>
MatchboxMEBase::colourGeometries(tcDiagPtr diag) const {
if ( matchboxAmplitude() ) {
if ( !matchboxAmplitude()->haveColourFlows() )
throw Exception() << "A colour flow implementation is not present."
<< Exception::abortnow;
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
return matchboxAmplitude()->colourGeometries(diag);
}
throw Exception()
<< "MatchboxMEBase::colourGeometries() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return Selector<const ColourLines *>();
}
unsigned int MatchboxMEBase::orderInAlphaS() const {
if ( matchboxAmplitude() ) {
return matchboxAmplitude()->orderInGs();
}
throw Exception()
<< "MatchboxMEBase::orderInAlphaS() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0;
}
unsigned int MatchboxMEBase::orderInAlphaEW() const {
if ( matchboxAmplitude() ) {
return matchboxAmplitude()->orderInGem();
}
throw Exception()
<< "MatchboxMEBase::orderInAlphaEW() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0;
}
void MatchboxMEBase::setXComb(tStdXCombPtr xc) {
MEBase::setXComb(xc);
lastMatchboxXComb(xc);
if ( phasespace() )
phasespace()->setXComb(xc);
if ( scaleChoice() )
scaleChoice()->setXComb(xc);
if ( matchboxAmplitude() )
matchboxAmplitude()->setXComb(xc);
}
double MatchboxMEBase::generateIncomingPartons(const double* r1, const double* r2) {
// shamelessly stolen from PartonExtractor.cc
Energy2 shmax = lastCuts().sHatMax();
Energy2 shmin = lastCuts().sHatMin();
Energy2 sh = shmin*pow(shmax/shmin, *r1);
double ymax = lastCuts().yHatMax();
double ymin = lastCuts().yHatMin();
double km = log(shmax/shmin);
ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/sh)));
ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/sh)));
double y = ymin + (*r2)*(ymax - ymin);
double x1 = exp(-0.5*log(lastS()/sh) + y);
double x2 = exp(-0.5*log(lastS()/sh) - y);
Lorentz5Momentum P1 = lastParticles().first->momentum();
LorentzMomentum p1 = lightCone((P1.rho() + P1.e())*x1, Energy());
p1.rotateY(P1.theta());
p1.rotateZ(P1.phi());
meMomenta()[0] = p1;
Lorentz5Momentum P2 = lastParticles().second->momentum();
LorentzMomentum p2 = lightCone((P2.rho() + P2.e())*x2, Energy());
p2.rotateY(P2.theta());
p2.rotateZ(P2.phi());
meMomenta()[1] = p2;
lastXCombPtr()->lastX1X2(make_pair(x1,x2));
lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
return km*(ymax - ymin);
}
bool MatchboxMEBase::generateKinematics(const double * r) {
if ( phasespace() ) {
jacobian(phasespace()->generateKinematics(r,meMomenta()));
if ( jacobian() == 0.0 )
return false;
setScale();
logGenerateKinematics(r);
assert(lastMatchboxXComb());
if ( nDimAmplitude() > 0 ) {
amplitudeRandomNumbers().resize(nDimAmplitude());
copy(r + nDimPhasespace(),
r + nDimPhasespace() + nDimAmplitude(),
amplitudeRandomNumbers().begin());
}
if ( nDimInsertions() > 0 ) {
insertionRandomNumbers().resize(nDimInsertions());
copy(r + nDimPhasespace() + nDimAmplitude(),
r + nDimPhasespace() + nDimAmplitude() + nDimInsertions(),
insertionRandomNumbers().begin());
}
return true;
}
throw Exception()
<< "MatchboxMEBase::generateKinematics() expects a MatchboxPhasespace object.\n"
<< "Please check your setup." << Exception::abortnow;
return false;
}
int MatchboxMEBase::nDim() const {
if ( lastMatchboxXComb() )
return nDimPhasespace() + nDimAmplitude() + nDimInsertions();
int ampAdd = 0;
if ( matchboxAmplitude() ) {
ampAdd = matchboxAmplitude()->nDimAdditional();
}
int insertionAdd = 0;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
insertionAdd = max(insertionAdd,(**v).nDimAdditional());
}
return nDimBorn() + ampAdd + insertionAdd;
}
int MatchboxMEBase::nDimBorn() const {
if ( lastMatchboxXComb() )
return nDimPhasespace();
if ( phasespace() ) {
size_t nout = diagrams().front()->partons().size()-2;
int n = phasespace()->nDim(nout);
return n;
}
throw Exception()
<< "MatchboxMEBase::nDim() expects a MatchboxPhasespace object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0;
}
void MatchboxMEBase::setScale() const {
if ( haveX1X2() ) {
lastXCombPtr()->lastSHat((meMomenta()[0]+meMomenta()[1]).m2());
}
Energy2 fscale = factorizationScale()*sqr(factorizationScaleFactor());
Energy2 rscale = renormalizationScale()*sqr(renormalizationScaleFactor());
Energy2 ewrscale = renormalizationScaleQED();
lastXCombPtr()->lastScale(fscale);
if ( !fixedCouplings() ) {
if ( rscale > lastCuts().scaleMin() )
lastXCombPtr()->lastAlphaS(SM().alphaS(rscale));
else
lastXCombPtr()->lastAlphaS(SM().alphaS(lastCuts().scaleMin()));
} else {
lastXCombPtr()->lastAlphaS(SM().alphaS());
}
if ( !fixedQEDCouplings() ) {
lastXCombPtr()->lastAlphaEM(SM().alphaEM(ewrscale));
} else {
lastXCombPtr()->lastAlphaEM(SM().alphaEMMZ());
}
logSetScale();
}
Energy2 MatchboxMEBase::factorizationScale() const {
if ( scaleChoice() ) {
return scaleChoice()->factorizationScale();
}
throw Exception()
<< "MatchboxMEBase::factorizationScale() expects a MatchboxScaleChoice object.\n"
<< "Please check your setup." << Exception::abortnow;
return ZERO;
}
Energy2 MatchboxMEBase::renormalizationScale() const {
if ( scaleChoice() ) {
return scaleChoice()->renormalizationScale();
}
throw Exception()
<< "MatchboxMEBase::renormalizationScale() expects a MatchboxScaleChoice object.\n"
<< "Please check your setup." << Exception::abortnow;
return ZERO;
}
Energy2 MatchboxMEBase::renormalizationScaleQED() const {
if ( scaleChoice() ) {
return scaleChoice()->renormalizationScaleQED();
}
return renormalizationScale();
}
void MatchboxMEBase::setVetoScales(tSubProPtr) const {}
void MatchboxMEBase::getPDFWeight(Energy2 factorizationScale) const {
if ( !mePartonData()[0]->coloured() &&
!mePartonData()[1]->coloured() ) {
lastMEPDFWeight(1.0);
logPDFWeight();
return;
}
double w = 1.;
if ( mePartonData()[0]->coloured() && havePDFWeight1() )
w *= pdf1(factorizationScale);
if ( mePartonData()[1]->coloured() && havePDFWeight2() )
w *= pdf2(factorizationScale);
lastMEPDFWeight(w);
logPDFWeight();
}
double MatchboxMEBase::pdf1(Energy2 fscale, double xEx) const {
assert(lastXCombPtr()->partonBins().first->pdf());
if ( xEx < 1. && lastX1() >= xEx ) {
return
( ( 1. - lastX1() ) / ( 1. - xEx ) ) *
lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
lastPartons().first->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
xEx)/xEx;
}
return lastXCombPtr()->partonBins().first->pdf()->xfx(lastParticles().first->dataPtr(),
lastPartons().first->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
lastX1())/lastX1();
}
double MatchboxMEBase::pdf2(Energy2 fscale, double xEx) const {
assert(lastXCombPtr()->partonBins().second->pdf());
if ( xEx < 1. && lastX2() >= xEx ) {
return
( ( 1. - lastX2() ) / ( 1. - xEx ) ) *
lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
lastPartons().second->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
xEx)/xEx;
}
return lastXCombPtr()->partonBins().second->pdf()->xfx(lastParticles().second->dataPtr(),
lastPartons().second->dataPtr(),
fscale == ZERO ? lastScale() : fscale,
lastX2())/lastX2();
}
double MatchboxMEBase::me2() const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->me2()*
crossingSign()*
me2Norm());
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::me2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::finalStateSymmetry() const {
if ( symmetryFactor() > 0.0 )
return symmetryFactor();
double sFactor = 1.;
map<long,int> counts;
cPDVector checkData;
copy(mePartonData().begin()+2,mePartonData().end(),back_inserter(checkData));
cPDVector::iterator p = checkData.begin();
while ( !checkData.empty() ) {
if ( counts.find((**p).id()) != counts.end() ) {
counts[(**p).id()] += 1;
} else {
counts[(**p).id()] = 1;
}
checkData.erase(p);
p = checkData.begin();
continue;
}
for ( map<long,int>::const_iterator c = counts.begin();
c != counts.end(); ++c ) {
if ( c->second == 1 )
continue;
if ( c->second == 2 )
sFactor /= 2.;
else if ( c->second == 3 )
sFactor /= 6.;
else if ( c->second == 4 )
sFactor /= 24.;
}
symmetryFactor(sFactor);
return symmetryFactor();
}
double MatchboxMEBase::me2Norm(unsigned int addAlphaS) const {
// assume that we always have incoming
// spin-1/2 or massless spin-1 particles
double fac = 1./4.;
if ( orderInAlphaS() > 0 || addAlphaS != 0 )
fac *= pow(lastAlphaS()/SM().alphaS(),double(orderInAlphaS()+addAlphaS));
if ( orderInAlphaEW() > 0 )
fac *= pow(lastAlphaEM()/SM().alphaEM(),double(orderInAlphaEW()));
if ( mePartonData()[0]->iColour() == PDT::Colour3 ||
mePartonData()[0]->iColour() == PDT::Colour3bar )
fac /= SM().Nc();
else if ( mePartonData()[0]->iColour() == PDT::Colour8 )
fac /= (SM().Nc()*SM().Nc()-1.);
if ( mePartonData()[1]->iColour() == PDT::Colour3 ||
mePartonData()[1]->iColour() == PDT::Colour3bar )
fac /= SM().Nc();
else if ( mePartonData()[1]->iColour() == PDT::Colour8 )
fac /= (SM().Nc()*SM().Nc()-1.);
return finalStateSymmetry()*fac;
}
CrossSection MatchboxMEBase::dSigHatDR() const {
getPDFWeight();
if ( !lastXCombPtr()->willPassCuts() ) {
lastME2(0.0);
lastMECrossSection(ZERO);
return lastMECrossSection();
}
double xme2 = me2();
lastME2(xme2);
if ( xme2 == 0. && !oneLoopNoBorn() ) {
lastME2(0.0);
lastMECrossSection(ZERO);
return lastMECrossSection();
}
double vme2 = 0.;
if ( oneLoop() )
vme2 = oneLoopInterference();
CrossSection res = ZERO;
if ( !oneLoopNoBorn() )
res +=
(sqr(hbarc)/(2.*lastSHat())) *
jacobian()* lastMEPDFWeight() * xme2;
if ( oneLoop() )
res +=
(sqr(hbarc)/(2.*lastSHat())) *
jacobian()* lastMEPDFWeight() * vme2;
if ( !onlyOneLoop() ) {
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).setXComb(lastXCombPtr());
res += (**v).dSigHatDR();
}
}
double weight = 0.0;
bool applied = false;
for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).setXComb(lastXCombPtr());
if ( !(**rw).apply() )
continue;
weight += (**rw).evaluate();
applied = true;
}
if ( applied )
res *= weight;
lastMECrossSection(res);
return lastMECrossSection();
}
double MatchboxMEBase::oneLoopInterference() const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->oneLoopAmplitudes() )
matchboxAmplitude()->prepareOneLoopAmplitudes(this);
lastME2(matchboxAmplitude()->oneLoopInterference()*
crossingSign()*
me2Norm(1));
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::oneLoopInterference() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
void MatchboxMEBase::logPoles() const {
double res2me = oneLoopDoublePole();
double res1me = oneLoopSinglePole();
double res2i = 0.;
double res1i = 0.;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
res2i += (**v).oneLoopDoublePole();
res1i += (**v).oneLoopSinglePole();
}
double diff2 = abs(res2me) != 0. ? 1.-abs(res2i/res2me) : abs(res2i)-abs(res2me);
double diff1 = abs(res1me) != 0. ? 1.-abs(res1i/res1me) : abs(res1i)-abs(res1me);
generator()->log()
<< "check "
<< log10(abs(diff2)) << " " << log10(abs(diff1)) << "\n"
<< flush;
}
bool MatchboxMEBase::haveOneLoop() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->haveOneLoop();
return false;
}
bool MatchboxMEBase::onlyOneLoop() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->onlyOneLoop();
return false;
}
bool MatchboxMEBase::isDR() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->isDR();
return false;
}
bool MatchboxMEBase::isCS() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->isCS();
return false;
}
bool MatchboxMEBase::isBDK() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->isBDK();
return false;
}
bool MatchboxMEBase::isExpanded() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->isExpanded();
return false;
}
Energy2 MatchboxMEBase::mu2() const {
if ( matchboxAmplitude() )
return matchboxAmplitude()->mu2();
return 0*GeV2;
}
double MatchboxMEBase::oneLoopDoublePole() const {
if ( matchboxAmplitude() ) {
return
matchboxAmplitude()->oneLoopDoublePole()*
crossingSign()*
me2Norm(1);
}
return 0.;
}
double MatchboxMEBase::oneLoopSinglePole() const {
if ( matchboxAmplitude() ) {
return
matchboxAmplitude()->oneLoopSinglePole()*
crossingSign()*
me2Norm(1);
}
return 0.;
}
vector<Ptr<SubtractionDipole>::ptr>
MatchboxMEBase::getDipoles(const vector<Ptr<SubtractionDipole>::ptr>& dipoles,
const vector<Ptr<MatchboxMEBase>::ptr>& borns) const {
vector<Ptr<SubtractionDipole>::ptr> res;
// keep track of the dipoles we already did set up
set<pair<pair<pair<int,int>,int>,pair<Ptr<MatchboxMEBase>::tptr,Ptr<SubtractionDipole>::tptr> > > done;
cPDVector rep = diagrams().front()->partons();
int nreal = rep.size();
// now loop over configs
for ( int emitter = 0; emitter < nreal; ++emitter ) {
for ( int spectator = 0; spectator < nreal; ++spectator ) {
if ( emitter == spectator )
continue;
for ( int emission = 2; emission < nreal; ++emission ) {
if ( emission == emitter || emission == spectator )
continue;
for ( vector<Ptr<MatchboxMEBase>::ptr>::const_iterator b =
borns.begin(); b != borns.end(); ++b ) {
if ( (**b).onlyOneLoop() )
continue;
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
dipoles.begin(); d != dipoles.end(); ++d ) {
if ( !rep[emitter]->coloured() ||
!rep[emission]->coloured() ||
!rep[spectator]->coloured() ) {
continue;
}
if ( noDipole(emitter,emission,spectator) ) {
continue;
}
if ( done.find(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(*b,*d)))
!= done.end() ) {
continue;
}
if ( !(**d).canHandle(rep,emitter,emission,spectator) ) {
continue;
}
// now get to work
Ptr<SubtractionDipole>::ptr nDipole = (**d).cloneMe();
nDipole->realEmitter(emitter);
nDipole->realEmission(emission);
nDipole->realSpectator(spectator);
nDipole->realEmissionME(const_cast<MatchboxMEBase*>(this));
nDipole->underlyingBornME(*b);
nDipole->setupBookkeeping();
if ( !(nDipole->empty()) ) {
res.push_back(nDipole);
done.insert(make_pair(make_pair(make_pair(emitter,emission),spectator),make_pair(*b,*d)));
if ( nDipole->isSymmetric() )
done.insert(make_pair(make_pair(make_pair(emission,emitter),spectator),make_pair(*b,*d)));
ostringstream dname;
dname << fullName() << "." << (**b).name() << "."
<< (**d).name() << ".[("
<< emitter << "," << emission << ")," << spectator << "]";
if ( ! (generator()->preinitRegister(nDipole,dname.str()) ) )
throw InitException() << "Dipole " << dname.str() << " already existing.";
nDipole->cloneDependencies(dname.str());
}
}
}
}
}
}
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d = res.begin();
d != res.end(); ++d )
(**d).partnerDipoles(res);
return res;
}
double MatchboxMEBase::colourCorrelatedME2(pair<int,int> ij) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->colourCorrelatedME2(ij)*
crossingSign()*
me2Norm());
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::colourCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::largeNColourCorrelatedME2(pair<int,int> ij,
Ptr<ColourBasis>::tptr largeNBasis) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->largeNColourCorrelatedME2(ij,largeNBasis)*
crossingSign()*
me2Norm());
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::largeNColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
double MatchboxMEBase::spinColourCorrelatedME2(pair<int,int> ij,
const SpinCorrelationTensor& c) const {
if ( matchboxAmplitude() ) {
if ( matchboxAmplitude()->treeAmplitudes() )
matchboxAmplitude()->prepareAmplitudes(this);
lastME2(matchboxAmplitude()->spinColourCorrelatedME2(ij,c)*
crossingSign()*
me2Norm());
logME2();
return lastME2();
}
throw Exception()
<< "MatchboxMEBase::spinColourCorrelatedME2() expects a MatchboxAmplitude object.\n"
<< "Please check your setup." << Exception::abortnow;
return 0.;
}
void MatchboxMEBase::flushCaches() {
MEBase::flushCaches();
if ( matchboxAmplitude() )
matchboxAmplitude()->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
reweights().begin(); r != reweights().end(); ++r ) {
(**r).flushCaches();
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
(**v).flushCaches();
}
}
void MatchboxMEBase::print(ostream& os) const {
os << "--- MatchboxMEBase setup -------------------------------------------------------\n";
os << " '" << name() << "' for subprocesses:\n";
for ( vector<PDVector>::const_iterator p = subProcesses().begin();
p != subProcesses().end(); ++p ) {
os << " ";
for ( PDVector::const_iterator pp = p->begin();
pp != p->end(); ++pp ) {
os << (**pp).PDGName() << " ";
if ( pp == p->begin() + 1 )
os << "-> ";
}
os << "\n";
}
os << " including " << (oneLoop() ? "" : "no ") << "virtual corrections";
if ( oneLoopNoBorn() )
os << " without Born contributions";
os << "\n";
if ( oneLoop() && !onlyOneLoop() ) {
os << " using insertion operators\n";
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
os << " '" << (**v).name() << "' with "
<< ((**v).isDR() ? "" : "C") << "DR/";
if ( (**v).isCS() )
os << "CS";
if ( (**v).isBDK() )
os << "BDK";
if ( (**v).isExpanded() )
os << "expanded";
os << " conventions\n";
}
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxMEBase::printLastEvent(ostream& os) const {
os << "--- MatchboxMEBase last event information --------------------------------------\n";
os << " for matrix element '" << name() << "'\n";
os << " process considered:\n ";
int in = 0;
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p ) {
os << (**p).PDGName() << " ";
if ( ++in == 2 )
os << " -> ";
}
os << " kinematic environment as set by the XComb " << lastXCombPtr() << ":\n"
<< " sqrt(shat)/GeV = " << sqrt(lastSHat()/GeV2)
<< " x1 = " << lastX1() << " x2 = " << lastX2()
<< " alphaS = " << lastAlphaS() << "\n";
os << " momenta/GeV generated from random numbers\n ";
copy(lastXComb().lastRandomNumbers().begin(),
lastXComb().lastRandomNumbers().end(),ostream_iterator<double>(os," "));
os << ":\n ";
for ( vector<Lorentz5Momentum>::const_iterator p = meMomenta().begin();
p != meMomenta().end(); ++p ) {
os << (*p/GeV) << "\n ";
}
os << "last cross section/nb calculated was:\n "
<< (lastMECrossSection()/nanobarn) << " (pdf weight " << lastMEPDFWeight() << ")\n";
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void MatchboxMEBase::logGenerateKinematics(const double * r) const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' generated kinematics\nfrom "
<< nDim() << " random numbers:\n";
copy(r,r+nDim(),ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "storing phase space information in XComb "
<< lastXCombPtr() << "\n";
generator()->log() << "generated phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
cPDVector::const_iterator dit = mePartonData().begin();
for ( ; pit != meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
<< "and Jacobian = " << jacobian() << " sHat/GeV2 = "
<< (lastSHat()/GeV2) << "\n" << flush;
}
void MatchboxMEBase::logSetScale() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' set scales using XComb " << lastXCombPtr() << ":\n"
<< "scale/GeV2 = " << (scale()/GeV2) << " xi_R = "
<< renormalizationScaleFactor() << " xi_F = "
<< factorizationScaleFactor() << "\n"
<< "alpha_s = " << lastAlphaS() << "\n" << flush;
}
void MatchboxMEBase::logPDFWeight() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' calculated pdf weight = "
<< lastMEPDFWeight() << " from XComb "
<< lastXCombPtr() << "\n"
<< "x1 = " << lastX1() << " (" << (mePartonData()[0]->coloured() ? "" : "not ") << "used) "
<< "x2 = " << lastX2() << " (" << (mePartonData()[1]->coloured() ? "" : "not ") << "used)\n"
<< flush;
}
void MatchboxMEBase::logME2() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' evaluated me2 using XComb "
<< lastXCombPtr() << "\n"
<< "and phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = meMomenta().begin();
cPDVector::const_iterator dit = mePartonData().begin();
for ( ; pit != meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << lastX1() << " x2 = " << lastX2() << "\n"
<< "sHat/GeV2 = " << (lastSHat()/GeV2)
<< " me2 = " << lastME2() << "\n" << flush;
}
void MatchboxMEBase::logDSigHatDR() const {
if ( !verbose() )
return;
generator()->log() << "'" << name() << "' evaluated cross section using XComb "
<< lastXCombPtr() << "\n"
<< "Jacobian = " << jacobian() << " sHat/GeV2 = "
<< (lastSHat()/GeV2) << " dsig/nb = "
<< (lastMECrossSection()/nanobarn) << "\n" << flush;
}
void MatchboxMEBase::cloneDependencies(const std::string& prefix) {
if ( phasespace() ) {
Ptr<MatchboxPhasespace>::ptr myPhasespace = phasespace()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myPhasespace->name();
if ( ! (generator()->preinitRegister(myPhasespace,pname.str()) ) )
throw InitException() << "Phasespace generator " << pname.str() << " already existing.";
myPhasespace->cloneDependencies(pname.str());
phasespace(myPhasespace);
}
theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
if ( matchboxAmplitude() ) {
Ptr<MatchboxAmplitude>::ptr myAmplitude = matchboxAmplitude()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myAmplitude->name();
if ( ! (generator()->preinitRegister(myAmplitude,pname.str()) ) )
throw InitException() << "Amplitude " << pname.str() << " already existing.";
myAmplitude->cloneDependencies(pname.str());
matchboxAmplitude(myAmplitude);
amplitude(myAmplitude);
matchboxAmplitude()->orderInGs(orderInAlphaS());
matchboxAmplitude()->orderInGem(orderInAlphaEW());
}
if ( scaleChoice() ) {
Ptr<MatchboxScaleChoice>::ptr myScaleChoice = scaleChoice()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myScaleChoice->name();
if ( ! (generator()->preinitRegister(myScaleChoice,pname.str()) ) )
throw InitException() << "Scale choice " << pname.str() << " already existing.";
scaleChoice(myScaleChoice);
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
Ptr<MatchboxReweightBase>::ptr myReweight = (**rw).cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << (**rw).name();
if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
throw InitException() << "Reweight " << pname.str() << " already existing.";
myReweight->cloneDependencies(pname.str());
*rw = myReweight;
}
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
Ptr<MatchboxInsertionOperator>::ptr myIOP = (**v).cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << (**v).name();
if ( ! (generator()->preinitRegister(myIOP,pname.str()) ) )
throw InitException() << "Insertion operator " << pname.str() << " already existing.";
*v = myIOP;
}
}
void MatchboxMEBase::prepareXComb(MatchboxXCombData& xc) const {
if ( phasespace() ) {
size_t nout = diagrams().front()->partons().size()-2;
xc.nDimPhasespace(phasespace()->nDim(nout));
}
if ( matchboxAmplitude() ) {
xc.nDimAmplitude(matchboxAmplitude()->nDimAdditional());
if ( matchboxAmplitude()->colourBasis() ) {
size_t cdim =
matchboxAmplitude()->colourBasis()->prepare(diagrams(),matchboxAmplitude()->noCorrelations());
xc.colourBasisDim(cdim);
}
}
int insertionAdd = 0;
for ( vector<Ptr<MatchboxInsertionOperator>::ptr>::const_iterator v =
virtuals().begin(); v != virtuals().end(); ++v ) {
insertionAdd = max(insertionAdd,(**v).nDimAdditional());
}
xc.nDimInsertions(insertionAdd);
xc.nLight(getNLight());
if ( verbose() ) {
string fname = name() + ".diagrams";
ifstream test(fname.c_str());
if ( !test ) {
test.close();
ofstream out(fname.c_str());
for ( vector<Ptr<DiagramBase>::ptr>::const_iterator d = diagrams().begin();
d != diagrams().end(); ++d ) {
DiagramDrawer::drawDiag(out,dynamic_cast<const Tree2toNDiagram&>(**d));
out << "\n";
}
}
}
}
StdXCombPtr MatchboxMEBase::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec&,
tStdXCombPtr newHead,
tMEPtr newME) {
if ( !newME )
newME = this;
Ptr<MatchboxXComb>::ptr xc =
new_ptr(MatchboxXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts, newME,
newDiagrams, mir,
newHead));
prepareXComb(*xc);
return xc;
}
StdXCombPtr MatchboxMEBase::makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME) {
if ( !newME )
newME = this;
Ptr<MatchboxXComb>::ptr xc =
new_ptr(MatchboxXComb(newHead, newPartonBins, newME, newDiagrams));
prepareXComb(*xc);
return xc;
}
void MatchboxMEBase::persistentOutput(PersistentOStream & os) const {
- os << theLastXComb
- << theFactory << thePhasespace << theAmplitude << theScaleChoice
- << theVirtuals << theReweights << theSubprocesses
- << theOneLoop << theOneLoopNoBorn;
+ os << theLastXComb << theFactory << thePhasespace
+ << theAmplitude << theScaleChoice << theVirtuals
+ << theReweights << theSubprocesses << theOneLoop
+ << theOneLoopNoBorn;
}
void MatchboxMEBase::persistentInput(PersistentIStream & is, int) {
- is >> theLastXComb
- >> theFactory >> thePhasespace >> theAmplitude >> theScaleChoice
- >> theVirtuals >> theReweights >> theSubprocesses
- >> theOneLoop >> theOneLoopNoBorn;
+ is >> theLastXComb >> theFactory >> thePhasespace
+ >> theAmplitude >> theScaleChoice >> theVirtuals
+ >> theReweights >> theSubprocesses >> theOneLoop
+ >> theOneLoopNoBorn;
lastMatchboxXComb(theLastXComb);
}
void MatchboxMEBase::Init() {
static ClassDocumentation<MatchboxMEBase> documentation
("MatchboxMEBase is the base class for matrix elements "
"in the context of the matchbox NLO interface.");
}
IBPtr MatchboxMEBase::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxMEBase::fullclone() const {
return new_ptr(*this);
}
void MatchboxMEBase::doinit() {
MEBase::doinit();
if ( !theAmplitude )
theAmplitude = dynamic_ptr_cast<Ptr<MatchboxAmplitude>::ptr>(amplitude());
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<MatchboxMEBase,MEBase>
describeHerwigMatchboxMEBase("Herwig::MatchboxMEBase", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Base/SubtractedME.cc b/MatrixElement/Matchbox/Base/SubtractedME.cc
--- a/MatrixElement/Matchbox/Base/SubtractedME.cc
+++ b/MatrixElement/Matchbox/Base/SubtractedME.cc
@@ -1,707 +1,717 @@
// -*- C++ -*-
//
// SubtractedME.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SubtractedME class.
//
#include "SubtractedME.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Handlers/StdXCombGroup.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/MatrixElement/Matchbox/Base/DipoleRepository.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/MatchboxXCombGroup.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.h"
using namespace Herwig;
SubtractedME::SubtractedME()
: MEGroup(),
theRealShowerSubtraction(false), theVirtualShowerSubtraction(false) {}
SubtractedME::~SubtractedME() {}
Ptr<MatchboxFactory>::tcptr SubtractedME::factory() const { return theFactory; }
void SubtractedME::factory(Ptr<MatchboxFactory>::tcptr f) { theFactory = f; }
bool SubtractedME::subProcessGroups() const {
return
(factory()->subProcessGroups() && !(showerApproximation() || inclusive())) ||
factory()->subtractionData() != "";
}
bool SubtractedME::inclusive() const { return factory()->inclusive(); }
Ptr<ShowerApproximation>::tptr SubtractedME::showerApproximation() const { return factory()->showerApproximation(); }
const vector<Ptr<MatchboxMEBase>::ptr>& SubtractedME::borns() const { return factory()->bornMEs(); }
bool SubtractedME::verbose() const { return factory()->verbose(); }
IBPtr SubtractedME::clone() const {
return new_ptr(*this);
}
IBPtr SubtractedME::fullclone() const {
return new_ptr(*this);
}
StdXCombPtr SubtractedME::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allPBins,
tStdXCombPtr newHead,
tMEPtr newME) {
tMEGroupPtr newMEGroup = dynamic_ptr_cast<tMEGroupPtr>(newME);
if ( !newMEGroup )
newMEGroup = this;
Ptr<MatchboxXCombGroup>::ptr res =
new_ptr(MatchboxXCombGroup(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts, newMEGroup,
newDiagrams, mir,
newHead));
res->build(allPBins);
theReal->prepareXComb(*res);
if ( factory()->subtractionData() != "" ) {
set<cPDVector> procs;
for ( DiagramVector::const_iterator d = head()->diagrams().begin();
d != head()->diagrams().end(); ++d ) {
if ( procs.find((**d).partons()) == procs.end() )
procs.insert((**d).partons());
}
for ( set<cPDVector>::const_iterator p = procs.begin();
p != procs.end(); ++p ) {
for ( size_t i = 0; i < (*p).size(); ++i ) {
if ( !(*p)[i]->coloured() )
continue;
if ( i > 1 &&
(*p)[i]->id() == ParticleID::g ) {
softHistograms[SoftSubtractionIndex(*p,i)] = SubtractionHistogram(0.001,10.);
if ( theReal->phasespace() )
- singularLimits().insert(make_pair(i,i));
+ res->singularLimits().insert(make_pair(i,i));
}
for ( size_t j = i+1; j < (*p).size(); ++j ) {
if ( !(*p)[j]->coloured() )
continue;
long iid = (*p)[i]->id();
long jid = (*p)[j]->id();
if ( i < 2 && j < 2 )
continue;
if ( i < 2 && j > 1 ) {
if ( abs(iid) < 7 && abs(jid) < 7 && iid != jid )
continue;
}
if ( i > 1 && j > 1 ) {
if ( abs(iid) < 7 && abs(jid) < 7 && iid + jid != 0 )
continue;
}
bool haveDipole = false;
for ( MEVector::const_iterator k = dependent().begin();
k != dependent().end(); ++k ) {
const SubtractionDipole& dip = dynamic_cast<const SubtractionDipole&>(**k);
if ( ( (size_t)(dip.realEmitter()) == i && (size_t)(dip.realEmission()) == j ) ||
( (size_t)(dip.realEmitter()) == j && (size_t)(dip.realEmission()) == i ) ) {
haveDipole = true;
break;
}
}
if ( !haveDipole )
continue;
collinearHistograms[CollinearSubtractionIndex(*p,make_pair(i,j))] = SubtractionHistogram(0.001,20.);
if ( theReal->phasespace() )
- singularLimits().insert(make_pair(i,j));
+ res->singularLimits().insert(make_pair(i,j));
}
}
}
}
return res;
}
void SubtractedME::setXComb(tStdXCombPtr xc) {
MEGroup::setXComb(xc);
lastMatchboxXComb(xc);
}
MEBase::DiagramVector SubtractedME::dependentDiagrams(const cPDVector& proc,
tMEPtr depME) const {
Ptr<SubtractionDipole>::tptr dipole =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(depME);
if ( !dipole ) {
Throw<InitException>() << "A dependent matrix element of SubtractedME "
<< "has not been derived from SubtractionDipole. "
<< "Please check the corresponding input file.";
}
return dipole->underlyingBornDiagrams(proc);
}
vector<Ptr<SubtractionDipole>::ptr> SubtractedME::dipoles() {
if ( dependent().empty() )
getDipoles();
vector<Ptr<SubtractionDipole>::ptr> res;
for ( MEVector::const_iterator k = dependent().begin();
k != dependent().end(); ++k )
res.push_back(dynamic_ptr_cast<Ptr<SubtractionDipole>::ptr>(*k));
return res;
}
void SubtractedME::getDipoles() {
if ( !dependent().empty() )
return;
Ptr<MatchboxMEBase>::tptr real =
dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( borns().empty() || !real )
Throw<InitException>() << "The SubtractedME '"
<< name() << "' could not generate "
<< "subtraction terms for the real emission "
<< "matrix element '" << real->name() << "'. "
<< "Please check the corresponding input file.";
Ptr<MatchboxMEBase>::ptr myRealEmissionME = real->cloneMe();
ostringstream pname;
pname << fullName() << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
myRealEmissionME->cloneDependencies(pname.str());
head(myRealEmissionME);
real = myRealEmissionME;
MEVector dipMEs;
vector<Ptr<SubtractionDipole>::ptr> genDipoles
= real->getDipoles(DipoleRepository::dipoles(),borns());
if ( factory()->subtractionData() != "" ) {
for ( vector<Ptr<SubtractionDipole>::ptr>::const_iterator d =
genDipoles.begin(); d != genDipoles.end(); ++d )
(**d).doTestSubtraction();
}
if ( genDipoles.empty() ) {
// probably finite real contribution, but warn
generator()->log() << "\nWarning: No subtraction dipoles could be found for the processes:\n";
for ( vector<PDVector>::const_iterator s = real->subProcesses().begin();
s != real->subProcesses().end(); ++s ) {
generator()->log() << (*s)[0]->PDGName() << " " << (*s)[1]->PDGName()
<< " -> ";
for ( PDVector::const_iterator p = s->begin() + 2; p != s->end(); ++p )
generator()->log() << (**p).PDGName() << " ";
generator()->log() << "\n" << flush;
}
generator()->log() << "Assuming finite tree-level O(alphaS) correction.\n";
}
dipMEs.resize(genDipoles.size());
copy(genDipoles.begin(),genDipoles.end(),dipMEs.begin());
dependent() = dipMEs;
}
void SubtractedME::cloneRealME(const string& prefix) {
theReal = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( theReal ) {
Ptr<MatchboxMEBase>::ptr myRealEmissionME = theReal->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
myRealEmissionME->cloneDependencies(pname.str());
theReal = myRealEmissionME;
}
head(theReal);
}
void SubtractedME::cloneDipoles(const string& prefix) {
MEVector dipMEs;
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
Ptr<SubtractionDipole>::ptr cloned = dip->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << cloned->name();
if ( ! (generator()->preinitRegister(cloned,pname.str()) ) )
throw InitException() << "Subtraction dipole " << pname.str() << " already existing.";
cloned->cloneDependencies(pname.str());
dipMEs.push_back(cloned);
}
dependent() = dipMEs;
}
vector<Ptr<SubtractionDipole>::ptr> SubtractedME::splitDipoles(const cPDVector& born) {
vector<Ptr<SubtractionDipole>::ptr> dips = dipoles();
vector<Ptr<SubtractionDipole>::ptr> res;
for ( vector<Ptr<SubtractionDipole>::ptr>::iterator d = dips.begin();
d != dips.end(); ++d ) {
for ( DiagramVector::const_iterator p = (**d).underlyingBornME()->diagrams().begin();
p != (**d).underlyingBornME()->diagrams().end(); ++p )
if ( born == (**p).partons() ) {
res.push_back(*d);
break;
}
}
return res;
}
void SubtractedME::doRealEmissionScales() {
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
dip->doRealEmissionScales();
}
}
void SubtractedME::doRealShowerSubtraction() {
theRealShowerSubtraction = true;
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
dip->showerApproximation(showerApproximation());
dip->doRealShowerSubtraction();
}
}
void SubtractedME::doVirtualShowerSubtraction() {
theVirtualShowerSubtraction = true;
for ( MEVector::const_iterator m = dependent().begin();
m != dependent().end(); ++m ) {
Ptr<SubtractionDipole>::tptr dip =
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*m);
assert(dip);
dip->showerApproximation(showerApproximation());
dip->doVirtualShowerSubtraction();
}
}
void SubtractedME::setVetoScales(tSubProPtr) const {}
void SubtractedME::fillProjectors() {
if ( !inclusive() && !virtualShowerSubtraction() )
return;
Ptr<StdXCombGroup>::tptr group =
dynamic_ptr_cast<Ptr<StdXCombGroup>::tptr>(lastXCombPtr());
for ( vector<StdXCombPtr>::const_iterator d = group->dependent().begin();
d != group->dependent().end(); ++d ) {
if ( !(**d).kinematicsGenerated() )
continue;
if ( (**d).willPassCuts() )
lastXCombPtr()->projectors().insert(1.,*d);
}
}
double SubtractedME::reweightHead(const vector<tStdXCombPtr>& dep) {
if ( inclusive() && !lastXComb().lastProjector() )
return 1.;
if ( virtualShowerSubtraction() && !lastXComb().lastProjector() ) {
return 0.;
}
if ( realShowerSubtraction() ) {
assert(showerApproximation());
bool below = showerApproximation()->belowCutoff();
bool haveDipole = false;
for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d ) {
if ( !(**d).kinematicsGenerated() )
continue;
if ( (**d).willPassCuts() ) {
haveDipole = true;
break;
}
}
if ( !haveDipole )
return 1.;
if ( below )
return 0.;
return 1.;
}
if ( virtualShowerSubtraction() || inclusive() ) {
if ( virtualShowerSubtraction() ) {
assert(showerApproximation());
bool above = !showerApproximation()->belowCutoff();
if ( above )
return 0.;
}
double sum = 0.;
size_t n = 0;
for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d ) {
if ( !(**d).kinematicsGenerated() )
continue;
if ( (**d).willPassCuts() ) {
sum += (**d).lastME2();
++n;
}
}
return
n * lastXComb().lastProjector()->lastME2() / sum;
}
return 1.;
}
double SubtractedME::reweightDependent(tStdXCombPtr xc, const vector<tStdXCombPtr>& dep) {
if ( inclusive() && !lastXComb().lastProjector() )
return 0.;
if ( virtualShowerSubtraction() && !lastXComb().lastProjector() ) {
return 0.;
}
if ( virtualShowerSubtraction() || inclusive() ) {
if ( xc != lastXComb().lastProjector() )
return 0.;
size_t n = 0;
for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d ) {
if ( !(**d).kinematicsGenerated() )
continue;
if ( (**d).willPassCuts() ) {
++n;
}
}
return n;
}
if ( realShowerSubtraction() ) {
bool haveDipole = false;
for ( vector<tStdXCombPtr>::const_iterator d = dep.begin(); d != dep.end(); ++d ) {
if ( !(**d).kinematicsGenerated() )
continue;
if ( (**d).willPassCuts() ) {
haveDipole = true;
break;
}
}
if ( !haveDipole )
return 0.;
}
return 1.;
}
void SubtractedME::doinit() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::doinit();
return;
}
theReal = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
if ( theReal ) {
getDipoles();
}
if ( verbose() )
print(Repository::clog());
MEGroup::doinit();
}
void SubtractedME::doinitrun() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::doinitrun();
return;
}
MEGroup::doinitrun();
theReal = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head());
}
void SubtractedME::dofinish() {
// has been deactivated by the factory
if ( !head() ) {
MEBase::dofinish();
return;
}
MEGroup::dofinish();
for ( map<CollinearSubtractionIndex,SubtractionHistogram>::
const_iterator b = collinearHistograms.begin();
b != collinearHistograms.end(); ++b ) {
b->second.dump(factory()->subtractionData(),
b->first.first,
b->first.second.first,
b->first.second.second);
}
for ( map<SoftSubtractionIndex,SubtractionHistogram>::
const_iterator b = softHistograms.begin();
b != softHistograms.end(); ++b ) {
b->second.dump(factory()->subtractionData(),
b->first.first,
b->first.second,
b->first.second);
}
}
void SubtractedME::print(ostream& os) const {
os << "--- SubtractedME setup ---------------------------------------------------------\n";
os << " '" << name() << "' subtracting real emission\n '"
<< head()->name() << "' using the dipoles:\n";
for ( MEVector::const_iterator d = dependent().begin();
d != dependent().end(); ++d )
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*d)->print(os);
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void SubtractedME::printLastEvent(ostream& os) const {
os << "--- SubtractedME last event information ----------------------------------------\n";
os << " for subtracted matrix element '" << name() << "'\n";
os << " real emission event information:\n";
dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(head())->printLastEvent(os);
os << " dipoles event information:\n";
for ( MEVector::const_iterator d = dependent().begin();
d != dependent().end(); ++d )
dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(*d)->printLastEvent(os);
os << "--- end SubtractedME last event information ------------------------------------\n\n\n";
os << flush;
}
void SubtractedME::lastEventStatistics() {
MEGroup::lastEventStatistics();
if ( !generator() )
return;
/*
if ( verbose() )
printLastEvent(generator()->log());
*/
if ( !collinearHistograms.empty() )
lastEventSubtraction();
}
SubtractedME::SubtractionHistogram::
SubtractionHistogram(double low,
double up,
unsigned int nbins)
: lower(low) {
nbins = nbins + 1;
double c = log10(up/low) / (nbins-1.);
for ( unsigned int k = 1; k < nbins; ++k ) {
bins[low*pow(10.0,k*c)] = make_pair(Constants::MaxDouble,0.);
}
}
+void SubtractedME::SubtractionHistogram::persistentOutput(PersistentOStream& os) const {
+ os << lower << bins;
+}
+
+void SubtractedME::SubtractionHistogram::persistentInput(PersistentIStream& is) {
+ is >> lower >> bins;
+}
+
void SubtractedME::SubtractionHistogram::
dump(const std::string& prefix,
const cPDVector& proc,
int i, int j) const {
ostringstream fname("");
for ( cPDVector::const_iterator p = proc.begin();
p != proc.end(); ++p )
fname << (**p).PDGName();
fname << "-" << i << "-" << j;
ofstream out((prefix+fname.str()+".dat").c_str());
for ( map<double,pair<double,double> >::const_iterator b = bins.begin();
b != bins.end(); ++b ) {
map<double,pair<double,double> >::const_iterator bp = b; --bp;
if ( b->second.first != Constants::MaxDouble &&
b->second.second != 0. ) {
if ( b != bins.begin() )
out << bp->first;
else
out << lower;
out << " " << b->first
<< " " << b->second.first
<< " " << b->second.second
<< "\n" << flush;
}
}
ofstream gpout((prefix+fname.str()+".gp").c_str());
gpout << "set terminal epslatex color solid;\n"
<< "set output '" << fname.str() << "-plot.tex';\n"
<< "set log x;\n"
<< "set size 0.5,0.6;\n"
<< "set yrange [0:2];\n"
<< "set xrange [0.001:10];\n";
if ( i != j ) {
gpout << "set xlabel '$\\sqrt{s_{" << i << j << "}}/{\\rm GeV}$'\n";
} else {
gpout << "set xlabel '$E_{" << i << "}/{\\rm GeV}$'\n";
}
gpout << "plot 1 w lines lc rgbcolor \"#DDDDDD\" notitle, '" << fname.str()
<< ".dat' u (($1+$2)/2.):3:($4 < 4. ? $4 : 4.) w filledcurves lc rgbcolor \"#00AACC\" t "
<< "'$";
for ( size_t k = 0; k < proc.size(); k++ ) {
if ( k == 2 )
gpout << "\\to ";
gpout << (proc[k]->id() < 0 ? "\\bar{" : "")
<< (proc[k]->id() < 0 ? proc[k]->CC()->PDGName() : proc[k]->PDGName())
<< (proc[k]->id() < 0 ? "}" : "") << " ";
}
gpout << "$';\n";
gpout << "reset;\n";
}
void SubtractedME::lastEventSubtraction() {
tStdXCombGroupPtr xc = dynamic_ptr_cast<tStdXCombGroupPtr>(lastXCombPtr());
CrossSection xcme2 = xc->lastHeadCrossSection();
CrossSection xcdip = ZERO;
if ( xcme2 == ZERO )
return;
for ( vector<StdXCombPtr>::const_iterator d = xc->dependent().begin();
d != xc->dependent().end(); ++d ) {
if ( !(*d) )
continue;
if ( !(**d).matrixElement()->apply() )
continue;
+ if ( !(**d).willPassCuts() )
+ continue;
xcdip += (**d).lastCrossSection();
}
if ( theReal->phasespace() ) {
size_t i = lastSingularLimit()->first;
size_t j = lastSingularLimit()->second;
if ( i == j &&
softHistograms.find(SoftSubtractionIndex(head()->mePartonData(),i))
!= softHistograms.end() ) {
softHistograms[SoftSubtractionIndex(head()->mePartonData(),i)].
book(meMomenta()[i].t()/GeV,abs(xcdip)/abs(xcme2));
}
if ( i != j &&
collinearHistograms.find(CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j)))
!= collinearHistograms.end() ) {
double s = sqrt(2.*meMomenta()[i]*meMomenta()[j])/GeV;
collinearHistograms[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))].
book(s,abs(xcdip)/abs(xcme2));
}
return;
}
for ( size_t i = 0; i < meMomenta().size(); ++i ) {
if ( i > 1 ) {
if ( softHistograms.find(SoftSubtractionIndex(head()->mePartonData(),i))
!= softHistograms.end() ) {
softHistograms[SoftSubtractionIndex(head()->mePartonData(),i)].
book(meMomenta()[i].t()/GeV,abs(xcdip)/abs(xcme2));
}
}
for ( size_t j = i+1; j < meMomenta().size(); ++j ) {
if ( collinearHistograms.find(CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j)))
== collinearHistograms.end() )
continue;
double s = sqrt(2.*meMomenta()[i]*meMomenta()[j])/GeV;
collinearHistograms[CollinearSubtractionIndex(head()->mePartonData(),make_pair(i,j))].
book(s,abs(xcdip)/abs(xcme2));
}
}
}
void SubtractedME::persistentOutput(PersistentOStream & os) const {
- os << theLastXComb
- << theFactory << theReal
+ os << theLastXComb << theFactory << theReal
+ << collinearHistograms << softHistograms
<< theRealShowerSubtraction << theVirtualShowerSubtraction;
}
void SubtractedME::persistentInput(PersistentIStream & is, int) {
- is >> theLastXComb
- >> theFactory >> theReal
+ is >> theLastXComb >> theFactory >> theReal
+ >> collinearHistograms >> softHistograms
>> theRealShowerSubtraction >> theVirtualShowerSubtraction;
lastMatchboxXComb(theLastXComb);
}
void SubtractedME::Init() {
static ClassDocumentation<SubtractedME> documentation
("SubtractedME represents a subtracted real emission matrix element.");
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<SubtractedME,MEGroup>
describeHerwigSubtractedME("Herwig::SubtractedME", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Base/SubtractedME.h b/MatrixElement/Matchbox/Base/SubtractedME.h
--- a/MatrixElement/Matchbox/Base/SubtractedME.h
+++ b/MatrixElement/Matchbox/Base/SubtractedME.h
@@ -1,469 +1,491 @@
// -*- C++ -*-
//
// SubtractedME.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SubtractedME_H
#define HERWIG_SubtractedME_H
//
// This is the declaration of the SubtractedME class.
//
#include "ThePEG/MatrixElement/MEGroup.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
#include "Herwig++/MatrixElement/Matchbox/Utility/LastMatchboxXCombInfo.h"
#include "Herwig++/MatrixElement/Matchbox/MatchboxFactory.fh"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief SubtractedME represents a subtracted real emission matrix element.
*
* @see \ref SubtractedMEInterfaces "The interfaces"
* defined for SubtractedME.
*/
class SubtractedME:
public MEGroup,
public LastMatchboxXCombInfo {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
SubtractedME();
/**
* The destructor.
*/
virtual ~SubtractedME();
//@}
public:
/**
* Return the factory which produced this matrix element
*/
Ptr<MatchboxFactory>::tcptr factory() const;
/**
* Set the factory which produced this matrix element
*/
void factory(Ptr<MatchboxFactory>::tcptr f);
/** @name Phasespace and subprocess information */
//@{
/**
* For the given event generation setup return an xcomb object
* appropriate to this matrix element.
*/
virtual StdXCombPtr makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allPBins,
tStdXCombPtr newHead = tStdXCombPtr(),
tMEPtr newME = tMEPtr());
/**
* Set the XComb object to be used in the next call to
* generateKinematics() and dSigHatDR().
*/
virtual void setXComb(tStdXCombPtr);
/**
* Return true, if the same additional random numbers
* should be presented to any of the dependent
* matrix elements.
*/
virtual bool uniformAdditional() const { return true; }
/**
* Return true, if the XComb steering this matrix element
* should keep track of the random numbers used to generate
* the last phase space point
*/
virtual bool keepRandomNumbers() const { return true; }
/**
* Given a process from the head matrix element,
* return a list of diagrams which should be considered for
* the given dependent matrix element.
*/
virtual MEBase::DiagramVector dependentDiagrams(const cPDVector& proc,
tMEPtr depME) const;
/**
* Return true, if SubProcessGroups should be
* setup from this MEGroup. If not, a single SubProcess
* is constructed from the data provided by the
* head matrix element.
*/
virtual bool subProcessGroups() const;
/**
* Return true, if one of the dependent subprocesses should be
* constructed in place of the one driven by the head matrix element
* or a full subprocess group.
*/
virtual bool selectDependentSubProcess() const { return inclusive(); }
/**
* Return true, if the integral over the unresolved emission should be
* calculated.
*/
bool inclusive() const;
/**
* Fill the projectors object of xcombs to choose subprocesses
* different than the one currently integrated.
*/
virtual void fillProjectors();
/**
* Return true, if this MEGroup will reweight the contributing cross
* sections.
*/
virtual bool groupReweighted() const { return inclusive() || showerApproximation(); }
/**
* Reweight the head cross section
*/
virtual double reweightHead(const vector<tStdXCombPtr>&);
/**
* Reweight the dependent cross section
*/
virtual double reweightDependent(tStdXCombPtr, const vector<tStdXCombPtr>&);
/**
* Switch on or off that scales should be calculated from real emission kinematics
*/
void doRealEmissionScales();
//@}
/** @name Methods relevant to matching */
//@{
/**
* Inform this matrix element that a new phase space
* point is about to be generated, so all caches should
* be flushed.
*/
virtual void flushCaches() {
MEGroup::flushCaches();
if ( showerApproximation() )
showerApproximation()->resetBelowCutoff();
}
/**
* Return the shower approximation.
*/
Ptr<ShowerApproximation>::tptr showerApproximation() const;
/**
* Indicate that the shower real emission contribution should be subtracted.
*/
void doRealShowerSubtraction();
/**
* Return true, if the shower real emission contribution should be subtracted.
*/
bool realShowerSubtraction() const { return theRealShowerSubtraction; }
/**
* Indicate that the shower virtual contribution should be subtracted.
*/
void doVirtualShowerSubtraction();
/**
* Return true, if the shower virtual contribution should be subtracted.
*/
bool virtualShowerSubtraction() const { return theVirtualShowerSubtraction; }
//@}
/** @name Matrix element and dipole information */
//@{
/**
* Return the subtraction dipoles.
*/
vector<Ptr<SubtractionDipole>::ptr> dipoles();
/**
* Return the underlying born matrix elements.
*/
const vector<Ptr<MatchboxMEBase>::ptr>& borns() const;
/**
* Build up dipoles needed.
*/
void getDipoles();
/**
* Clone all dipoles.
*/
void cloneDipoles(const string& prefix = "");
/**
* Clone the real emission matrix element.
*/
void cloneRealME(const string& prefix = "");
/**
* Clone all dependencies.
*/
void cloneDependencies(const string& prefix = "") {
cloneDipoles(prefix);
cloneRealME(prefix);
}
/**
* Return all dipoles matching the given Born process
*/
vector<Ptr<SubtractionDipole>::ptr> splitDipoles(const cPDVector&);
//@}
/** @name Veto scale settings */
//@{
/**
* Set veto scales on the particles at the given
* SubProcess which has been generated using this
* matrix element.
*/
virtual void setVetoScales(tSubProPtr) const;
//@}
/** @name Diagnostic information */
//@{
/**
* Dump the setup to an ostream
*/
void print(ostream&) const;
/**
* Collect information on the last evaluated phasespace
* point for verification or debugging purposes. This
* only called, if the StdXCombGroup did accumulate
* a non-zero cross section from this ME group.
*/
virtual void lastEventStatistics();
/**
* Print debug information on the last event
*/
void printLastEvent(ostream&) const;
/**
* Check the subtraction for the last event
*/
void lastEventSubtraction();
/**
* Return true, if verbose
*/
bool verbose() const;
//@}
/** @name Setup of Subtracted ME objects */
//@{
/**
* Return true if this object needs to be initialized before all
* other objects (except those for which this function also returns
* true). This default version always returns false, but subclasses
* may override it to return true.
*/
virtual bool preInitialize() const { return true; }
+ /**
+ * Simple envelope histogram to keep track of subtraction
+ */
+ struct SubtractionHistogram {
+
+ /**
+ * The lower bound
+ */
+ double lower;
+
+ /**
+ * The bins, indexed by upper bound.
+ */
+ map<double,pair<double,double> > bins;
+
+ /**
+ * Constructor
+ */
+ SubtractionHistogram(double low = 0.001,
+ double up = 10.,
+ unsigned int nbins = 100);
+
+ /**
+ * Book an event.
+ */
+ void book(double inv, double ratio) {
+ map<double,pair<double,double> >::iterator b =
+ bins.upper_bound(inv);
+ if ( b == bins.end() ) return;
+ b->second.first = min(b->second.first,abs(ratio));
+ b->second.second = max(b->second.second,abs(ratio));
+ }
+
+ /**
+ * Write to file given name and invariant.
+ */
+ void dump(const std::string& prefix,
+ const cPDVector& proc,
+ int i, int j) const;
+
+ /**
+ * Write to persistent ostream
+ */
+ void persistentOutput(PersistentOStream&) const;
+
+ /**
+ * Read from persistent istream
+ */
+ void persistentInput(PersistentIStream&);
+
+ };
+
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
private:
/**
* The factory which produced this matrix element
*/
Ptr<MatchboxFactory>::tcptr theFactory;
/**
* Pointer to the head real emission ME casted to a MatchboxMEBase
* object.
*/
Ptr<MatchboxMEBase>::ptr theReal;
/**
- * Simple envelope histogram to keep track of subtraction
- */
- struct SubtractionHistogram {
-
- /**
- * The lower bound
- */
- double lower;
-
- /**
- * The bins, indexed by upper bound.
- */
- map<double,pair<double,double> > bins;
-
- /**
- * Constructor
- */
- SubtractionHistogram(double low = 0.001,
- double up = 10.,
- unsigned int nbins = 100);
-
- /**
- * Book an event.
- */
- void book(double inv, double ratio) {
- map<double,pair<double,double> >::iterator b =
- bins.upper_bound(inv);
- if ( b == bins.end() ) return;
- b->second.first = min(b->second.first,abs(ratio));
- b->second.second = max(b->second.second,abs(ratio));
- }
-
- /**
- * Write to file given name and invariant.
- */
- void dump(const std::string& prefix,
- const cPDVector& proc,
- int i, int j) const;
-
- };
-
- /**
* Define the key for the collinear subtraction data.
*/
typedef pair<cPDVector,pair<size_t, size_t> > CollinearSubtractionIndex;
/**
* subtraction data for collinear limits.
*/
map<CollinearSubtractionIndex,SubtractionHistogram> collinearHistograms;
/**
* Define the key for the soft subtraction data.
*/
typedef pair<cPDVector,size_t> SoftSubtractionIndex;
/**
* subtraction data for soft limits.
*/
map<SoftSubtractionIndex,SubtractionHistogram> softHistograms;
/**
* True, if the shower real emission contribution should be subtracted.
*/
bool theRealShowerSubtraction;
/**
* True, if the shower virtual contribution should be subtracted.
*/
bool theVirtualShowerSubtraction;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SubtractedME & operator=(const SubtractedME &);
};
+inline PersistentOStream& operator<<(PersistentOStream& os,
+ const SubtractedME::SubtractionHistogram& h) {
+ h.persistentOutput(os);
+ return os;
+}
+
+inline PersistentIStream& operator>>(PersistentIStream& is,
+ SubtractedME::SubtractionHistogram& h) {
+ h.persistentInput(is);
+ return is;
+}
+
}
#endif /* HERWIG_SubtractedME_H */
diff --git a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
--- a/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
+++ b/MatrixElement/Matchbox/Dipoles/SubtractionDipole.cc
@@ -1,1109 +1,1114 @@
// -*- C++ -*-
//
// SubtractionDipole.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SubtractionDipole class.
//
#include "SubtractionDipole.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDF/PartonBin.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/TildeKinematics.h"
#include "Herwig++/MatrixElement/Matchbox/Phasespace/InvertedTildeKinematics.h"
#include <iterator>
using std::ostream_iterator;
using namespace Herwig;
SubtractionDipole::SubtractionDipole()
: MEBase(), theSplitting(false), theApply(true), theSubtractionTest(false),
theIgnoreCuts(false), theShowerKernel(false),
theRealEmitter(-1), theRealEmission(-1), theRealSpectator(-1),
theBornEmitter(-1), theBornSpectator(-1),
theRealShowerSubtraction(false), theVirtualShowerSubtraction(false),
theRealEmissionScales(false) {}
SubtractionDipole::~SubtractionDipole() {}
void SubtractionDipole::setupBookkeeping() {
theMergingMap.clear();
theSplittingMap.clear();
int xemitter = -1;
int xspectator = -1;
map<int,int> mergeLegs;
map<int,int> remapLegs;
map<int,int> realBornMap;
map<int,int> bornRealMap;
for ( DiagramVector::const_iterator d =
theRealEmissionME->diagrams().begin();
d != theRealEmissionME->diagrams().end(); ++d ) {
Ptr<Tree2toNDiagram>::ptr check =
new_ptr(*dynamic_ptr_cast<Ptr<Tree2toNDiagram>::ptr>(*d));
map<int,int> theMergeLegs;
map<int,int> theRemapLegs;
for ( unsigned int i = 0; i < check->external().size(); ++i )
theMergeLegs[i] = -1;
int theEmitter = check->mergeEmission(realEmitter(),realEmission(),theMergeLegs);
// no underlying Born
if ( theEmitter == -1 )
continue;
DiagramVector::const_iterator bd =
theUnderlyingBornME->diagrams().end();
for ( DiagramVector::const_iterator b =
theUnderlyingBornME->diagrams().begin();
b != theUnderlyingBornME->diagrams().end(); ++b )
if ( check->isSame(*b,theRemapLegs) ) {
bd = b; break;
}
// no underlying Born
if ( bd == theUnderlyingBornME->diagrams().end() )
continue;
if ( xemitter == -1 ) {
xemitter = theEmitter;
mergeLegs = theMergeLegs;
remapLegs = theRemapLegs;
assert(remapLegs.find(xemitter) != remapLegs.end());
xemitter = remapLegs[xemitter];
// work out the leg remapping real -> born
for ( map<int,int>::const_iterator k = mergeLegs.begin();
k != mergeLegs.end(); ++k ) {
assert(remapLegs.find(k->second) != remapLegs.end());
realBornMap[k->first] = remapLegs[k->second];
}
// work out the leg remapping born -> real
for ( map<int,int>::const_iterator k = realBornMap.begin();
k != realBornMap.end(); ++k ) {
bornRealMap[k->second] = k->first;
}
// work out the spectator
assert(mergeLegs.find(realSpectator()) != mergeLegs.end());
assert(remapLegs.find(mergeLegs[realSpectator()]) != remapLegs.end());
xspectator = realBornMap[realSpectator()];
}
RealEmissionKey realKey = realEmissionKey((**d).partons(),realEmitter(),realEmission(),realSpectator());
UnderlyingBornKey bornKey = underlyingBornKey((**bd).partons(),xemitter,xspectator);
if ( theMergingMap.find(realKey) == theMergingMap.end() )
theMergingMap.insert(make_pair(realKey,make_pair(bornKey,realBornMap)));
RealEmissionInfo realInfo = make_pair(realKey,bornRealMap);
bool gotit = false;
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spIterator;
pair<spIterator,spIterator> range = theSplittingMap.equal_range(bornKey);
for ( ; range.first != range.second; ++range.first )
if ( range.first->second == realInfo ) {
gotit = true;
break;
}
if ( !gotit )
theSplittingMap.insert(make_pair(bornKey,realInfo));
}
if ( theSplittingMap.empty() )
return;
theIndexMap.clear();
for ( multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator s =
theSplittingMap.begin(); s != theSplittingMap.end(); ++s ) {
theIndexMap[process(s->first)] = make_pair(emitter(s->first),spectator(s->first));
}
theUnderlyingBornDiagrams.clear();
for ( map<RealEmissionKey,UnderlyingBornInfo>::const_iterator r =
theMergingMap.begin(); r != theMergingMap.end(); ++r ) {
DiagramVector borns;
for ( DiagramVector::const_iterator d = theUnderlyingBornME->diagrams().begin();
d != theUnderlyingBornME->diagrams().end(); ++d ) {
if ( (**d).partons() == process(r->second.first) )
borns.push_back(*d);
}
theUnderlyingBornDiagrams[process(r->first)] = borns;
}
theRealEmissionDiagrams.clear();
for ( multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator r =
theSplittingMap.begin(); r != theSplittingMap.end(); ++r ) {
DiagramVector reals;
for ( DiagramVector::const_iterator d = theRealEmissionME->diagrams().begin();
d != theRealEmissionME->diagrams().end(); ++d ) {
if ( (**d).partons() == process(r->second.first) )
reals.push_back(*d);
}
theRealEmissionDiagrams[process(r->first)] = reals;
}
}
void SubtractionDipole::subtractionBookkeeping() {
if ( theMergingMap.empty() )
setupBookkeeping();
lastRealEmissionKey =
realEmissionKey(lastHeadXComb().mePartonData(),realEmitter(),realEmission(),realSpectator());
map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
theMergingMap.find(lastRealEmissionKey);
if ( k == theMergingMap.end() ) {
theApply = false;
return;
}
theApply = true;
lastUnderlyingBornKey = k->second.first;
bornEmitter(emitter(lastUnderlyingBornKey));
bornSpectator(spectator(lastUnderlyingBornKey));
}
void SubtractionDipole::splittingBookkeeping() {
if ( theMergingMap.empty() )
setupBookkeeping();
map<cPDVector,pair<int,int> >::const_iterator esit =
theIndexMap.find(lastHeadXComb().mePartonData());
if ( esit == theIndexMap.end() ) {
theApply = false;
return;
}
theApply = true;
pair<int,int> es = esit->second;
bornEmitter(es.first);
bornSpectator(es.second);
lastUnderlyingBornKey = underlyingBornKey(lastHeadXComb().mePartonData(),bornEmitter(),bornSpectator());
typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
assert(kr.first != kr.second);
lastRealEmissionInfo = kr.first;
for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
break;
assert(lastRealEmissionInfo != kr.second);
lastRealEmissionKey = lastRealEmissionInfo->second.first;
realEmitter(emitter(lastRealEmissionKey));
realEmission(emission(lastRealEmissionKey));
realSpectator(spectator(lastRealEmissionKey));
}
StdXCombPtr SubtractionDipole::makeXComb(Energy newMaxEnergy, const cPDPair & inc,
tEHPtr newEventHandler,tSubHdlPtr newSubProcessHandler,
tPExtrPtr newExtractor, tCascHdlPtr newCKKW,
const PBPair & newPartonBins, tCutsPtr newCuts,
const DiagramVector & newDiagrams, bool mir,
const PartonPairVec& allBins,
tStdXCombPtr newHead,
tMEPtr newME) {
if ( !newME )
newME = this;
if ( !splitting() ) {
return
underlyingBornME()->makeXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts,
newDiagrams, mir, allBins,
newHead, newME);
}
return
realEmissionME()->makeXComb(newMaxEnergy, inc,
newEventHandler, newSubProcessHandler,
newExtractor, newCKKW,
newPartonBins, newCuts,
newDiagrams, mir, allBins,
newHead, newME);
}
StdXCombPtr SubtractionDipole::makeXComb(tStdXCombPtr newHead,
const PBPair & newPartonBins,
const DiagramVector & newDiagrams,
tMEPtr newME) {
if ( !newME )
newME = this;
if ( !splitting() ) {
return
underlyingBornME()->makeXComb(newHead, newPartonBins,
newDiagrams, newME);
}
return
realEmissionME()->makeXComb(newHead, newPartonBins,
newDiagrams, newME);
}
StdXCombPtr SubtractionDipole::makeBornXComb(tStdXCombPtr realXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*realXC).mePartonData();
lastRealEmissionKey =
realEmissionKey(proc,realEmitter(),realEmission(),realSpectator());
map<RealEmissionKey,UnderlyingBornInfo>::const_iterator k =
theMergingMap.find(lastRealEmissionKey);
if ( k == theMergingMap.end() )
return StdXCombPtr();
PartonPairVec pbs = realXC->pExtractor()->getPartons(realXC->maxEnergy(),
realXC->particles(),
*(realXC->cuts()));
DiagramVector bornDiags = underlyingBornDiagrams(proc);
assert(!bornDiags.empty());
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == bornDiags.front()->partons()[0] &&
ppit->second->parton() == bornDiags.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
return
underlyingBornME()->makeXComb(realXC,*ppit,bornDiags,this);
}
vector<StdXCombPtr> SubtractionDipole::makeRealXCombs(tStdXCombPtr bornXC) {
const cPDVector& proc = const_cast<const StandardXComb&>(*bornXC).mePartonData();
map<cPDVector,pair<int,int> >::const_iterator esit = theIndexMap.find(proc);
if ( esit == theIndexMap.end() )
return vector<StdXCombPtr>();
pair<int,int> es = esit->second;
bornEmitter(es.first);
bornSpectator(es.second);
lastUnderlyingBornKey = underlyingBornKey(proc,bornEmitter(),bornSpectator());
if ( theSplittingMap.find(lastUnderlyingBornKey) == theSplittingMap.end() )
return vector<StdXCombPtr>();
PartonPairVec pbs = bornXC->pExtractor()->getPartons(bornXC->maxEnergy(),
bornXC->particles(),
*(bornXC->cuts()));
DiagramVector realDiags = realEmissionDiagrams(proc);
assert(!realDiags.empty());
vector<StdXCombPtr> res;
map<cPDVector,DiagramVector> realProcs;
for ( MEBase::DiagramVector::const_iterator d = realDiags.begin();
d != realDiags.end(); ++d ) {
realProcs[(**d).partons()].push_back(*d);
}
for ( map<cPDVector,DiagramVector>::const_iterator pr =
realProcs.begin(); pr != realProcs.end(); ++pr ) {
PartonPairVec::iterator ppit = pbs.begin();
for ( ; ppit != pbs.end(); ++ppit ) {
if ( ppit->first->parton() == pr->second.front()->partons()[0] &&
ppit->second->parton() == pr->second.front()->partons()[1] )
break;
}
assert(ppit != pbs.end());
StdXCombPtr rxc =
realEmissionME()->makeXComb(bornXC,*ppit,pr->second,this);
res.push_back(rxc);
}
return res;
}
const MEBase::DiagramVector& SubtractionDipole::underlyingBornDiagrams(const cPDVector& real) const {
static DiagramVector empty;
map<cPDVector,DiagramVector>::const_iterator k = theUnderlyingBornDiagrams.find(real);
if (k == theUnderlyingBornDiagrams.end() )
return empty;
return k->second;
}
const MEBase::DiagramVector& SubtractionDipole::realEmissionDiagrams(const cPDVector& born) const {
static DiagramVector empty;
map<cPDVector,DiagramVector>::const_iterator k = theRealEmissionDiagrams.find(born);
if ( k == theRealEmissionDiagrams.end() )
return empty;
return k->second;
}
void SubtractionDipole::getDiagrams() const {
if ( splitting() ) {
realEmissionME()->diagrams();
useDiagrams(realEmissionME());
} else {
underlyingBornME()->diagrams();
useDiagrams(underlyingBornME());
}
}
Selector<MEBase::DiagramIndex> SubtractionDipole::diagrams(const DiagramVector & dv) const {
Ptr<MatchboxMEBase>::tcptr me =
splitting() ?
realEmissionME() :
underlyingBornME();
if ( me->phasespace() ) {
me->phasespace()->setXComb(lastXCombPtr());
me->phasespace()->fillDiagramWeights();
}
return
me->diagrams(dv);
}
MEBase::DiagramIndex SubtractionDipole::diagram(const DiagramVector & dv) const {
Ptr<MatchboxMEBase>::tcptr me =
splitting() ?
realEmissionME() :
underlyingBornME();
if ( me->phasespace() ) {
me->phasespace()->setXComb(lastXCombPtr());
me->phasespace()->fillDiagramWeights();
}
return
me->diagram(dv);
}
Selector<const ColourLines *>
SubtractionDipole::colourGeometries(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->colourGeometries(diag) :
underlyingBornME()->colourGeometries(diag);
}
const ColourLines &
SubtractionDipole::selectColourGeometry(tcDiagPtr diag) const {
return
splitting() ?
realEmissionME()->selectColourGeometry(diag) :
underlyingBornME()->selectColourGeometry(diag);
}
void SubtractionDipole::flushCaches() {
theUnderlyingBornME->flushCaches();
theRealEmissionME->flushCaches();
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator r =
reweights().begin(); r != reweights().end(); ++r ) {
(**r).flushCaches();
}
}
void SubtractionDipole::setXComb(tStdXCombPtr xc) {
if ( !xc ) {
theApply = false;
return;
} else {
theApply = true;
}
lastMatchboxXComb(xc);
MEBase::setXComb(xc);
if ( splitting() ) {
realEmissionME()->setXComb(xc);
underlyingBornME()->setXComb(xc->head());
splittingBookkeeping();
} else {
realEmissionME()->setXComb(xc->head());
underlyingBornME()->setXComb(xc);
subtractionBookkeeping();
}
if ( !apply() )
return;
}
void SubtractionDipole::setKinematics() {
MEBase::setKinematics();
if ( splitting() )
realEmissionME()->setKinematics();
else
underlyingBornME()->setKinematics();
}
bool SubtractionDipole::generateKinematics(const double * r) {
if ( lastXCombPtr()->kinematicsGenerated() )
return true;
if ( splitting() ) {
if ( !generateRadiationKinematics(r) )
return false;
realEmissionME()->lastXCombPtr()->setIncomingPartons();
realEmissionME()->setScale();
double jac = jacobian();
jac *= pow(underlyingBornME()->lastXComb().lastSHat() / realEmissionME()->lastXComb().lastSHat(),
realEmissionME()->lastXComb().mePartonData().size()-4.);
jacobian(jac);
assert(lastXCombPtr() == realEmissionME()->lastXCombPtr());
lastXCombPtr()->didGenerateKinematics();
return true;
}
if ( !generateTildeKinematics() )
return false;
underlyingBornME()->setScale();
assert(lastXCombPtr() == underlyingBornME()->lastXCombPtr());
underlyingBornME()->lastXCombPtr()->setIncomingPartons();
lastXCombPtr()->didGenerateKinematics();
return true;
}
int SubtractionDipole::nDim() const {
if ( !splitting() )
return underlyingBornME()->nDim();
return underlyingBornME()->nDim() + nDimRadiation();
}
void SubtractionDipole::clearKinematics() {
MEBase::clearKinematics();
if ( splitting() )
realEmissionME()->clearKinematics();
else
underlyingBornME()->clearKinematics();
}
void SubtractionDipole::tildeKinematics(Ptr<TildeKinematics>::tptr tk) {
theTildeKinematics = tk;
}
bool SubtractionDipole::generateTildeKinematics() {
assert(!splitting());
if ( !tildeKinematics() ) {
jacobian(0.0);
return false;
}
theTildeKinematics->prepare(lastHeadXCombPtr(),lastXCombPtr());
if ( !theTildeKinematics->doMap() ) {
jacobian(0.0);
return false;
}
theLastSubtractionScale = theTildeKinematics->lastScale();
theLastSubtractionPt = theTildeKinematics->lastPt();
meMomenta().resize(lastHeadXComb().meMomenta().size() - 1);
assert(mergingMap().find(lastRealEmissionKey) != mergingMap().end());
map<int,int>& trans = theMergingMap[lastRealEmissionKey].second;
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == realEmitter() || k == realEmission() || k == realSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( theTildeKinematics->doesTransform() && k > 1 )
meMomenta()[trans[k]] = theTildeKinematics->transform(meMomenta()[trans[k]]);
}
meMomenta()[bornEmitter()] = tildeKinematics()->bornEmitterMomentum();
meMomenta()[bornSpectator()] = tildeKinematics()->bornSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).mass());
p->rescaleRho();
}
jacobian(realEmissionME()->lastXComb().jacobian());
logGenerateTildeKinematics();
return true;
}
void SubtractionDipole::invertedTildeKinematics(Ptr<InvertedTildeKinematics>::tptr itk) {
theInvertedTildeKinematics = itk;
}
int SubtractionDipole::nDimRadiation() const {
return invertedTildeKinematics() ?
invertedTildeKinematics()->nDimRadiation() :
0;
}
bool SubtractionDipole::generateRadiationKinematics(const double * r) {
assert(splitting());
if ( !invertedTildeKinematics() ) {
jacobian(0.0);
return false;
}
theInvertedTildeKinematics->prepare(lastXCombPtr(),lastHeadXCombPtr());
if ( !theInvertedTildeKinematics->doMap(r) ) {
jacobian(0.0);
return false;
}
theLastSplittingScale = theInvertedTildeKinematics->lastScale();
theLastSplittingPt = theInvertedTildeKinematics->lastPt();
meMomenta().resize(lastHeadXComb().meMomenta().size() + 1);
assert(splittingMap().find(lastUnderlyingBornKey) != splittingMap().end());
map<int,int>& trans = const_cast<map<int,int>&>(lastRealEmissionInfo->second.second);
int n = lastHeadXComb().meMomenta().size();
for ( int k = 0; k < n; ++k ) {
if ( k == bornEmitter() || k == bornSpectator() )
continue;
meMomenta()[trans[k]] = lastHeadXComb().meMomenta()[k];
if ( invertedTildeKinematics()->doesTransform() && k > 1 )
meMomenta()[trans[k]] = invertedTildeKinematics()->transform(meMomenta()[trans[k]]);
}
meMomenta()[realEmitter()] = invertedTildeKinematics()->realEmitterMomentum();
meMomenta()[realEmission()] = invertedTildeKinematics()->realEmissionMomentum();
meMomenta()[realSpectator()] = invertedTildeKinematics()->realSpectatorMomentum();
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = meMomenta().begin();
for ( ; pd != mePartonData().end(); ++pd, ++p ) {
p->setMass((**pd).mass());
p->rescaleRho();
}
jacobian(underlyingBornME()->lastXComb().jacobian() *
invertedTildeKinematics()->jacobian());
logGenerateRadiationKinematics(r);
return true;
}
void SubtractionDipole::ptCut(Energy cut) {
theInvertedTildeKinematics->ptCut(cut);
}
CrossSection SubtractionDipole::dSigHatDR(Energy2 factorizationScale) const {
double pdfweight = 1.;
double jac = jacobian();
if ( splitting() && jac == 0.0 ) {
lastMECrossSection(ZERO);
lastME2(0.0);
return ZERO;
}
if ( factorizationScale == ZERO ) {
if ( realEmissionScales() )
factorizationScale = realEmissionME()->lastScale();
else
factorizationScale = underlyingBornME()->lastScale();
}
if ( havePDFWeight1() ) {
pdfweight *= realEmissionME()->pdf1(factorizationScale);
}
if ( havePDFWeight2() ) {
pdfweight *= realEmissionME()->pdf2(factorizationScale);
}
lastMEPDFWeight(pdfweight);
if ( showerApproximation() && realShowerSubtraction() ) {
assert(!splitting());
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(this);
if ( !showerApproximation()->isAboveCutoff() )
showerApproximation()->wasBelowCutoff();
if ( !showerApproximation()->isInShowerPhasespace() ) {
lastMECrossSection(ZERO);
lastME2(0.0);
return ZERO;
}
lastMECrossSection(-showerApproximation()->dSigHatDR());
return lastMECrossSection();
}
double xme2 = 0.0;
if ( lastME2() == 0.0 ) {
if ( !showerKernel() )
xme2 = me2();
else
xme2 = me2Avg(-underlyingBornME()->me2());
} else {
xme2 = lastME2();
}
if ( xme2 == 0.0 ) {
lastMECrossSection(ZERO);
lastME2(0.0);
return ZERO;
}
if ( realEmissionScales() ) {
xme2 *=
pow(realEmissionME()->lastXComb().lastAlphaS()/
underlyingBornME()->lastXComb().lastAlphaS(),
realEmissionME()->orderInAlphaS());
xme2 *=
pow(realEmissionME()->lastXComb().lastAlphaEM()/
underlyingBornME()->lastXComb().lastAlphaEM(),
underlyingBornME()->orderInAlphaEW());
}
lastME2(xme2);
CrossSection res =
sqr(hbarc) * jac * pdfweight * xme2 /
(2. * realEmissionME()->lastXComb().lastSHat());
if ( !showerApproximation() ) {
double weight = 0.0;
bool applied = false;
for ( vector<Ptr<MatchboxReweightBase>::ptr>::const_iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
(**rw).setXComb(theRealEmissionME->lastXCombPtr());
if ( !(**rw).apply() )
continue;
weight += (**rw).evaluate();
applied = true;
}
if ( applied )
res *= weight;
}
if ( showerApproximation() && virtualShowerSubtraction() ) {
assert(!splitting());
showerApproximation()->setBornXComb(lastXCombPtr());
showerApproximation()->setRealXComb(realEmissionME()->lastXCombPtr());
showerApproximation()->setDipole(this);
if ( !showerApproximation()->isAboveCutoff() )
showerApproximation()->wasBelowCutoff();
else
lastME2(0.0);
CrossSection shower = ZERO;
if ( showerApproximation()->isInShowerPhasespace() )
shower = showerApproximation()->dSigHatDR();
res -= shower;
}
lastMECrossSection(-res);
logDSigHatDR(jac);
return lastMECrossSection();
}
void SubtractionDipole::print(ostream& os) const {
os << "--- SubtractionDipole setup ----------------------------------------------------\n";
os << " subtraction '" << name() << "'\n for real emission '"
<< theRealEmissionME->name() << "'\n using underlying Born '"
<< theUnderlyingBornME->name() << "'\n";
os << " tilde kinematics are '"
<< (theTildeKinematics ? theTildeKinematics->name() : "")
<< " '\n inverted tilde kinematics are '"
<< (theInvertedTildeKinematics ? theInvertedTildeKinematics->name() : "") << "'\n";
os << " the following subtraction mappings have been found:\n";
for ( map<RealEmissionKey,UnderlyingBornInfo>::const_iterator m =
theMergingMap.begin(); m != theMergingMap.end(); ++m ) {
os << " " << process(m->second.first)[0]->PDGName() << " "
<< process(m->second.first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->second.first).begin() + 2;
p != process(m->second.first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[" << emitter(m->second.first) << "," << spectator(m->second.first) << "] <=> ";
os << process(m->first)[0]->PDGName() << " "
<< process(m->first)[1]->PDGName() << " -> ";
for ( cPDVector::const_iterator p = process(m->first).begin() + 2;
p != process(m->first).end(); ++p ) {
os << (**p).PDGName() << " ";
}
os << "[(" << emitter(m->first) << "," << emission(m->first) << ")," << spectator(m->first) << "]\n"
<< " non-dipole momenta ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->second << " ";
}
os << ") <=> ( ";
for ( map<int,int>::const_iterator k = m->second.second.begin();
k != m->second.second.end(); ++k ) {
if ( k->first == spectator(m->first) )
continue;
os << k->first << " ";
}
os << ")\n";
}
os << "--------------------------------------------------------------------------------\n";
os << flush;
}
void SubtractionDipole::printLastEvent(ostream& os) const {
os << "--- SubtractionDipole last event information -----------------------------------\n";
os << " for dipole '" << name() << "' applying ["
<< bornEmitter() << "," << bornSpectator() << "] <=> [("
<< realEmitter() << "," << realEmission() << ")," << realSpectator() << "]\n"
<< " evaluated the cross section/nb " << (lastMECrossSection()/nanobarn) << "\n"
<< " with subtraction parameters x[0] = " << subtractionParameters()[0]
<< " x[1] = " << subtractionParameters()[1] << "\n";
os << " the last real emission event was:\n";
realEmissionME()->printLastEvent(os);
os << " the last underlying Born event was:\n";
underlyingBornME()->printLastEvent(os);
os << "--- end SubtractionDipole last event information -------------------------------\n";
os << flush;
}
void SubtractionDipole::logME2() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated me2 using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "Born phase space point (in GeV):\n";
vector<Lorentz5Momentum>::const_iterator pit = bornxc->meMomenta().begin();
cPDVector::const_iterator dit = bornxc->mePartonData().begin();
for ( ; pit != bornxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << bornxc->lastX1() << " x2 = " << bornxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (bornxc->lastSHat()/GeV2) << "\n";
generator()->log() << "Real emission phase space point (in GeV):\n";
pit = realxc->meMomenta().begin();
dit = realxc->mePartonData().begin();
for ( ; pit != realxc->meMomenta().end() ; ++pit, ++dit )
generator()->log() << (**dit).PDGName() << " : "
<< (*pit/GeV) << "\n";
generator()->log() << "with x1 = " << realxc->lastX1() << " x2 = " << realxc->lastX2() << "\n"
<< "sHat/GeV2 = " << (realxc->lastSHat()/GeV2) << "\n";
generator()->log() << "me2 = " << lastME2() << "\n" << flush;
}
void SubtractionDipole::logDSigHatDR(double effectiveJac) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
tcStdXCombPtr bornxc = splitting() ? lastHeadXCombPtr() : lastXCombPtr();
tcStdXCombPtr realxc = splitting() ? lastXCombPtr() : lastHeadXCombPtr();
generator()->log() << "'" << name() << "' evaluated cross section using\n"
<< "Born XComb " << bornxc << " real XComb " << realxc << "\n"
<< "Jacobian = " << jacobian()
<< " effective Jacobian = " << effectiveJac << "\n"
<< "Born sHat/GeV2 = " << (bornxc->lastSHat()/GeV2)
<< " real sHat/GeV2 = " << (realxc->lastSHat()/GeV2)
<< " dsig/nb = "
<< (lastMECrossSection()/nanobarn) << "\n" << flush;
}
void SubtractionDipole::logGenerateTildeKinematics() const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating tilde kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with real xcomb " << lastHeadXCombPtr() << " born xcomb "
<< lastXCombPtr() << "\n"
<< "from real emission phase space point:\n";
Lorentz5Momentum rSum;
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
size_t count = 0;
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
rSum -= *pr;
} else {
rSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (rSum/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n"
<< "with scale/GeV = " << (theLastSubtractionScale/GeV)
<< "and pt/GeV = " << (theLastSubtractionPt/GeV) << "\n";
generator()->log() << "generated tilde kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
count = 0;
Lorentz5Momentum bSum;
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr ) {
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
if ( count < 2 ) {
bSum -= *pr;
} else {
bSum += *pr;
}
++count;
}
generator()->log() << "sum : " << (bSum/GeV) << "\n";
generator()->log() << "Jacobian = " << jacobian() << "\n" << flush;
}
void SubtractionDipole::logGenerateRadiationKinematics(const double * r) const {
if ( !realEmissionME()->verbose() &&
!underlyingBornME()->verbose() )
return;
generator()->log() << "'" << name() << "' generating radiation kinematics.\n"
<< "configuration: [" << bornEmitter() << ","
<< bornSpectator() << "] => "
<< "[(" << realEmitter() << "," << realEmission() << "),"
<< realSpectator() << "]\n"
<< "with born xcomb " << lastHeadXCombPtr() << " real xcomb "
<< lastXCombPtr() << "\n"
<< "from random numbers:\n";
copy(r,r+nDimRadiation(),ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n";
generator()->log() << "and born phase space point:\n";
vector<Lorentz5Momentum>::const_iterator pr = lastHeadXComb().meMomenta().begin();
cPDVector::const_iterator dr = lastHeadXComb().mePartonData().begin();
for ( ; pr != lastHeadXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "subtraction parameters: ";
copy(subtractionParameters().begin(),subtractionParameters().end(),
ostream_iterator<double>(generator()->log()," "));
generator()->log() << "\n" << flush;
generator()->log() << "scales: scale/GeV = " << (theLastSplittingScale/GeV)
<< " pt/GeV = " << (theLastSplittingPt/GeV) << "\n" << flush;
generator()->log() << "generated real emission kinematics:\n";
pr = lastXComb().meMomenta().begin();
dr = lastXComb().mePartonData().begin();
for ( ; pr != lastXComb().meMomenta().end(); ++pr,++dr )
generator()->log() << (**dr).PDGName() << " : "
<< (*pr/GeV) << "\n";
generator()->log() << "Jacobian = "
<< jacobian() << " = "
<< underlyingBornME()->lastXComb().jacobian()
<< "|Born * "
<< invertedTildeKinematics()->jacobian()
<< "|Radiation\n" << flush;
}
void SubtractionDipole::cloneDependencies(const std::string& prefix) {
if ( underlyingBornME() ) {
Ptr<MatchboxMEBase>::ptr myUnderlyingBornME = underlyingBornME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myUnderlyingBornME->name();
if ( ! (generator()->preinitRegister(myUnderlyingBornME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
myUnderlyingBornME->cloneDependencies(pname.str());
underlyingBornME(myUnderlyingBornME);
}
if ( realEmissionME() ) {
Ptr<MatchboxMEBase>::ptr myRealEmissionME = realEmissionME()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myRealEmissionME->name();
if ( ! (generator()->preinitRegister(myRealEmissionME,pname.str()) ) )
throw InitException() << "Matrix element " << pname.str() << " already existing.";
myRealEmissionME->cloneDependencies(pname.str());
realEmissionME(myRealEmissionME);
}
if ( tildeKinematics() ) {
Ptr<TildeKinematics>::ptr myTildeKinematics = tildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myTildeKinematics->name();
if ( ! (generator()->preinitRegister(myTildeKinematics,pname.str()) ) )
throw InitException() << "Tilde kinematics " << pname.str() << " already existing.";
myTildeKinematics->dipole(this);
tildeKinematics(myTildeKinematics);
}
if ( invertedTildeKinematics() ) {
Ptr<InvertedTildeKinematics>::ptr myInvertedTildeKinematics = invertedTildeKinematics()->cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << myInvertedTildeKinematics->name();
if ( ! (generator()->preinitRegister(myInvertedTildeKinematics,pname.str()) ) )
throw InitException() << "Inverted tilde kinematics " << pname.str() << " already existing.";
myInvertedTildeKinematics->dipole(this);
invertedTildeKinematics(myInvertedTildeKinematics);
}
for ( vector<Ptr<MatchboxReweightBase>::ptr>::iterator rw =
theReweights.begin(); rw != theReweights.end(); ++rw ) {
Ptr<MatchboxReweightBase>::ptr myReweight = (**rw).cloneMe();
ostringstream pname;
pname << (prefix == "" ? fullName() : prefix) << "/" << (**rw).name();
if ( ! (generator()->preinitRegister(myReweight,pname.str()) ) )
throw InitException() << "Reweight " << pname.str() << " already existing.";
myReweight->cloneDependencies(pname.str());
*rw = myReweight;
}
}
void SubtractionDipole::constructVertex(tSubProPtr sub) {
if ( splitting() )
realEmissionME()->constructVertex(sub);
else
underlyingBornME()->constructVertex(sub);
}
void SubtractionDipole::generateSubCollision(SubProcess & sub) {
if ( splitting() )
realEmissionME()->generateSubCollision(sub);
else
underlyingBornME()->generateSubCollision(sub);
}
void SubtractionDipole::persistentOutput(PersistentOStream & os) const {
- os << theSplitting << theApply << theSubtractionTest << theIgnoreCuts
- << theShowerKernel << theRealEmissionME << theUnderlyingBornME
- << theTildeKinematics << theInvertedTildeKinematics << theReweights
- << theRealEmitter << theRealEmission << theRealSpectator
- << theSubtractionParameters
- << theMergingMap << theSplittingMap << theIndexMap
- << theUnderlyingBornDiagrams << theRealEmissionDiagrams
- << lastRealEmissionKey << lastUnderlyingBornKey
- << theBornEmitter << theBornSpectator
- << theShowerApproximation
- << theRealShowerSubtraction << theVirtualShowerSubtraction
- << thePartners << theRealEmissionScales
- << theLastXComb;
+ os << theLastXComb << theSplitting << theApply << theSubtractionTest
+ << theIgnoreCuts << theShowerKernel << theRealEmissionME << theUnderlyingBornME
+ << thePartners << theTildeKinematics << theInvertedTildeKinematics
+ << theReweights << theRealEmitter << theRealEmission << theRealSpectator
+ << theSubtractionParameters << theMergingMap << theSplittingMap
+ << theIndexMap << theUnderlyingBornDiagrams << theRealEmissionDiagrams
+ << lastRealEmissionKey << lastUnderlyingBornKey
+ << theBornEmitter << theBornSpectator << ounit(theLastSubtractionScale,GeV)
+ << ounit(theLastSplittingScale,GeV) << ounit(theLastSubtractionPt,GeV)
+ << ounit(theLastSplittingPt,GeV) << theShowerApproximation
+ << theRealShowerSubtraction << theVirtualShowerSubtraction
+ << theRealEmissionScales;
}
void SubtractionDipole::persistentInput(PersistentIStream & is, int) {
- is >> theSplitting >> theApply >> theSubtractionTest >> theIgnoreCuts
- >> theShowerKernel >> theRealEmissionME >> theUnderlyingBornME
- >> theTildeKinematics >> theInvertedTildeKinematics >> theReweights
- >> theRealEmitter >> theRealEmission >> theRealSpectator
- >> theSubtractionParameters
- >> theMergingMap >> theSplittingMap >> theIndexMap
- >> theUnderlyingBornDiagrams >> theRealEmissionDiagrams
- >> lastRealEmissionKey >> lastUnderlyingBornKey
- >> theBornEmitter >> theBornSpectator
- >> theShowerApproximation
- >> theRealShowerSubtraction >> theVirtualShowerSubtraction
- >> thePartners >> theRealEmissionScales
- >> theLastXComb;
+ is >> theLastXComb >> theSplitting >> theApply >> theSubtractionTest
+ >> theIgnoreCuts >> theShowerKernel >> theRealEmissionME >> theUnderlyingBornME
+ >> thePartners >> theTildeKinematics >> theInvertedTildeKinematics
+ >> theReweights >> theRealEmitter >> theRealEmission >> theRealSpectator
+ >> theSubtractionParameters >> theMergingMap >> theSplittingMap
+ >> theIndexMap >> theUnderlyingBornDiagrams >> theRealEmissionDiagrams
+ >> lastRealEmissionKey >> lastUnderlyingBornKey
+ >> theBornEmitter >> theBornSpectator >> iunit(theLastSubtractionScale,GeV)
+ >> iunit(theLastSplittingScale,GeV) >> iunit(theLastSubtractionPt,GeV)
+ >> iunit(theLastSplittingPt,GeV) >> theShowerApproximation
+ >> theRealShowerSubtraction >> theVirtualShowerSubtraction
+ >> theRealEmissionScales;
lastMatchboxXComb(theLastXComb);
+ typedef multimap<UnderlyingBornKey,RealEmissionInfo>::const_iterator spit;
+ pair<spit,spit> kr = theSplittingMap.equal_range(lastUnderlyingBornKey);
+ assert(kr.first != kr.second);
+ lastRealEmissionInfo = kr.first;
+ for ( ; lastRealEmissionInfo != kr.second; ++lastRealEmissionInfo )
+ if ( process(lastRealEmissionInfo->second.first) == lastXComb().mePartonData() )
+ break;
}
void SubtractionDipole::Init() {
static ClassDocumentation<SubtractionDipole> documentation
("SubtractionDipole represents a dipole subtraction "
"term in the formalism of Catani and Seymour.");
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<SubtractionDipole,MEBase>
describeSubtractionDipole("Herwig::SubtractionDipole", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc b/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc
--- a/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc
+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationGenerator.cc
@@ -1,400 +1,404 @@
// -*- C++ -*-
//
// ShowerApproximationGenerator.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ShowerApproximationGenerator class.
//
#include "ShowerApproximationGenerator.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/PDF/PartonExtractor.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
ShowerApproximationGenerator::ShowerApproximationGenerator()
: thePresamplingPoints(10000), theMaxTry(100000) {}
ShowerApproximationGenerator::~ShowerApproximationGenerator() {}
double ShowerApproximationGenerator::generateFraction(tcPDPtr pd, double r, double xmin) const {
if ( pd->coloured() || pd->id() == ParticleID::gamma ) {
return pow(xmin,r);
}
double x0 = 1.e-5;
return 1. + x0 - x0*pow((1.+x0)/x0,r);
}
double ShowerApproximationGenerator::invertFraction(tcPDPtr pd, double x, double xmin) const {
if ( pd->coloured() || pd->id() == ParticleID::gamma ) {
return log(x)/log(xmin);
}
double x0 = 1.e-5;
return log((1.-x+x0)/x0)/log((1.+x0)/x0);
}
bool ShowerApproximationGenerator::prepare() {
tSubProPtr oldSub = lastIncomingXComb->subProcess();
tcStdXCombPtr cIncomingXC = lastIncomingXComb;
theLastMomenta.resize(cIncomingXC->mePartonData().size());
theLastRandomNumbers.resize(thePhasespace->nDim(theLastMomenta.size()-2) + 2);
double x1 =
oldSub->incoming().first->momentum().plus() /
lastIncomingXComb->lastParticles().first->momentum().plus();
theLastRandomNumbers[0] = invertFraction(oldSub->incoming().first->dataPtr(),x1,
lastIncomingXComb->cuts()->x1Min());
double x2 =
oldSub->incoming().second->momentum().minus() /
lastIncomingXComb->lastParticles().second->momentum().minus();
theLastRandomNumbers[1] = invertFraction(oldSub->incoming().second->dataPtr(),x2,
lastIncomingXComb->cuts()->x2Min());
Boost toCMS =
(oldSub->incoming().first->momentum() +
oldSub->incoming().second->momentum()).findBoostToCM();
theLastMomenta[0] = oldSub->incoming().first->momentum();
theLastMomenta[0].boost(toCMS);
theLastMomenta[1] = oldSub->incoming().second->momentum();
theLastMomenta[1].boost(toCMS);
ParticleVector::const_iterator out = oldSub->outgoing().begin();
vector<Lorentz5Momentum>::iterator p = theLastMomenta.begin() + 2;
for ( ; out != oldSub->outgoing().end(); ++out, ++p ) {
*p = (**out).momentum();
p->boost(toCMS);
}
theLastPartons.first =
oldSub->incoming().first->data().produceParticle(oldSub->incoming().first->momentum());
theLastPartons.second =
oldSub->incoming().second->data().produceParticle(oldSub->incoming().second->momentum());
thePhasespace->invertKinematics(theLastMomenta,&theLastRandomNumbers[2]);
theLastBornXComb->clean();
theLastBornXComb->fill(lastIncomingXComb->lastParticles(),theLastPartons,
theLastMomenta,theLastRandomNumbers);
if ( !theLastBornXComb->cuts()->initSubProcess(theLastBornXComb->lastSHat(),
theLastBornXComb->lastY(),
theLastBornXComb->mirror()) )
return false;
theLastBornME->setXComb(theLastBornXComb);
if ( !theLastBornME->generateKinematics(&theLastRandomNumbers[2]) )
return false;
CrossSection bornXS = theLastBornME->dSigHatDR();
if ( bornXS == ZERO )
return false;
return true;
}
bool ShowerApproximationGenerator::generate(const vector<double>& r) {
theLastBornXComb->clean();
double x = generateFraction(theLastPartons.first->dataPtr(),r[0],
lastIncomingXComb->cuts()->x1Min());
Energy Q = lastIncomingXComb->lastParticles().first->momentum().plus();
Energy mass = theLastPartons.first->dataPtr()->mass();
double xi = (4.*sqr(x*Q) - sqr(mass))/(4.*sqr(Q)*x);
Lorentz5Momentum p1(ZERO,ZERO,xi*Q);
p1.setMass(mass); p1.rescaleEnergy();
theLastPartons.first->set5Momentum(p1);
x = generateFraction(theLastPartons.second->dataPtr(),r[1],
lastIncomingXComb->cuts()->x1Min());
Q = lastIncomingXComb->lastParticles().second->momentum().minus();
mass = theLastPartons.second->dataPtr()->mass();
xi = (4.*sqr(x*Q) - sqr(mass))/(4.*sqr(Q)*x);
Lorentz5Momentum p2(ZERO,ZERO,-xi*Q);
p2.setMass(mass); p2.rescaleEnergy();
theLastPartons.second->set5Momentum(p2);
theLastPresamplingMomenta.resize(theLastMomenta.size());
Boost toCMS =
(theLastPartons.first->momentum() +
theLastPartons.second->momentum()).findBoostToCM();
theLastPresamplingMomenta[0] = theLastPartons.first->momentum();
theLastPresamplingMomenta[0].boost(toCMS);
theLastPresamplingMomenta[1] = theLastPartons.second->momentum();
theLastPresamplingMomenta[1].boost(toCMS);
theLastBornXComb->fill(lastIncomingXComb->lastParticles(),theLastPartons,
theLastPresamplingMomenta,r);
if ( !theLastBornXComb->cuts()->initSubProcess(theLastBornXComb->lastSHat(),
theLastBornXComb->lastY(),
theLastBornXComb->mirror()) )
return false;
theLastBornME->setXComb(theLastBornXComb);
if ( !theLastBornME->generateKinematics(&r[2]) )
return false;
CrossSection bornXS = theLastBornME->dSigHatDR();
if ( bornXS == ZERO )
return false;
return true;
}
void ShowerApproximationGenerator::restore() {
tSubProPtr oldSub = lastIncomingXComb->subProcess();
theLastPartons.first->set5Momentum(oldSub->incoming().first->momentum());
theLastPartons.second->set5Momentum(oldSub->incoming().second->momentum());
theLastBornXComb->clean();
theLastBornXComb->fill(lastIncomingXComb->lastParticles(),theLastPartons,
theLastMomenta,theLastRandomNumbers);
theLastBornME->setXComb(theLastBornXComb);
theLastBornME->generateKinematics(&theLastRandomNumbers[2]);
theLastBornME->dSigHatDR();
}
void ShowerApproximationGenerator::
handle(EventHandler & eh, const tPVector &,
const Hint &) {
lastIncomingXComb = dynamic_ptr_cast<tStdXCombPtr>(eh.lastXCombPtr());
if ( !lastIncomingXComb )
throw Exception() << "expecting a standard event handler"
<< Exception::abortnow;
if ( lastIncomingXComb->lastProjector() )
lastIncomingXComb = lastIncomingXComb->lastProjector();
const StandardXComb& xc = *lastIncomingXComb;
map<cPDVector,set<Ptr<ShowerApproximationKernel>::ptr> >::const_iterator
kernelit = theKernelMap.find(xc.mePartonData());
if ( kernelit == theKernelMap.end() ) {
list<MatchboxFactory::SplittingChannel> channels =
theFactory->getSplittingChannels(lastIncomingXComb);
set<Ptr<ShowerApproximationKernel>::ptr> newKernels;
for ( list<MatchboxFactory::SplittingChannel>::const_iterator c =
channels.begin(); c != channels.end(); ++c ) {
Ptr<ShowerApproximationKernel>::ptr kernel =
new_ptr(ShowerApproximationKernel());
kernel->setBornXComb(c->bornXComb);
kernel->setRealXComb(c->realXComb);
kernel->setTildeXCombs(c->tildeXCombs);
kernel->setDipole(c->dipole);
kernel->showerApproximation(theShowerApproximation);
kernel->presamplingPoints(thePresamplingPoints);
kernel->maxtry(theMaxTry);
kernel->showerApproximationGenerator(this);
if ( kernel->dipole()->bornEmitter() > 1 &&
kernel->dipole()->bornSpectator() > 1 ) {
kernel->ptCut(ffPtCut());
} else if ( ( kernel->dipole()->bornEmitter() > 1 &&
kernel->dipole()->bornSpectator() < 2 ) ||
( kernel->dipole()->bornEmitter() < 2 &&
kernel->dipole()->bornSpectator() > 1 ) ) {
kernel->ptCut(fiPtCut());
} else {
assert(kernel->dipole()->bornEmitter() < 2 &&
kernel->dipole()->bornSpectator() < 2);
kernel->ptCut(iiPtCut());
}
newKernels.insert(kernel);
}
theKernelMap[xc.mePartonData()] = newKernels;
kernelit = theKernelMap.find(xc.mePartonData());
}
if ( kernelit->second.empty() )
return;
const set<Ptr<ShowerApproximationKernel>::ptr>& kernels = kernelit->second;
theLastBornME = (**kernels.begin()).dipole()->underlyingBornME();
theLastBornME->phasespace(thePhasespace);
theLastBornXComb = (**kernels.begin()).bornXComb();
if ( !prepare() )
return;
Energy winnerPt = ZERO;
Ptr<ShowerApproximationKernel>::ptr winnerKernel;
for ( set<Ptr<ShowerApproximationKernel>::ptr>::const_iterator k =
kernels.begin(); k != kernels.end(); ++k ) {
if ( (**k).generate() != 0. ) {
winnerKernel = *k;
winnerPt = max(winnerPt,winnerKernel->dipole()->lastPt());
}
}
if ( !winnerKernel || winnerPt == ZERO )
return;
SubProPtr oldSub = lastIncomingXComb->subProcess();
SubProPtr newSub;
try {
newSub = winnerKernel->realXComb()->construct();
} catch(Veto&) {
return;
}
tParticleSet firstS = oldSub->incoming().first->siblings();
assert(firstS.empty() || firstS.size() == 1);
if ( !firstS.empty() ) {
eh.currentStep()->removeParticle(*firstS.begin());
}
tParticleSet secondS = oldSub->incoming().second->siblings();
assert(secondS.empty() || secondS.size() == 1);
if ( !secondS.empty() ) {
eh.currentStep()->removeParticle(*secondS.begin());
}
// prevent the colliding particles from disappearing
// in the initial state and appearing
// in the final state when we've cut off all their
// (physical) children; only applies to the case
// where we have a parton extractor not build from
// noPDF, so check wether the incoming particle
// doesnt equal the incoming parton -- this needs fixing in ThePEG
PPtr dummy = new_ptr(Particle(getParticleData(ParticleID::gamma)));
bool usedDummy = false;
if ( eh.currentStep()->incoming().first != oldSub->incoming().first ) {
eh.currentStep()->addDecayProduct(eh.currentStep()->incoming().first,dummy);
usedDummy = true;
}
if ( eh.currentStep()->incoming().second != oldSub->incoming().second ) {
eh.currentStep()->addDecayProduct(eh.currentStep()->incoming().second,dummy);
usedDummy = true;
}
eh.currentStep()->removeSubProcess(oldSub);
eh.currentStep()->addSubProcess(newSub);
// get rid of the dummy
if ( usedDummy ) {
eh.currentStep()->removeParticle(dummy);
}
eh.select(winnerKernel->realXComb());
winnerKernel->realXComb()->recreatePartonBinInstances(winnerKernel->realXComb()->lastScale());
winnerKernel->realXComb()->refillPartonBinInstances(&(xc.lastRandomNumbers()[0]));
winnerKernel->realXComb()->pExtractor()->constructRemnants(winnerKernel->realXComb()->partonBinInstances(),
newSub, eh.currentStep());
}
IBPtr ShowerApproximationGenerator::clone() const {
return new_ptr(*this);
}
IBPtr ShowerApproximationGenerator::fullclone() const {
return new_ptr(*this);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void ShowerApproximationGenerator::persistentOutput(PersistentOStream & os) const {
- os << theShowerApproximation << thePhasespace
- << theFactory << theKernelMap
- << thePresamplingPoints << theMaxTry;
+ os << theShowerApproximation << thePhasespace << theFactory
+ << theKernelMap << thePresamplingPoints << theMaxTry
+ << lastIncomingXComb << theLastBornME << ounit(theLastMomenta,GeV)
+ << ounit(theLastPresamplingMomenta,GeV) << theLastRandomNumbers
+ << theLastBornXComb << theLastPartons;
}
void ShowerApproximationGenerator::persistentInput(PersistentIStream & is, int) {
- is >> theShowerApproximation >> thePhasespace
- >> theFactory >> theKernelMap
- >> thePresamplingPoints >> theMaxTry;
+ is >> theShowerApproximation >> thePhasespace >> theFactory
+ >> theKernelMap >> thePresamplingPoints >> theMaxTry
+ >> lastIncomingXComb >> theLastBornME >> iunit(theLastMomenta,GeV)
+ >> iunit(theLastPresamplingMomenta,GeV) >> theLastRandomNumbers
+ >> theLastBornXComb >> theLastPartons;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<ShowerApproximationGenerator,StepHandler>
describeHerwigShowerApproximationGenerator("Herwig::ShowerApproximationGenerator", "HwMatchbox.so");
void ShowerApproximationGenerator::Init() {
static ClassDocumentation<ShowerApproximationGenerator> documentation
("ShowerApproximationGenerator generates emissions according to a "
"shower approximation entering a NLO matching.");
static Reference<ShowerApproximationGenerator,ShowerApproximation> interfaceShowerApproximation
("ShowerApproximation",
"Set the shower approximation to sample.",
&ShowerApproximationGenerator::theShowerApproximation, false, false, true, false, false);
static Reference<ShowerApproximationGenerator,MatchboxPhasespace> interfacePhasespace
("Phasespace",
"The phase space generator to use.",
&ShowerApproximationGenerator::thePhasespace, false, false, true, false, false);
static Reference<ShowerApproximationGenerator,MatchboxFactory> interfaceFactory
("Factory",
"The factory object to use.",
&ShowerApproximationGenerator::theFactory, false, false, true, false, false);
static Parameter<ShowerApproximationGenerator,unsigned long> interfacePresamplingPoints
("PresamplingPoints",
"Set the number of presampling points.",
&ShowerApproximationGenerator::thePresamplingPoints, 10000, 1, 0,
false, false, Interface::lowerlim);
static Parameter<ShowerApproximationGenerator,unsigned long> interfaceMaxTry
("MaxTry",
"Set the number of maximum attempts.",
&ShowerApproximationGenerator::theMaxTry, 100000, 1, 0,
false, false, Interface::lowerlim);
}
diff --git a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc
--- a/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc
+++ b/MatrixElement/Matchbox/Matching/ShowerApproximationKernel.cc
@@ -1,190 +1,201 @@
// -*- C++ -*-
//
// ShowerApproximationKernel.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ShowerApproximationKernel class.
//
#include "ShowerApproximationKernel.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ShowerApproximationGenerator.h"
using namespace Herwig;
ShowerApproximationKernel::ShowerApproximationKernel()
: thePresampling(false), thePresamplingPoints(10000),
theMaxTry(100000), sampler(0) {}
ShowerApproximationKernel::~ShowerApproximationKernel() {}
IBPtr ShowerApproximationKernel::clone() const {
return new_ptr(*this);
}
IBPtr ShowerApproximationKernel::fullclone() const {
return new_ptr(*this);
}
void ShowerApproximationKernel::showerApproximationGenerator(Ptr<ShowerApproximationGenerator>::tptr gen) {
theShowerApproximationGenerator = gen;
}
Ptr<ShowerApproximationGenerator>::tptr ShowerApproximationKernel::showerApproximationGenerator() const {
return theShowerApproximationGenerator;
}
const vector<bool>& ShowerApproximationKernel::sampleFlags() {
if ( !theFlags.empty() )
return theFlags;
theFlags.resize(nDim(),false);
for ( int k = nDimBorn();
k < nDimBorn() + dipole()->nDimRadiation(); ++k )
theFlags[k] = true;
return theFlags;
}
const pair<vector<double>,vector<double> >& ShowerApproximationKernel::support() {
if ( !theSupport.first.empty() )
return theSupport;
vector<double> l(nDim(),0.0);
vector<double> u(nDim(),1.0);
theSupport.first = l;
theSupport.second = u;
return theSupport;
}
const vector<double>& ShowerApproximationKernel::parameterPoint() {
theLastParameterPoint.resize(nDim());
copy(bornCXComb()->lastRandomNumbers().begin(),
bornCXComb()->lastRandomNumbers().end(),
theLastParameterPoint.begin());
theLastParameterPoint[evolutionVariable()] = 1.;
return theLastParameterPoint;
}
void ShowerApproximationKernel::startPresampling() {
thePresampling = true;
}
void ShowerApproximationKernel::stopPresampling() {
showerApproximationGenerator()->restore();
thePresampling = false;
}
double ShowerApproximationKernel::evaluate(const vector<double>& r) {
if ( presampling() ) {
theLastBornPoint.resize(nDimBorn());
copy(r.begin(),r.begin()+nDimBorn(),theLastBornPoint.begin());
if ( !showerApproximationGenerator()->generate(theLastBornPoint) )
return 0.;
}
assert(dipole()->splitting());
realXComb()->clean();
dipole()->setXComb(realXComb());
for ( vector<StdXCombPtr>::const_iterator t = tildeXCombs().begin();
t != tildeXCombs().end(); ++t ) {
(**t).clean();
(**t).matrixElement()->setXComb(*t);
}
if ( !dipole()->generateKinematics(&r[nDimBorn()]) )
return 0.;
double jac = dipole()->invertedTildeKinematics()->jacobian();
showerApproximation()->setBornXComb(bornXComb());
showerApproximation()->setRealXComb(realXComb());
showerApproximation()->setTildeXCombs(tildeXCombs());
showerApproximation()->setDipole(dipole());
if ( !showerApproximation()->isInShowerPhasespace() )
return 0.;
return showerApproximation()->me2() * jac;
}
double ShowerApproximationKernel::generate() {
if ( !sampler ) {
sampler = new ExponentialGenerator();
sampler->sampling_parameters().maxtry = maxtry();
sampler->sampling_parameters().presampling_points = presamplingPoints();
sampler->function(this);
sampler->initialize();
}
double res = 0.;
while (true) {
try {
res = sampler->generate();
} catch (exsample::exponential_regenerate&) {
continue;
} catch (exsample::hit_and_miss_maxtry&) {
throw Veto();
} catch (exsample::selection_maxtry&) {
throw Veto();
}
break;
}
return res;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void ShowerApproximationKernel::persistentOutput(PersistentOStream & os) const {
- os << theDipole << theShowerApproximation
- << theBornXComb << theRealXComb << theTildeXCombs
- << thePresamplingPoints << theMaxTry << theFlags << theSupport
- << theLastParameterPoint << theShowerApproximationGenerator;
-
+ os << theDipole << theShowerApproximation << theBornXComb
+ << theRealXComb << theTildeXCombs << thePresampling
+ << thePresamplingPoints << theMaxTry << theFlags
+ << theSupport << theShowerApproximationGenerator
+ << theLastParameterPoint << theLastBornPoint
+ << (sampler ? true : false);
+ if ( sampler )
+ sampler->put(os);
}
void ShowerApproximationKernel::persistentInput(PersistentIStream & is, int) {
- is >> theDipole >> theShowerApproximation
- >> theBornXComb >> theRealXComb >> theTildeXCombs
- >> thePresamplingPoints >> theMaxTry >> theFlags >> theSupport
- >> theLastParameterPoint >> theShowerApproximationGenerator;
+ bool haveSampler;
+ is >> theDipole >> theShowerApproximation >> theBornXComb
+ >> theRealXComb >> theTildeXCombs >> thePresampling
+ >> thePresamplingPoints >> theMaxTry >> theFlags
+ >> theSupport >> theShowerApproximationGenerator
+ >> theLastParameterPoint >> theLastBornPoint
+ >> haveSampler;
+ if ( haveSampler ) {
+ sampler = new ExponentialGenerator();
+ sampler->get(is);
+ sampler->function(this);
+ sampler->initialize();
+ }
}
-
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<ShowerApproximationKernel,HandlerBase>
describeHerwigShowerApproximationKernel("Herwig::ShowerApproximationKernel", "HwMatchbox.so");
void ShowerApproximationKernel::Init() {
static ClassDocumentation<ShowerApproximationKernel> documentation
("ShowerApproximationKernel generates emissions according to a "
"shower approximation entering a NLO matching.");
}
diff --git a/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IILightInvertedTildeKinematics.cc
@@ -1,144 +1,148 @@
// -*- C++ -*-
//
// IILightInvertedTildeKinematics.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the IILightInvertedTildeKinematics class.
//
#include "IILightInvertedTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IILightInvertedTildeKinematics::IILightInvertedTildeKinematics() {}
IILightInvertedTildeKinematics::~IILightInvertedTildeKinematics() {}
IBPtr IILightInvertedTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IILightInvertedTildeKinematics::fullclone() const {
return new_ptr(*this);
}
bool IILightInvertedTildeKinematics::doMap(const double * r) {
if ( ptMax() < ptCut() ) {
jacobian(0.0);
return false;
}
Lorentz5Momentum emitter = bornEmitterMomentum();
Lorentz5Momentum spectator = bornSpectatorMomentum();
double mapping = 1.0;
pair<Energy,double> ptz = generatePtZ(mapping,r);
if ( mapping == 0.0 ) {
jacobian(0.0);
return false;
}
Energy pt = ptz.first;
double z = ptz.second;
double ratio = sqr(pt/lastScale());
double x = z*(1.-z) / ( 1. - z + ratio );
double v = ratio * z / ( 1. - z + ratio );
if ( x < emitterX() || x > 1. ||
v > 1. || v > 1.-x ) {
jacobian(0.0);
return false;
}
mapping /= z*(1.-z);
jacobian(mapping*(sqr(lastScale())/sHat())/(16.*sqr(Constants::pi)));
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = v;
double phi = 2.*Constants::pi*r[2];
Lorentz5Momentum kt = getKt(emitter,spectator,pt,phi);
realEmitterMomentum() = (1./x)*emitter;
realEmissionMomentum() = ((1.-x-v)/x)*emitter+v*spectator+kt;
realSpectatorMomentum() = spectator;
K = realEmitterMomentum() + realSpectatorMomentum() - realEmissionMomentum();
K2 = K.m2();
Ktilde = emitter + spectator;
KplusKtilde = K + Ktilde;
KplusKtilde2 = KplusKtilde.m2();
realEmitterMomentum().setMass(ZERO);
realEmitterMomentum().rescaleEnergy();
realEmissionMomentum().setMass(ZERO);
realEmissionMomentum().rescaleEnergy();
realSpectatorMomentum().setMass(ZERO);
realSpectatorMomentum().rescaleEnergy();
return true;
}
Energy IILightInvertedTildeKinematics::lastPt() const {
Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
double x = subtractionParameters()[0];
double v = subtractionParameters()[1];
return scale * sqrt(v*(1.-x-v)/x);
}
Energy IILightInvertedTildeKinematics::ptMax() const {
double tau = emitterX()*spectatorX();
return (1.-tau)*lastScale()/(2.*sqrt(tau));
}
pair<double,double> IILightInvertedTildeKinematics::zBounds(Energy pt) const {
double tau = emitterX()*spectatorX();
double s = sqrt(1.-sqr(pt/ptMax()));
return make_pair(0.5*(1.+tau-(1.-tau)*s),0.5*(1.+tau+(1.-tau)*s));
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
-void IILightInvertedTildeKinematics::persistentOutput(PersistentOStream &) const {
+void IILightInvertedTildeKinematics::persistentOutput(PersistentOStream & os) const {
+ os << ounit(K,GeV) << ounit(K2,GeV2) << ounit(Ktilde,GeV)
+ << ounit(KplusKtilde,GeV) << ounit(KplusKtilde2,GeV2);
}
-void IILightInvertedTildeKinematics::persistentInput(PersistentIStream &, int) {
+void IILightInvertedTildeKinematics::persistentInput(PersistentIStream & is, int) {
+ is >> iunit(K,GeV) >> iunit(K2,GeV2) >> iunit(Ktilde,GeV)
+ >> iunit(KplusKtilde,GeV) >> iunit(KplusKtilde2,GeV2);
}
void IILightInvertedTildeKinematics::Init() {
static ClassDocumentation<IILightInvertedTildeKinematics> documentation
("IILightInvertedTildeKinematics inverts the initial-initial tilde "
"kinematics.");
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<IILightInvertedTildeKinematics,InvertedTildeKinematics>
describeHerwigIILightInvertedTildeKinematics("Herwig::IILightInvertedTildeKinematics", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.cc b/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.cc
--- a/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.cc
+++ b/MatrixElement/Matchbox/Phasespace/IILightTildeKinematics.cc
@@ -1,106 +1,108 @@
// -*- C++ -*-
//
// IILightTildeKinematics.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the IILightTildeKinematics class.
//
#include "IILightTildeKinematics.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
IILightTildeKinematics::IILightTildeKinematics() {}
IILightTildeKinematics::~IILightTildeKinematics() {}
IBPtr IILightTildeKinematics::clone() const {
return new_ptr(*this);
}
IBPtr IILightTildeKinematics::fullclone() const {
return new_ptr(*this);
}
Lorentz5Momentum IILightTildeKinematics::transform(const Lorentz5Momentum& k) const {
LorentzMomentum res =
k - 2.*((k*(K+Ktilde)/(K+Ktilde).m2())*(K+Ktilde)-((k*K)/(K.m2()))*Ktilde);
return res;
}
bool IILightTildeKinematics::doMap() {
Lorentz5Momentum emitter = realEmitterMomentum();
Lorentz5Momentum emission = realEmissionMomentum();
Lorentz5Momentum spectator = realSpectatorMomentum();
double x = (emitter*spectator - emitter*emission - spectator*emission)/(emitter*spectator);
double v = (emitter*emission)/(emitter*spectator);
subtractionParameters().resize(2);
subtractionParameters()[0] = x;
subtractionParameters()[1] = v;
bornEmitterMomentum() = x * emitter;
bornSpectatorMomentum() = spectator;
bornEmitterMomentum().setMass(ZERO);
bornEmitterMomentum().rescaleEnergy();
bornSpectatorMomentum().setMass(ZERO);
bornSpectatorMomentum().rescaleEnergy();
K = emitter + spectator - emission;
Ktilde = x * emitter + spectator;
return true;
}
Energy IILightTildeKinematics::lastPt() const {
Energy scale = sqrt(2.*(bornEmitterMomentum()*bornSpectatorMomentum()));
double x = subtractionParameters()[0];
double v = subtractionParameters()[1];
return scale * sqrt(v*(1.-x-v)/x);
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
-void IILightTildeKinematics::persistentOutput(PersistentOStream &) const {
+void IILightTildeKinematics::persistentOutput(PersistentOStream & os) const {
+ os << ounit(K,GeV) << ounit(Ktilde,GeV);
}
-void IILightTildeKinematics::persistentInput(PersistentIStream &, int) {
+void IILightTildeKinematics::persistentInput(PersistentIStream & is, int) {
+ is >> iunit(K,GeV) >> iunit(Ktilde,GeV);
}
void IILightTildeKinematics::Init() {
static ClassDocumentation<IILightTildeKinematics> documentation
("IILightTildeKinematics implements the 'tilde' kinematics for "
"a initial-initial subtraction dipole.");
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<IILightTildeKinematics,TildeKinematics>
describeHerwigIILightTildeKinematics("Herwig::IILightTildeKinematics", "HwMatchbox.so");
diff --git a/MatrixElement/Matchbox/Phasespace/MatchboxRambo.cc b/MatrixElement/Matchbox/Phasespace/MatchboxRambo.cc
--- a/MatrixElement/Matchbox/Phasespace/MatchboxRambo.cc
+++ b/MatrixElement/Matchbox/Phasespace/MatchboxRambo.cc
@@ -1,268 +1,268 @@
// -*- C++ -*-
//
// MatchboxRambo.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MatchboxRambo class.
//
#include "MatchboxRambo.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "Herwig++/Utilities/GSLBisection.h"
#include "ThePEG/Cuts/Cuts.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
MatchboxRambo::MatchboxRambo()
: needToReshuffle(false), theMakeReferenceSample(false),
referenceSample(0) {}
MatchboxRambo::~MatchboxRambo() {}
IBPtr MatchboxRambo::clone() const {
return new_ptr(*this);
}
IBPtr MatchboxRambo::fullclone() const {
return new_ptr(*this);
}
static double weights[7] = {
-1.,-1.,
0.039788735772973833942,
0.00012598255637968550463,
1.3296564302788840628E-7,
7.0167897579949011130E-11,
2.2217170114046130768E-14
};
void MatchboxRambo::setXComb(tStdXCombPtr xc) {
MatchboxPhasespace::setXComb(xc);
needToReshuffle = false;
if ( xc ) {
for ( cPDVector::const_iterator d = mePartonData().begin();
d != mePartonData().end(); ++d ) {
if ( (**d).mass() != ZERO ) {
needToReshuffle = true;
break;
}
}
}
}
void MatchboxRambo::dumpReference(const vector<Lorentz5Momentum>& momenta, double weight) const {
*referenceSample << lastX1() << " " << lastX2() << " ";
Boost toLab = (lastPartons().first->momentum() +
lastPartons().second->momentum()).boostVector();
for ( vector<Lorentz5Momentum>::const_iterator p = momenta.begin();
p != momenta.end(); ++p ) {
Lorentz5Momentum pl = *p;
if ( toLab.mag2() > Constants::epsilon )
pl.boost(toLab);
*referenceSample
<< (pl.x()/GeV) << " "
<< (pl.y()/GeV) << " "
<< (pl.z()/GeV) << " "
<< (pl.t()/GeV) << " "
<< (pl.mass()/GeV) << " ";
}
double ymax = lastCuts().yHatMax();
double ymin = lastCuts().yHatMin();
double km = log(lastCuts().sHatMax()/lastCuts().sHatMin());
ymax = min(ymax, log(lastCuts().x1Max()*sqrt(lastS()/lastSHat())));
ymin = max(ymin, -log(lastCuts().x2Max()*sqrt(lastS()/lastSHat())));
*referenceSample << weight*km*(ymax-ymin)/(lastX1()*lastX2()) << "\n" << flush;
}
double MatchboxRambo::generateTwoToNKinematics(const double* r,
vector<Lorentz5Momentum>& momenta) {
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = momenta.begin();
for ( ; pd != mePartonData().end(); ++pd, ++p )
p->setMass((**pd).mass());
if ( theMakeReferenceSample ) {
map<cPDVector,ofstream*>::iterator ref =
referenceSamples.find(mePartonData());
if ( ref == referenceSamples.end() ) {
ostringstream refname;
for ( cPDVector::const_iterator p = mePartonData().begin();
p != mePartonData().end(); ++p ) {
refname << (**p).PDGName();
}
refname << ".rambo";
- referenceSamples[mePartonData()] = new ofstream(refname.str().c_str());
+ referenceSamples[mePartonData()] = new ofstream(refname.str().c_str(),std::ios_base::app);
ref = referenceSamples.find(mePartonData());
*(ref->second) << setprecision(26);
}
assert(ref != referenceSamples.end());
referenceSample = ref->second;
}
size_t offset = dynamic_cast<const Tree2toNDiagram&>(*lastXComb().diagrams().front()).nSpace() > 0 ? 2 : 1;
Energy w = sqrt(lastSHat());
size_t count = 0;
Lorentz5Momentum Q;
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin() + offset;
k != momenta.end(); ++k ) {
Energy q = -w*log(r[count]*r[count+1]);
double ct = 2.*r[count+2]-1.;
double st = sqrt(1.-sqr(ct));
double phi = 2.*Constants::pi*r[count+3];
double cphi = cos(phi);
double sphi = sqrt(1.-sqr(cphi));
if ( phi > Constants::pi )
sphi = -sphi;
(*k).setMass(ZERO);
(*k).setT(q);
(*k).setX(q*cphi*st);
(*k).setY(q*sphi*st);
(*k).setZ(q*ct);
count += 4;
Q += *k;
}
Energy M = sqrt(Q.m2());
double x = w/M;
Boost beta = -(Q.vect() * (1./M));
double gamma = Q.t()/M;
double a = 1./(1.+gamma);
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin() + offset;
k != momenta.end(); ++k ) {
Energy q = (*k).t();
Energy bq = beta*(*k).vect();
(*k).setT(x*(gamma*q+bq));
(*k).setVect(x*((*k).vect()+(q+a*bq)*beta));
}
size_t n = momenta.size()-offset;
double weight = weights[n];
if ( !needToReshuffle ) {
if ( !matchConstraints(momenta) )
return 0.;
fillDiagramWeights();
if ( theMakeReferenceSample )
dumpReference(momenta, weight);
return weight;
}
double xi;
ReshuffleEquation solve(w,mePartonData().begin()+offset,mePartonData().end(),
momenta.begin()+2,momenta.end());
GSLBisection solver(1e-10,1e-8,10000);
try {
xi = solver.value(solve,0.0,1.1);
} catch (GSLBisection::GSLerror) {
return 0.;
} catch (GSLBisection::IntervalError) {
return 0.;
}
weight *= pow(xi,3.*(n-1.));
Energy num = ZERO;
Energy den = ZERO;
cPDVector::const_iterator d = mePartonData().begin()+offset;
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin()+offset;
k != momenta.end(); ++k, ++d ) {
num += (*k).vect().mag2()/(*k).t();
Energy q = (*k).t();
(*k).setT(sqrt(sqr((**d).mass())+xi*xi*sqr((*k).t())));
(*k).setVect(xi*(*k).vect());
weight *= q/(*k).t();
den += (*k).vect().mag2()/(*k).t();
(*k).setMass((**d).mass());
}
if ( !matchConstraints(momenta) )
return 0.;
weight *= num/den;
fillDiagramWeights();
if ( theMakeReferenceSample )
dumpReference(momenta, weight);
return weight;
}
Energy MatchboxRambo::ReshuffleEquation::operator() (double xi) const {
cPDVector::const_iterator d = dataBegin;
vector<Lorentz5Momentum>::const_iterator p = momentaBegin;
Energy res = -w;
for ( ; d != dataEnd; ++d, ++p ) {
res += sqrt(sqr((**d).mass()) +
xi*xi*sqr(p->t()));
}
return res;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void MatchboxRambo::persistentOutput(PersistentOStream & os) const {
- os << theMakeReferenceSample;
+ os << needToReshuffle << theMakeReferenceSample;
}
void MatchboxRambo::persistentInput(PersistentIStream & is, int) {
- is >> theMakeReferenceSample;
+ is >> needToReshuffle >> theMakeReferenceSample;
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<MatchboxRambo,MatchboxPhasespace>
describeHerwigMatchboxRambo("Herwig::MatchboxRambo", "HwMatchbox.so");
void MatchboxRambo::Init() {
static ClassDocumentation<MatchboxRambo> documentation
("MatchboxRambo implements RAMBO phase space generation.");
static Switch<MatchboxRambo,bool> interfaceMakeReferenceSample
("MakeReferenceSample",
"Switch on generation of a reference sample of phasespace points.",
&MatchboxRambo::theMakeReferenceSample, false, false, false);
static SwitchOption interfaceMakeReferenceSampleOn
(interfaceMakeReferenceSample,
"On",
"Generate a reference sample.",
true);
static SwitchOption interfaceMakeReferenceSampleOff
(interfaceMakeReferenceSample,
"Off",
"Do not generate a reference sample.",
false);
}
diff --git a/MatrixElement/Matchbox/Phasespace/TreePhasespace.cc b/MatrixElement/Matchbox/Phasespace/TreePhasespace.cc
--- a/MatrixElement/Matchbox/Phasespace/TreePhasespace.cc
+++ b/MatrixElement/Matchbox/Phasespace/TreePhasespace.cc
@@ -1,259 +1,266 @@
// -*- C++ -*-
//
// TreePhasespace.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the TreePhasespace class.
//
#include "TreePhasespace.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
using namespace Herwig::PhasespaceHelpers;
TreePhasespace::TreePhasespace()
: x0(0.01), xc(1e-4), M0(ZERO), Mc(ZERO) {
lastPhasespaceInfo.x0 = x0;
lastPhasespaceInfo.xc = xc;
lastPhasespaceInfo.M0 = M0;
lastPhasespaceInfo.Mc = Mc;
theIncludeMirrored = true;
}
TreePhasespace::~TreePhasespace() {}
IBPtr TreePhasespace::clone() const {
return new_ptr(*this);
}
IBPtr TreePhasespace::fullclone() const {
return new_ptr(*this);
}
void TreePhasespace::setXComb(tStdXCombPtr xco) {
MatchboxPhasespace::setXComb(xco);
lastChannelsIterator = channelMap().find(lastXCombPtr());
if ( lastChannelsIterator == channelMap().end() ) {
map<Ptr<Tree2toNDiagram>::ptr,pair<PhasespaceTree, PhasespaceTree> > channels;
for ( StandardXComb::DiagramVector::const_iterator d =
lastXComb().diagrams().begin(); d != lastXComb().diagrams().end(); ++d ) {
PhasespaceTree tree;
Ptr<Tree2toNDiagram>::ptr diag =
dynamic_ptr_cast<Ptr<Tree2toNDiagram>::ptr>(*d);
tree.setup(*diag);
PhasespaceTree treeMirror;
treeMirror.setupMirrored(*diag, diag->nSpace() - 1);
channels[diag] = make_pair(tree,treeMirror);
}
channelMap()[lastXCombPtr()] = channels;
lastChannelsIterator = channelMap().find(lastXCombPtr());
}
lastPhasespaceInfo.sHat = lastXComb().lastSHat();
lastPhasespaceInfo.sqrtSHat = sqrt(lastXComb().lastSHat());
lastPhasespaceInfo.weight = 1.;
}
double TreePhasespace::generateTwoToNKinematics(const double* random,
vector<Lorentz5Momentum>& momenta) {
cPDVector::const_iterator pd = mePartonData().begin();
vector<Lorentz5Momentum>::iterator p = momenta.begin();
for ( ; pd != mePartonData().end(); ++pd, ++p )
p->setMass((**pd).mass());
if ( momenta.size() > 3 ) {
size_t nchannels = lastXComb().diagrams().size();
bool doMirror = (UseRandom::rnd() < 0.5) && theIncludeMirrored;
map<Ptr<Tree2toNDiagram>::ptr,
pair <PhasespaceHelpers::PhasespaceTree, PhasespaceHelpers::PhasespaceTree> >::iterator ds =
lastChannels().begin();
advance(ds,(size_t)(random[0]*nchannels));
Ptr<Tree2toNDiagram>::ptr channel = ds->first;
++random;
lastPhasespaceInfo.rnd.numbers = random;
lastPhasespaceInfo.rnd.nRnd = 3*momenta.size() - 10;
try {
if ( !doMirror )
lastChannels()[channel].first.generateKinematics(lastPhasespaceInfo,momenta);
else
lastChannels()[channel].second.generateKinematics(lastPhasespaceInfo,momenta);
} catch (Veto) {
return 0.;
}
if ( !matchConstraints(momenta) )
return 0.;
double flatCut = x0;
if ( M0 != ZERO )
flatCut = M0/sqrt(lastSHat());
fillDiagramWeights(flatCut);
double sum = 0.;
for ( map<Ptr<Tree2toNDiagram>::ptr,
pair <PhasespaceHelpers::PhasespaceTree, PhasespaceHelpers::PhasespaceTree> >::const_iterator d
= lastChannels().begin(); d != lastChannels().end(); ++d )
sum += diagramWeight(*(d->first));
double piWeight = pow(2.*Constants::pi,(double)(3*(momenta.size()-2)-4));
for ( vector<Lorentz5Momentum>::iterator k = momenta.begin();
k != momenta.end(); ++k )
k->rescaleRho();
return nchannels*lastPhasespaceInfo.weight*diagramWeight(*channel)/(sum*piWeight);
}
double tau = momenta[2].mass2()/lastXCombPtr()->lastS();
double ltau = log(tau)/2.;
double y = ltau - 2.*random[0]*ltau;
double x1 = sqrt(tau)*exp(y);
double x2 = sqrt(tau)*exp(-y);
ThreeVector<Energy> p1 =
x1*(lastXCombPtr()->lastParticles().first->momentum().vect());
ThreeVector<Energy> p2 =
x2*(lastXCombPtr()->lastParticles().second->momentum().vect());
ThreeVector<Energy> q = p1 + p2;
momenta[0] = Lorentz5Momentum(momenta[0].mass(),p1);
momenta[1] = Lorentz5Momentum(momenta[1].mass(),p2);
momenta[2] = Lorentz5Momentum(momenta[2].mass(),q);
lastXCombPtr()->lastX1X2(make_pair(x1,x2));
lastXCombPtr()->lastSHat((momenta[0]+momenta[1]).m2());
fillDiagramWeights();
double weight = -4.*Constants::pi*ltau;
return weight;
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
void TreePhasespace::doinit() {
MatchboxPhasespace::doinit();
lastPhasespaceInfo.x0 = x0;
lastPhasespaceInfo.xc = xc;
lastPhasespaceInfo.M0 = M0;
lastPhasespaceInfo.Mc = Mc;
}
void TreePhasespace::doinitrun() {
MatchboxPhasespace::doinitrun();
lastPhasespaceInfo.x0 = x0;
lastPhasespaceInfo.xc = xc;
lastPhasespaceInfo.M0 = M0;
lastPhasespaceInfo.Mc = Mc;
}
void TreePhasespace::persistentOutput(PersistentOStream & os) const {
os << theChannelMap << x0 << xc
<< ounit(M0,GeV) << ounit(Mc,GeV)
- << theIncludeMirrored;
+ << theIncludeMirrored
+ << theLastXComb;
}
void TreePhasespace::persistentInput(PersistentIStream & is, int) {
is >> theChannelMap >> x0 >> xc
>> iunit(M0,GeV) >> iunit(Mc,GeV)
- >> theIncludeMirrored;
+ >> theIncludeMirrored
+ >> theLastXComb;
+ lastPhasespaceInfo.x0 = x0;
+ lastPhasespaceInfo.xc = xc;
+ lastPhasespaceInfo.M0 = M0;
+ lastPhasespaceInfo.Mc = Mc;
+ lastChannelsIterator = channelMap().find(lastXCombPtr());
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeClass<TreePhasespace,MatchboxPhasespace>
describeHerwigTreePhasespace("Herwig::TreePhasespace", "HwMatchbox.so");
void TreePhasespace::Init() {
static ClassDocumentation<TreePhasespace> documentation
("TreePhasespace is a multichannel phasespace generator "
"adapting to singularity structures as determined from the matrix "
"elements diagrams.");
static Reference<TreePhasespace,TreePhasespaceChannels> interfaceChannelMap
("ChannelMap",
"Set the object storing the channels.",
&TreePhasespace::theChannelMap, false, false, true, false, false);
static Parameter<TreePhasespace,double> interfaceX0
("X0",
"Set the cut below which flat virtuality sampling is imposed.",
&TreePhasespace::x0, 0.01, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<TreePhasespace,double> interfaceXC
("XC",
"Set the cut below which no virtualities are generated.",
&TreePhasespace::xc, 1e-4, 0.0, 0,
false, false, Interface::lowerlim);
static Parameter<TreePhasespace,Energy> interfaceM0
("M0",
"Set the cut below which flat virtuality sammpling is imposed.",
&TreePhasespace::M0, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Parameter<TreePhasespace,Energy> interfaceMC
("MC",
"Set the cut below which no virtualities are generated.",
&TreePhasespace::Mc, GeV, 0.0*GeV, 0.0*GeV, 0*GeV,
false, false, Interface::lowerlim);
static Switch<TreePhasespace,bool> interfaceIncludeMirrored
("IncludeMirrored",
"Choose whether to include mirrored diagrams for PS generation",
&TreePhasespace::theIncludeMirrored, true, true, false);
static SwitchOption interfaceIncludeMirroredTrue
(interfaceIncludeMirrored,
"True",
"Use unmirrored and mirrored diagrams",
true);
static SwitchOption interfaceIncludeMirroredFalse
(interfaceIncludeMirrored,
"False",
"Use only unmirrored diagrams",
false);
}
diff --git a/MatrixElement/Matchbox/Utility/ColourBasis.cc b/MatrixElement/Matchbox/Utility/ColourBasis.cc
--- a/MatrixElement/Matchbox/Utility/ColourBasis.cc
+++ b/MatrixElement/Matchbox/Utility/ColourBasis.cc
@@ -1,1189 +1,1190 @@
// -*- C++ -*-
//
// ColourBasis.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ColourBasis class.
//
#include "ColourBasis.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <iterator>
using std::ostream_iterator;
#include "DiagramDrawer.h"
using namespace Herwig;
using boost::numeric::ublas::trans;
using boost::numeric::ublas::conj;
using boost::numeric::ublas::row;
using boost::numeric::ublas::column;
using boost::numeric::ublas::prod;
ColourBasis::ColourBasis()
: theLargeN(false), theSearchPath("."), didRead(false), didWrite(false) {}
ColourBasis::~ColourBasis() {
for ( map<Ptr<Tree2toNDiagram>::tcptr,vector<ColourLines*> >::iterator cl =
theColourLineMap.begin(); cl != theColourLineMap.end(); ++cl ) {
for ( vector<ColourLines*>::iterator c = cl->second.begin();
c != cl->second.end(); ++c ) {
if ( *c )
delete *c;
}
}
theColourLineMap.clear();
}
// If needed, insert default implementations of virtual function defined
// in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs).
bool ColourBasis::colourConnected(const cPDVector& sub,
const vector<PDT::Colour>& basis,
const pair<int,bool>& i,
const pair<int,bool>& j,
size_t a) const {
// translate process to basis ids
map<cPDVector,map<size_t,size_t> >::const_iterator trans
= indexMap().find(sub);
assert(trans != indexMap().end());
int idColoured = i.second ? j.first : i.first;
idColoured = trans->second.find(idColoured)->second;
int idAntiColoured = i.second ? i.first : j.first;
idAntiColoured = trans->second.find(idAntiColoured)->second;
return colourConnected(basis,idColoured,idAntiColoured,a);
}
const string& ColourBasis::orderingString(const cPDVector& sub,
const map<size_t,size_t>& colourToAmplitude,
size_t tensorId) {
const vector<PDT::Colour>& basis = normalOrderedLegs(sub);
map<size_t,string>& orderings = theOrderingStringIdentifiers[basis][colourToAmplitude];
if ( orderings.empty() ) {
map<size_t,vector<vector<size_t> > > tensors =
basisList(basis);
for ( map<size_t,vector<vector<size_t> > >::const_iterator t =
tensors.begin(); t != tensors.end(); ++t ) {
ostringstream oid;
for ( vector<vector<size_t> >::const_iterator s = t->second.begin();
s != t->second.end(); ++s ) {
oid << "[";
for ( vector<size_t>::const_iterator l = s->begin();
l != s->end(); ++l ) {
map<size_t,size_t>::const_iterator trans =
colourToAmplitude.find(*l);
assert(trans != colourToAmplitude.end());
oid << trans->second << (l != --(s->end()) ? "," : "");
}
oid << "]";
}
orderings[t->first] = oid.str();
}
}
return orderings[tensorId];
}
const vector<vector<size_t> >& ColourBasis::ordering(const cPDVector& sub,
const map<size_t,size_t>& colourToAmplitude,
size_t tensorId) {
const vector<PDT::Colour>& basis = normalOrderedLegs(sub);
map<size_t,vector<vector<size_t> > >& orderings = theOrderingIdentifiers[basis][colourToAmplitude];
if ( orderings.empty() ) {
map<size_t,vector<vector<size_t> > > tensors =
basisList(basis);
for ( map<size_t,vector<vector<size_t> > >::const_iterator t =
tensors.begin(); t != tensors.end(); ++t ) {
orderings[t->first] = t->second;
}
}
return orderings[tensorId];
}
vector<PDT::Colour> ColourBasis::normalOrderMap(const cPDVector& sub) {
vector<PDT::Colour> allLegs = projectColour(sub);
vector<PDT::Colour> legs = normalOrder(allLegs);
if ( allLegs[0] == PDT::Colour3 )
allLegs[0] = PDT::Colour3bar;
else if ( allLegs[0] == PDT::Colour3bar )
allLegs[0] = PDT::Colour3;
if ( allLegs[1] == PDT::Colour3 )
allLegs[1] = PDT::Colour3bar;
else if ( allLegs[1] == PDT::Colour3bar )
allLegs[1] = PDT::Colour3;
if ( theIndexMap.find(sub) == theIndexMap.end() ) {
map<size_t,size_t> trans;
vector<PDT::Colour> checkLegs = legs;
size_t n = checkLegs.size();
for ( size_t i = 0; i < allLegs.size(); ++i ) {
size_t j = 0;
while ( checkLegs[j] != allLegs[i] ) {
++j; if ( j == n ) break;
}
if ( j == n ) continue;
trans[i] = j;
checkLegs[j] = PDT::ColourUndefined;
}
theIndexMap[sub] = trans;
}
return legs;
}
const vector<PDT::Colour>& ColourBasis::normalOrderedLegs(const cPDVector& sub) const {
static vector<PDT::Colour> empty;
map<cPDVector,vector<PDT::Colour> >::const_iterator n =
theNormalOrderedLegs.find(sub);
if ( n != theNormalOrderedLegs.end() )
return n->second;
return empty;
}
size_t ColourBasis::prepare(const cPDVector& sub,
bool noCorrelations) {
vector<PDT::Colour> legs = normalOrderMap(sub);
bool doPrepare = false;
if ( theNormalOrderedLegs.find(sub) == theNormalOrderedLegs.end() )
theNormalOrderedLegs[sub] = legs;
if ( theScalarProducts.find(legs) == theScalarProducts.end() )
doPrepare = true;
if ( doPrepare )
doPrepare = !readBasis(legs);
size_t dim = doPrepare ? prepareBasis(legs) : theScalarProducts[legs].size1();
if ( theCharges.find(legs) != theCharges.end() )
return dim;
if ( !doPrepare && noCorrelations )
return dim;
symmetric_matrix<double,upper>& sp =
theScalarProducts.insert(make_pair(legs,symmetric_matrix<double,upper>(dim,dim))).first->second;
for ( size_t a = 0; a < dim; ++a )
for ( size_t b = a; b < dim; ++b )
sp(a,b) = scalarProduct(a,b,legs);
if ( noCorrelations )
return dim;
vector<PDT::Colour> legsPlus = legs;
legsPlus.push_back(PDT::Colour8);
legsPlus = normalOrder(legsPlus);
bool doPreparePlus = theScalarProducts.find(legsPlus) == theScalarProducts.end();
size_t dimPlus = doPreparePlus ? prepareBasis(legsPlus) : theScalarProducts[legsPlus].size1();
symmetric_matrix<double,upper>& spPlus =
doPreparePlus ?
theScalarProducts.insert(make_pair(legsPlus,symmetric_matrix<double,upper>(dimPlus,dimPlus))).first->second :
theScalarProducts[legsPlus];
if ( doPreparePlus ) {
for ( size_t a = 0; a < dimPlus; ++a )
for ( size_t b = a; b < dimPlus; ++b )
spPlus(a,b) = scalarProduct(a,b,legsPlus);
}
typedef map<size_t,compressed_matrix<double> > cMap;
cMap& cm = theCharges.insert(make_pair(legs,cMap())).first->second;
typedef map<size_t,vector<pair<size_t,size_t> > > ccMap;
ccMap& ccm = theChargeNonZeros.insert(make_pair(legs,ccMap())).first->second;
tmp.resize(dimPlus,dim);
for ( size_t i = 0; i < legs.size(); ++i ) {
size_t nonZero = 0;
vector<pair<size_t,size_t> > nonZeros;
for ( size_t a = 0; a < dimPlus; ++a )
for ( size_t b = 0; b < dim; ++b ) {
tmp(a,b) = tMatrixElement(i,a,b,legsPlus,legs);
if ( tmp(a,b) != 0. ) {
++nonZero;
nonZeros.push_back(make_pair(a,b));
}
}
ccm.insert(make_pair(i,nonZeros));
compressed_matrix<double>& tm =
cm.insert(make_pair(i,compressed_matrix<double>(dimPlus,dim,nonZero))).first->second;
for ( size_t a = 0; a < dimPlus; ++a )
for ( size_t b = 0; b < dim; ++b ) {
if ( tmp(a,b) != 0. )
tm(a,b) = tmp(a,b);
}
}
map<pair<size_t,size_t>,symmetric_matrix<double,upper> >& xm = theCorrelators[legs];
for ( size_t i = 0; i < legs.size(); ++i )
for ( size_t j = i+1; j < legs.size(); ++j ) {
symmetric_matrix<double,upper>& mm =
xm.insert(make_pair(make_pair(i,j),symmetric_matrix<double,upper>(dim,dim))).first->second;
chargeProduct(cm[i],ccm[i],spPlus,cm[j],ccm[j],mm);
}
return dim;
}
void ColourBasis::chargeProduct(const compressed_matrix<double>& ti,
const vector<pair<size_t,size_t> >& tiNonZero,
const symmetric_matrix<double,upper>& X,
const compressed_matrix<double>& tj,
const vector<pair<size_t,size_t> >& tjNonZero,
symmetric_matrix<double,upper>& result) const {
for ( size_t i = 0; i < result.size1(); ++i )
for ( size_t j = i; j < result.size1(); ++j )
result(i,j) = 0.;
for ( vector<pair<size_t,size_t> >::const_iterator i = tiNonZero.begin();
i != tiNonZero.end(); ++i )
for ( vector<pair<size_t,size_t> >::const_iterator j = tjNonZero.begin();
j != tjNonZero.end(); ++j ) {
if ( j->second < i->second )
continue;
result(i->second,j->second) +=
ti(i->first,i->second)*tj(j->first,j->second)*X(i->first,j->first);
}
}
void ColourBasis::chargeProductAdd(const compressed_matrix<double>& ti,
const vector<pair<size_t,size_t> >& tiNonZero,
const matrix<Complex>& X,
const compressed_matrix<double>& tj,
const vector<pair<size_t,size_t> >& tjNonZero,
matrix<Complex>& result,
double factor) const {
for ( vector<pair<size_t,size_t> >::const_iterator i = tiNonZero.begin();
i != tiNonZero.end(); ++i )
for ( vector<pair<size_t,size_t> >::const_iterator j = tjNonZero.begin();
j != tjNonZero.end(); ++j ) {
result(i->first,j->first) += factor*
ti(i->first,i->second)*tj(j->first,j->second)*X(i->second,j->second);
}
}
string ColourBasis::cfstring(const list<list<pair<int,bool> > >& flow) {
ostringstream out("");
for ( list<list<pair<int,bool> > >::const_iterator line =
flow.begin(); line != flow.end(); ++line ) {
for ( list<pair<int,bool> >::const_iterator node =
line->begin(); node != line->end(); ++node ) {
out << (node->second ? "-" : "") << (node->first+1) << " ";
}
if ( line != --(flow.end()) )
out << ", ";
}
return out.str();
}
vector<string> ColourBasis::makeFlows(Ptr<Tree2toNDiagram>::tcptr diag,
size_t dim) const {
vector<string> res(dim);
list<list<list<pair<int,bool> > > > fdata =
colourFlows(diag);
cPDVector ext;
tcPDVector dext = diag->external();
copy(dext.begin(),dext.end(),back_inserter(ext));
vector<PDT::Colour> colouredLegs =
normalOrder(projectColour(ext));
for ( list<list<list<pair<int,bool> > > >::const_iterator flow =
fdata.begin(); flow != fdata.end(); ++flow ) {
for ( size_t i = 0; i < dim; ++i ) {
bool matches = true;
for ( list<list<pair<int,bool> > >::const_iterator line =
flow->begin(); line != flow->end(); ++line ) {
pair<int,bool> front(diag->externalId(line->front().first),line->front().second);
if ( front.first < 2 )
front.second = !front.second;
pair<int,bool> back(diag->externalId(line->back().first),line->back().second);
if ( back.first < 2 )
back.second = !back.second;
if ( !colourConnected(ext,colouredLegs,front,back,i) ) {
matches = false;
break;
}
}
if ( matches ) {
res[i] = cfstring(*flow);
}
}
}
bool gotone = false;
for ( vector<string>::const_iterator f = res.begin();
f != res.end(); ++f ) {
if ( *f != "" ) {
gotone = true;
break;
}
}
if ( !gotone ) {
generator()->log() << "warning no color flow found for diagram\n";
DiagramDrawer::drawDiag(generator()->log(),*diag);
}
return res;
}
size_t ColourBasis::prepare(const MEBase::DiagramVector& diags,
bool noCorrelations) {
size_t dim = 0;
for ( MEBase::DiagramVector::const_iterator d = diags.begin();
d != diags.end(); ++d ) {
Ptr<Tree2toNDiagram>::tcptr dd = dynamic_ptr_cast<Ptr<Tree2toNDiagram>::ptr>(*d);
assert(dd);
dim = prepare(dd->partons(),noCorrelations);
if ( !haveColourFlows() || theFlowMap.find(dd) != theFlowMap.end() )
continue;
theFlowMap[dd] = makeFlows(dd,dim);
}
return dim;
}
bool matchEnd(int a, pair<int,bool> b,
Ptr<Tree2toNDiagram>::tcptr diag) {
if ( a != b.first )
return false;
if ( b.first != diag->nSpace()-1 ) {
return
!b.second ?
diag->allPartons()[b.first]->hasColour() :
diag->allPartons()[b.first]->hasAntiColour();
} else {
return
!b.second ?
diag->allPartons()[b.first]->hasAntiColour() :
diag->allPartons()[b.first]->hasColour();
}
return false;
}
bool findPath(pair<int,bool> a, pair<int,bool> b,
Ptr<Tree2toNDiagram>::tcptr diag,
list<pair<int,bool> >& path,
bool backward) {
assert(a.first==0 ? !backward : true);
if ( path.empty() )
path.push_back(a);
if ( !backward ) {
if ( diag->children(a.first).first == -1 )
return matchEnd(a.first,b,diag);
pair<int,int> children = diag->children(a.first);
bool cc = (children.first == diag->nSpace()-1);
if ( diag->allPartons()[children.first]->coloured() )
if ( !cc ?
(!a.second ?
diag->allPartons()[children.first]->hasColour() :
diag->allPartons()[children.first]->hasAntiColour()) :
(!a.second ?
diag->allPartons()[children.first]->hasAntiColour() :
diag->allPartons()[children.first]->hasColour()) ) {
pair<int,bool> next(children.first,a.second);
path.push_back(next);
if ( !findPath(next,b,diag,path,false) ) {
path.pop_back();
} else return true;
}
cc = (children.second == diag->nSpace()-1);
if ( diag->allPartons()[children.second]->coloured() )
if ( !cc ?
(!a.second ?
diag->allPartons()[children.second]->hasColour() :
diag->allPartons()[children.second]->hasAntiColour()) :
(!a.second ?
diag->allPartons()[children.second]->hasAntiColour() :
diag->allPartons()[children.second]->hasColour()) ) {
pair<int,bool> next(children.second,a.second);
path.push_back(next);
if ( !findPath(next,b,diag,path,false) ) {
path.pop_back();
} else return true;
}
if ( path.size() == 1 )
path.pop_back();
return false;
} else {
int parent = diag->parent(a.first);
pair<int,int> neighbours = diag->children(parent);
int neighbour = a.first == neighbours.first ? neighbours.second : neighbours.first;
if ( matchEnd(parent,b,diag) ) {
path.push_back(b);
return true;
}
if ( matchEnd(neighbour,b,diag) ) {
path.push_back(b);
return true;
}
if ( diag->allPartons()[neighbour]->coloured() )
if ( a.second ?
diag->allPartons()[neighbour]->hasColour() :
diag->allPartons()[neighbour]->hasAntiColour() ) {
pair<int,bool> next(neighbour,!a.second);
path.push_back(next);
if ( !findPath(next,b,diag,path,false) ) {
path.pop_back();
} else return true;
}
if ( parent == 0 ) {
if ( path.size() == 1 )
path.pop_back();
return false;
}
if ( diag->allPartons()[parent]->coloured() )
if ( !a.second ?
diag->allPartons()[parent]->hasColour() :
diag->allPartons()[parent]->hasAntiColour() ) {
pair<int,bool> next(parent,a.second);
path.push_back(next);
if ( !findPath(next,b,diag,path,true) ) {
path.pop_back();
} else return true;
}
if ( path.size() == 1 )
path.pop_back();
return false;
}
return false;
}
list<pair<int,bool> > ColourBasis::colouredPath(pair<int,bool> a, pair<int,bool> b,
Ptr<Tree2toNDiagram>::tcptr diag) {
list<pair<int,bool> > res;
if ( a.first == b.first )
return res;
bool aIn = (a.first < 2);
bool bIn = (b.first < 2);
if ( (aIn && bIn) || (!aIn && !bIn) )
if ( (a.second && b.second) ||
(!a.second && !b.second) )
return res;
if ( (aIn && !bIn) || (!aIn && bIn) )
if ( (!a.second && b.second) ||
(a.second && !b.second) )
return res;
if ( a.first > b.first )
swap(a,b);
a.first = diag->diagramId(a.first);
b.first = diag->diagramId(b.first);
if ( a.first == diag->nSpace()-1 )
a.second = !a.second;
if ( b.first == diag->nSpace()-1 )
b.second = !b.second;
if ( !findPath(a,b,diag,res,a.first != 0) )
return res;
if ( b.first == diag->nSpace()-1 ) {
res.back().second = !res.back().second;
}
if ( a.first == diag->nSpace()-1 ) {
res.front().second = !res.front().second;
}
return res;
}
list<list<list<pair<int,bool> > > >
ColourBasis::colourFlows(Ptr<Tree2toNDiagram>::tcptr diag) {
vector<pair<int,bool> > connectSource;
vector<pair<int,bool> > connectSink;
for ( size_t i = 0; i != diag->partons().size(); ++i ) {
if ( i < 2 && diag->partons()[i]->hasAntiColour() )
connectSource.push_back(make_pair(i,true));
if ( i < 2 && diag->partons()[i]->hasColour() )
connectSink.push_back(make_pair(i,false));
if ( i > 1 && diag->partons()[i]->hasColour() )
connectSource.push_back(make_pair(i,false));
if ( i > 1 && diag->partons()[i]->hasAntiColour() )
connectSink.push_back(make_pair(i,true));
}
assert(connectSource.size() == connectSink.size());
list<list<list<pair<int,bool> > > > ret;
do {
vector<pair<int,bool> >::iterator source =
connectSource.begin();
vector<pair<int,bool> >::iterator sink =
connectSink.begin();
list<list<pair<int,bool> > > res;
for ( ; source != connectSource.end(); ++source, ++sink ) {
if ( source->first == sink->first ) {
res.clear();
break;
}
list<pair<int,bool> > line =
colouredPath(*source,*sink,diag);
if ( line.empty() ) {
res.clear();
break;
}
res.push_back(line);
}
if ( !res.empty() ) {
// check, if all dressed properly
vector<pair<int,int> > dressed((*diag).allPartons().size(),make_pair(0,0));
for ( size_t p = 0; p < diag->allPartons().size(); ++p ) {
if ( diag->allPartons()[p]->hasColour() &&
!diag->allPartons()[p]->hasAntiColour() )
dressed[p].first = 1;
if ( diag->allPartons()[p]->hasAntiColour() &&
!diag->allPartons()[p]->hasColour() )
dressed[p].second = 1;
if ( diag->allPartons()[p]->hasAntiColour() &&
diag->allPartons()[p]->hasColour() ) {
dressed[p].first = 1; dressed[p].second = 1;
}
}
for ( list<list<pair<int,bool> > >::const_iterator l = res.begin();
l != res.end(); ++l ) {
for ( list<pair<int,bool> >::const_iterator n = l->begin();
n != l->end(); ++n ) {
if ( !(n->second) )
dressed[n->first].first -= 1;
else
dressed[n->first].second -= 1;
}
}
for ( vector<pair<int,int> >::const_iterator d = dressed.begin();
d != dressed.end(); ++d ) {
if ( d->first != 0 || d->second != 0 ) {
res.clear();
break;
}
}
if ( !res.empty() )
ret.push_back(res);
}
} while ( std::next_permutation(connectSink.begin(),connectSink.end()) );
return ret;
}
map<Ptr<Tree2toNDiagram>::tcptr,vector<ColourLines*> >&
ColourBasis::colourLineMap() {
if ( !theColourLineMap.empty() )
return theColourLineMap;
for ( map<Ptr<Tree2toNDiagram>::tcptr,vector<string> >::const_iterator cl =
theFlowMap.begin(); cl != theFlowMap.end(); ++cl ) {
vector<ColourLines*> clines(cl->second.size());
for ( size_t k = 0; k < cl->second.size(); ++k ) {
if ( cl->second[k] == "" ) {
clines[k] = 0;
continue;
}
clines[k] = new ColourLines(cl->second[k]);
}
theColourLineMap[cl->first] = clines;
}
return theColourLineMap;
}
Selector<const ColourLines *> ColourBasis::colourGeometries(tcDiagPtr diag,
const map<vector<int>,CVector>& amps) {
Ptr<Tree2toNDiagram>::tcptr dd =
dynamic_ptr_cast<Ptr<Tree2toNDiagram>::tcptr>(diag);
assert(dd && theFlowMap.find(dd) != theFlowMap.end());
const vector<ColourLines*>& cl = colourLineMap()[dd];
Selector<const ColourLines *> sel;
size_t dim = amps.begin()->second.size();
assert(dim == cl.size());
double w = 0.;
for ( size_t i = 0; i < dim; ++i ) {
if ( !cl[i] )
continue;
w = 0.;
for ( map<vector<int>,CVector>::const_iterator a = amps.begin();
a != amps.end(); ++a )
w += real(conj((a->second)(i))*((a->second)(i)));
if ( w > 0. )
sel.insert(w,cl[i]);
}
assert(!sel.empty());
return sel;
}
const symmetric_matrix<double,upper>& ColourBasis::scalarProducts(const cPDVector& sub) const {
map<cPDVector,vector<PDT::Colour> >::const_iterator lit =
theNormalOrderedLegs.find(sub);
assert(lit != theNormalOrderedLegs.end());
ScalarProductMap::const_iterator spit =
theScalarProducts.find(lit->second);
assert(spit != theScalarProducts.end());
return spit->second;
}
const compressed_matrix<double>& ColourBasis::charge(const cPDVector& sub, size_t iIn) const {
map<cPDVector,vector<PDT::Colour> >::const_iterator lit =
theNormalOrderedLegs.find(sub);
assert(lit != theNormalOrderedLegs.end());
ChargeMap::const_iterator ct =
theCharges.find(lit->second);
assert(ct != theCharges.end());
map<cPDVector,map<size_t,size_t> >::const_iterator trans
= theIndexMap.find(sub);
assert(trans != theIndexMap.end());
size_t i = trans->second.find(iIn)->second;
map<size_t,compressed_matrix<double> >::const_iterator cit
= ct->second.find(i);
assert(cit != ct->second.end());
return cit->second;
}
const vector<pair<size_t,size_t> >& ColourBasis::chargeNonZero(const cPDVector& sub, size_t iIn) const {
map<cPDVector,vector<PDT::Colour> >::const_iterator lit =
theNormalOrderedLegs.find(sub);
assert(lit != theNormalOrderedLegs.end());
ChargeNonZeroMap::const_iterator ct =
theChargeNonZeros.find(lit->second);
assert(ct != theChargeNonZeros.end());
map<cPDVector,map<size_t,size_t> >::const_iterator trans
= theIndexMap.find(sub);
assert(trans != theIndexMap.end());
size_t i = trans->second.find(iIn)->second;
map<size_t,vector<pair<size_t,size_t> > >::const_iterator cit
= ct->second.find(i);
assert(cit != ct->second.end());
return cit->second;
}
const symmetric_matrix<double,upper>& ColourBasis::correlator(const cPDVector& sub,
const pair<size_t,size_t>& ijIn) const {
map<cPDVector,vector<PDT::Colour> >::const_iterator lit =
theNormalOrderedLegs.find(sub);
assert(lit != theNormalOrderedLegs.end());
CorrelatorMap::const_iterator cit =
theCorrelators.find(lit->second);
assert(cit != theCorrelators.end());
map<cPDVector,map<size_t,size_t> >::const_iterator trans
= theIndexMap.find(sub);
assert(trans != theIndexMap.end());
pair<size_t,size_t> ij(trans->second.find(ijIn.first)->second,
trans->second.find(ijIn.second)->second);
if ( ij.first > ij.second )
swap(ij.first,ij.second);
map<pair<size_t,size_t>,symmetric_matrix<double,upper> >::const_iterator cijit
= cit->second.find(ij);
assert(cijit != cit->second.end());
return cijit->second;
}
double ColourBasis::me2(const cPDVector& sub,
const map<vector<int>,CVector>& amps) const {
const symmetric_matrix<double,upper>& sp = scalarProducts(sub);
double res = 0.;
for ( map<vector<int>,CVector>::const_iterator a = amps.begin();
a != amps.end(); ++a ) {
res += real(inner_prod(conj(a->second),prod(sp,a->second)));
}
return res;
}
double ColourBasis::interference(const cPDVector& sub,
const map<vector<int>,CVector>& amps1,
const map<vector<int>,CVector>& amps2) const {
const symmetric_matrix<double,upper>& sp = scalarProducts(sub);
double res = 0.;
map<vector<int>,CVector>::const_iterator a = amps1.begin();
map<vector<int>,CVector>::const_iterator b = amps2.begin();
for ( ; a != amps1.end(); ++a, ++b ) {
assert(a->first == b->first);
res += 2.*real(inner_prod(conj(a->second),prod(sp,b->second)));
}
assert(!isnan(res));
return res;
}
double ColourBasis::colourCorrelatedME2(const pair<size_t,size_t>& ij,
const cPDVector& sub,
const map<vector<int>,CVector>& amps) const {
const symmetric_matrix<double,upper>& cij = correlator(sub,ij);
double res = 0.;
for ( map<vector<int>,CVector>::const_iterator a = amps.begin();
a != amps.end(); ++a ) {
res += real(inner_prod(conj(a->second),prod(cij,a->second)));
}
return res;
}
Complex ColourBasis::interference(const cPDVector& sub,
const CVector& left,
const CVector& right) const {
const symmetric_matrix<double,upper>& sp = scalarProducts(sub);
return inner_prod(conj(left),prod(sp,right));
}
Complex ColourBasis::colourCorrelatedInterference(const pair<size_t,size_t>& ij,
const cPDVector& sub,
const CVector& left,
const CVector& right) const {
const symmetric_matrix<double,upper>& cij = correlator(sub,ij);
return inner_prod(conj(left),prod(cij,right));
}
double ColourBasis::me2(const cPDVector& sub,
const matrix<Complex>& amp) const {
const symmetric_matrix<double,upper>& sp = scalarProducts(sub);
double tr = 0;
size_t n = amp.size1();
for ( size_t i = 0; i < n; ++i ) {
tr += real(inner_prod(row(sp,i),column(amp,i)));
}
return tr;
}
double ColourBasis::colourCorrelatedME2(const pair<size_t,size_t>& ij,
const cPDVector& sub,
const matrix<Complex>& amp) const {
const symmetric_matrix<double,upper>& cij = correlator(sub,ij);
double tr = 0;
size_t n = amp.size1();
for ( size_t i = 0; i < n; ++i ) {
tr += real(inner_prod(row(cij,i),column(amp,i)));
}
return tr;
}
struct pickColour {
PDT::Colour operator()(tcPDPtr p) const {
return p->iColour();
}
};
vector<PDT::Colour> ColourBasis::projectColour(const cPDVector& sub) const {
vector<PDT::Colour> res(sub.size());
transform(sub.begin(),sub.end(),res.begin(),pickColour());
return res;
}
vector<PDT::Colour> ColourBasis::normalOrder(const vector<PDT::Colour>& legs) const {
vector<PDT::Colour> crosslegs = legs;
if ( crosslegs[0] == PDT::Colour3 )
crosslegs[0] = PDT::Colour3bar;
else if ( crosslegs[0] == PDT::Colour3bar )
crosslegs[0] = PDT::Colour3;
if ( crosslegs[1] == PDT::Colour3 )
crosslegs[1] = PDT::Colour3bar;
else if ( crosslegs[1] == PDT::Colour3bar )
crosslegs[1] = PDT::Colour3;
int n3 = count_if(crosslegs.begin(),crosslegs.end(),matchRep(PDT::Colour3));
int n8 = count_if(crosslegs.begin(),crosslegs.end(),matchRep(PDT::Colour8));
vector<PDT::Colour> ordered(2*n3+n8,PDT::Colour8);
int i = 0;
while ( i < 2*n3 ) {
ordered[i] = PDT::Colour3;
ordered[i+1] = PDT::Colour3bar;
i+=2;
}
return ordered;
}
string ColourBasis::file(const vector<PDT::Colour>& sub) const {
string res = "";
for ( vector<PDT::Colour>::const_iterator lit = sub.begin();
lit != sub.end(); ++lit ) {
if ( *lit == PDT::Colour3 )
res += "3";
if ( *lit == PDT::Colour3bar )
res += "3bar";
if ( *lit == PDT::Colour8 )
res += "8";
}
if ( largeN() )
res += "largeN";
return res;
}
void ColourBasis::writeBasis(const string& prefix) const {
if ( didWrite )
return;
set<vector<PDT::Colour> > legs;
for ( map<cPDVector,vector<PDT::Colour> >::const_iterator lit
= theNormalOrderedLegs.begin(); lit != theNormalOrderedLegs.end(); ++lit ) {
legs.insert(lit->second);
}
string searchPath = theSearchPath;
if ( searchPath != "" )
if ( *(--searchPath.end()) != '/' )
searchPath += "/";
for ( set<vector<PDT::Colour> >::const_iterator known = legs.begin();
known != legs.end(); ++known ) {
string fname = searchPath + prefix + file(*known) + ".cdat";
ifstream check(fname.c_str());
if ( check ) continue;
ofstream out(fname.c_str());
if ( !out )
throw Exception() << "ColourBasis failed to open "
<< fname << " for storing colour basis information."
<< Exception::abortnow;
out << setprecision(18);
const symmetric_matrix<double,upper>& sp =
theScalarProducts.find(*known)->second;
write(sp,out);
if ( theCharges.find(*known) != theCharges.end() ) {
out << "#charges\n";
const map<size_t,compressed_matrix<double> >& tm =
theCharges.find(*known)->second;
const map<size_t,vector<pair<size_t,size_t> > >& tc =
theChargeNonZeros.find(*known)->second;
map<size_t,vector<pair<size_t,size_t> > >::const_iterator kc =
tc.begin();
for ( map<size_t,compressed_matrix<double> >::const_iterator k = tm.begin();
k != tm.end(); ++k, ++kc ) {
out << k->first << "\n";
write(k->second,out,kc->second);
}
const map<pair<size_t,size_t>,symmetric_matrix<double,upper> >& cm =
theCorrelators.find(*known)->second;
for ( map<pair<size_t,size_t>,symmetric_matrix<double,upper> >::const_iterator k =
cm.begin(); k != cm.end(); ++k ) {
out << k->first.first << "\n" << k->first.second << "\n";
write(k->second,out);
}
} else {
out << "#nocharges\n";
}
out << flush;
}
didWrite = true;
}
bool ColourBasis::readBasis(const vector<PDT::Colour>& legs) {
string searchPath = theSearchPath;
if ( searchPath != "" )
if ( *(--searchPath.end()) != '/' )
searchPath += "/";
string fname = searchPath + file(legs) + ".cdat";
ifstream in(fname.c_str());
if ( !in )
return false;
read(theScalarProducts[legs],in);
string tag; in >> tag;
if ( tag != "#nocharges" ) {
for ( size_t k = 0; k < legs.size(); ++k ) {
size_t i; in >> i;
read(theCharges[legs][i],in,theChargeNonZeros[legs][i]);
}
for ( size_t k = 0; k < legs.size()*(legs.size()-1)/2; ++k ) {
size_t i,j; in >> i >> j;
read(theCorrelators[legs][make_pair(i,j)],in);
}
}
readBasisDetails(legs);
return true;
}
void ColourBasis::readBasis() {
if ( didRead )
return;
string searchPath = theSearchPath;
if ( searchPath != "" )
if ( *(--searchPath.end()) != '/' )
searchPath += "/";
set<vector<PDT::Colour> > legs;
for ( map<cPDVector,vector<PDT::Colour> >::const_iterator lit
= theNormalOrderedLegs.begin(); lit != theNormalOrderedLegs.end(); ++lit )
legs.insert(lit->second);
for ( set<vector<PDT::Colour> >::const_iterator known = legs.begin();
known != legs.end(); ++known ) {
if ( theScalarProducts.find(*known) != theScalarProducts.end() )
continue;
string fname = searchPath + file(*known) + ".cdat";
if ( !readBasis(*known) )
throw Exception() << "ColourBasis failed to open "
<< fname << " for reading colour basis information."
<< Exception::abortnow;
}
didRead = true;
}
void ColourBasis::write(const symmetric_matrix<double,upper>& m, ostream& os) const {
os << m.size1() << "\n";
for ( size_t i = 0; i < m.size1(); ++i )
for ( size_t j = i; j < m.size1(); ++j )
os << m(i,j) << "\n";
os << flush;
}
void ColourBasis::read(symmetric_matrix<double,upper>& m, istream& is) {
size_t s; is >> s;
m.resize(s);
for ( size_t i = 0; i < m.size1(); ++i )
for ( size_t j = i; j < m.size1(); ++j )
is >> m(i,j);
}
void ColourBasis::write(const compressed_matrix<double>& m, ostream& os,
const vector<pair<size_t,size_t> >& nonZeros) const {
os << nonZeros.size() << "\n"
<< m.size1() << "\n"
<< m.size2() << "\n";
for ( vector<pair<size_t,size_t> >::const_iterator nz = nonZeros.begin();
nz != nonZeros.end(); ++nz )
os << nz->first << "\n" << nz->second << "\n"
<< m(nz->first,nz->second) << "\n";
os << flush;
}
void ColourBasis::read(compressed_matrix<double>& m, istream& is,
vector<pair<size_t,size_t> >& nonZeros) {
size_t nonZero, size1, size2;
is >> nonZero >> size1 >> size2;
nonZeros.resize(nonZero);
m = compressed_matrix<double>(size1,size2,nonZero);
for ( size_t k = 0; k < nonZero; ++k ) {
size_t i,j; double val;
is >> i >> j >> val;
nonZeros[k] = make_pair(i,j);
m(i,j) = val;
}
}
void ColourBasis::doinit() {
HandlerBase::doinit();
readBasis();
}
void ColourBasis::dofinish() {
HandlerBase::dofinish();
writeBasis();
}
void ColourBasis::doinitrun() {
HandlerBase::doinitrun();
readBasis();
}
void ColourBasis::persistentOutput(PersistentOStream & os) const {
os << theLargeN << theSearchPath << theNormalOrderedLegs
<< theIndexMap << theFlowMap << theOrderingStringIdentifiers << theOrderingIdentifiers;
writeBasis();
}
void ColourBasis::persistentInput(PersistentIStream & is, int) {
is >> theLargeN >> theSearchPath >> theNormalOrderedLegs
>> theIndexMap >> theFlowMap >> theOrderingStringIdentifiers >> theOrderingIdentifiers;
+ readBasis();
}
// *** Attention *** The following static variable is needed for the type
// description system in ThePEG. Please check that the template arguments
// are correct (the class and its base class), and that the constructor
// arguments are correct (the class name and the name of the dynamically
// loadable library where the class implementation can be found).
DescribeAbstractClass<ColourBasis,HandlerBase>
describeColourBasis("Herwig::ColourBasis", "HwMatchbox.so");
void ColourBasis::Init() {
static ClassDocumentation<ColourBasis> documentation
("ColourBasis is an interface to a colour basis "
"implementation.");
static Parameter<ColourBasis,string> interfaceSearchPath
("SearchPath",
"Set the search path for pre-computed colour basis data.",
&ColourBasis::theSearchPath, ".",
false, false);
static Switch<ColourBasis,bool> interfaceLargeN
("LargeN",
"Switch on or off large-N evaluation.",
&ColourBasis::theLargeN, false, false, false);
static SwitchOption interfaceLargeNOn
(interfaceLargeN,
"On",
"Work in N=infinity",
true);
static SwitchOption interfaceLargeNOff
(interfaceLargeN,
"Off",
"Work in N=3",
false);
}
diff --git a/MatrixElement/Matchbox/Utility/MatchboxXCombData.cc b/MatrixElement/Matchbox/Utility/MatchboxXCombData.cc
--- a/MatrixElement/Matchbox/Utility/MatchboxXCombData.cc
+++ b/MatrixElement/Matchbox/Utility/MatchboxXCombData.cc
@@ -1,103 +1,145 @@
// -*- C++ -*-
//
// MatchboxXCombData.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MatchboxXCombData class.
//
#include "MatchboxXCombData.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+
using namespace Herwig;
MatchboxXCombData::MatchboxXCombData()
: theCrossingSign(0.0), theCalculateTreeAmplitudes(true),
theCalculateOneLoopAmplitudes(true), theCalculateTreeME2(true),
theLastTreeME2(0.0), theCalculateOneLoopInterference(true),
theLastOneLoopInterference(0.0), theNLight(0),
theColourBasisDim(0), theNDimPhasespace(0),
theNDimAmplitude(0), theNDimInsertions(0),
theSymmetryFactor(0.0) {
flushCaches();
}
MatchboxXCombData::~MatchboxXCombData() {}
MatchboxXCombData::MatchboxXCombData(tMEPtr newME)
: theCrossingSign(0.0), theCalculateTreeAmplitudes(true),
theCalculateOneLoopAmplitudes(true), theCalculateTreeME2(true),
theLastTreeME2(0.0), theCalculateOneLoopInterference(true),
theLastOneLoopInterference(0.0), theNLight(0),
theColourBasisDim(0), theNDimPhasespace(0),
theNDimAmplitude(0), theNDimInsertions(0),
theSymmetryFactor(0.0){
flushCaches();
theMatchboxME = dynamic_ptr_cast<Ptr<MatchboxMEBase>::tptr>(newME);
theSubtractionDipole = dynamic_ptr_cast<Ptr<SubtractionDipole>::tptr>(newME);
}
Ptr<MatchboxMEBase>::tptr MatchboxXCombData::matchboxME() const {
return theMatchboxME;
}
Ptr<SubtractionDipole>::tptr MatchboxXCombData::subtractionDipole() const {
return theSubtractionDipole;
}
void MatchboxXCombData::flushCaches() {
theCalculateTreeAmplitudes = true;
theCalculateOneLoopAmplitudes = true;
theCalculateTreeME2 = true;
theCalculateOneLoopInterference = true;
for ( map<pair<int,int>,bool>::iterator f = theCalculateColourCorrelators.begin();
f != theCalculateColourCorrelators.end(); ++f )
f->second = true;
for ( map<pair<int,int>,bool>::iterator f = theCalculateLargeNColourCorrelators.begin();
f != theCalculateLargeNColourCorrelators.end(); ++f )
f->second = true;
for ( map<pair<int,int>,bool>::iterator f = theCalculateColourSpinCorrelators.begin();
f != theCalculateColourSpinCorrelators.end(); ++f )
f->second = true;
}
+void MatchboxXCombData::putCVector(PersistentOStream& os, const CVector& v) {
+ size_t n = v.size();
+ os << n;
+ for ( size_t k = 0; k < n; k++ )
+ os << v(k);
+}
+
+void MatchboxXCombData::getCVector(PersistentIStream& is, CVector& v) {
+ size_t n; is >> n;
+ v.resize(n);
+ Complex value;
+ for ( size_t k = 0; k < n; k++ ) {
+ is >> value; v(k) = value;
+ }
+}
+
+void MatchboxXCombData::putAmplitudeMap(PersistentOStream& os, const map<vector<int>,CVector>& amps) {
+ os << amps.size();
+ for ( map<vector<int>,CVector>::const_iterator a = amps.begin();
+ a != amps.end(); ++a ) {
+ os << a->first;
+ putCVector(os,a->second);
+ }
+}
+
+void MatchboxXCombData::getAmplitudeMap(PersistentIStream& is, map<vector<int>,CVector>& amps) {
+ size_t n; is >> n;
+ for ( size_t k = 0; k < n; ++k ) {
+ vector<int> hel; is >> hel;
+ getCVector(is,amps[hel]);
+ }
+}
+
void MatchboxXCombData::persistentOutput(PersistentOStream & os) const {
os << theMatchboxME << theSubtractionDipole << theCrossingMap
<< theAmplitudeToColourMap << theColourToAmplitudeMap
<< theCrossingSign << theAmplitudePartonData << ounit(theAmplitudeMomenta,GeV)
- << theCalculateTreeAmplitudes /*<< theLastAmplitudes << theLastLargeNAmplitudes*/
- << theCalculateOneLoopAmplitudes /*<< theLastOneLoopAmplitudes*/
+ << theCalculateTreeAmplitudes /* << theLastAmplitudes << theLastLargeNAmplitudes */
+ << theCalculateOneLoopAmplitudes /* << theLastOneLoopAmplitudes */
<< theCalculateTreeME2 << theLastTreeME2 << theCalculateOneLoopInterference
<< theLastOneLoopInterference << theCalculateColourCorrelators
<< theColourCorrelators << theCalculateLargeNColourCorrelators
<< theLargeNColourCorrelators << theCalculateColourSpinCorrelators
<< theColourSpinCorrelators << theNLight << theColourBasisDim
<< theNDimPhasespace << theNDimAmplitude << theNDimInsertions
<< theAmplitudeRandomNumbers << theInsertionRandomNumbers
<< theDiagramWeights << theSingularLimits// << theLastSingularLimit
<< theStandardModel << theSymmetryFactor;
+ putAmplitudeMap(os,theLastAmplitudes);
+ putAmplitudeMap(os,theLastLargeNAmplitudes);
+ putAmplitudeMap(os,theLastOneLoopAmplitudes);
}
void MatchboxXCombData::persistentInput(PersistentIStream & is, int) {
is >> theMatchboxME >> theSubtractionDipole >> theCrossingMap
>> theAmplitudeToColourMap >> theColourToAmplitudeMap
>> theCrossingSign >> theAmplitudePartonData >> iunit(theAmplitudeMomenta,GeV)
- >> theCalculateTreeAmplitudes /*>> theLastAmplitudes >> theLastLargeNAmplitudes*/
- >> theCalculateOneLoopAmplitudes /*>> theLastOneLoopAmplitudes*/
+ >> theCalculateTreeAmplitudes /* >> theLastAmplitudes >> theLastLargeNAmplitudes */
+ >> theCalculateOneLoopAmplitudes /* >> theLastOneLoopAmplitudes */
>> theCalculateTreeME2 >> theLastTreeME2 >> theCalculateOneLoopInterference
>> theLastOneLoopInterference >> theCalculateColourCorrelators
>> theColourCorrelators >> theCalculateLargeNColourCorrelators
>> theLargeNColourCorrelators >> theCalculateColourSpinCorrelators
>> theColourSpinCorrelators >> theNLight >> theColourBasisDim
>> theNDimPhasespace >> theNDimAmplitude >> theNDimInsertions
>> theAmplitudeRandomNumbers >> theInsertionRandomNumbers
>> theDiagramWeights >> theSingularLimits// >> theLastSingularLimit
>> theStandardModel >> theSymmetryFactor;
+ getAmplitudeMap(is,theLastAmplitudes);
+ getAmplitudeMap(is,theLastLargeNAmplitudes);
+ getAmplitudeMap(is,theLastOneLoopAmplitudes);
}
diff --git a/MatrixElement/Matchbox/Utility/MatchboxXCombData.h b/MatrixElement/Matchbox/Utility/MatchboxXCombData.h
--- a/MatrixElement/Matchbox/Utility/MatchboxXCombData.h
+++ b/MatrixElement/Matchbox/Utility/MatchboxXCombData.h
@@ -1,746 +1,746 @@
// -*- C++ -*-
//
// MatchboxXCombData.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2012 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_MatchboxXCombData_H
#define Herwig_MatchboxXCombData_H
//
// This is the declaration of the MatchboxXCombData class.
//
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include "ThePEG/MatrixElement/MEBase.h"
#include "Herwig++/MatrixElement/Matchbox/Base/MatchboxMEBase.fh"
#include "Herwig++/MatrixElement/Matchbox/Dipoles/SubtractionDipole.fh"
#include "Herwig++/Models/StandardModel/StandardModel.h"
-#include "ThePEG/Persistency/PersistentOStream.h"
-#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/Persistency/PersistentOStream.fh"
+#include "ThePEG/Persistency/PersistentIStream.fh"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Define complex vector from boost::uBLAS
*/
typedef boost::numeric::ublas::vector<Complex> CVector;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Define how amplitudes are stored
*/
typedef map<vector<int>,CVector> AmplitudeMap;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Define amplitude iterators
*/
typedef map<vector<int>,CVector>::iterator AmplitudeIterator;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Define amplitude const iterators
*/
typedef map<vector<int>,CVector>::const_iterator AmplitudeConstIterator;
/**
* \ingroup Matchbox
* \author Simon Platzer
*
* \brief Matchbox extensions to StandardXComb
*/
class MatchboxXCombData {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* Standard constructor.
*/
explicit MatchboxXCombData(tMEPtr newME);
/**
* Default constructor.
*/
MatchboxXCombData();
/**
* Destructor.
*/
virtual ~MatchboxXCombData();
//@}
public:
/**
* Reset all cache flags
*/
void flushCaches();
public:
/**
* Get the matrix element; may return null
*/
Ptr<MatchboxMEBase>::tptr matchboxME() const;
/**
* Get the dipole; may return null
*/
Ptr<SubtractionDipole>::tptr subtractionDipole() const;
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
const vector<int>& crossingMap() const { return theCrossingMap; }
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
vector<int>& crossingMap() { return theCrossingMap; }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
const map<size_t,size_t>& amplitudeToColourMap() const { return theAmplitudeToColourMap; }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t>& amplitudeToColourMap() { return theAmplitudeToColourMap; }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
const map<size_t,size_t>& colourToAmplitudeMap() const { return theColourToAmplitudeMap; }
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t>& colourToAmplitudeMap() { return theColourToAmplitudeMap; }
/**
* The crossing sign as filled by the last call to
* fillCrossingMap()
*/
double crossingSign() const { return theCrossingSign; }
/**
* The crossing sign as filled by the last call to
* fillCrossingMap()
*/
void crossingSign(double c) { theCrossingSign = c; }
/**
* The amplitude parton data.
*/
const cPDVector& amplitudePartonData() const { return theAmplitudePartonData; }
/**
* The amplitude parton data.
*/
cPDVector& amplitudePartonData() { return theAmplitudePartonData; }
/**
* The crossed momenta
*/
const vector<Lorentz5Momentum>& amplitudeMomenta() const { return theAmplitudeMomenta; }
/**
* The crossed momenta
*/
vector<Lorentz5Momentum>& amplitudeMomenta() { return theAmplitudeMomenta; }
/**
* True, if the the tree level amplitudes need to be calculated
*/
bool calculateTreeAmplitudes() const { return theCalculateTreeAmplitudes; }
/**
* The amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
const map<vector<int>,CVector>& lastAmplitudes() const { return theLastAmplitudes; }
/**
* True, if the the tree level amplitudes need to be calculated
*/
void haveTreeAmplitudes(bool f = true) { theCalculateTreeAmplitudes = !f; }
/**
* The amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector>& lastAmplitudes() { return theLastAmplitudes; }
/**
* The leading N amplitude values which have been
* contributing to the last call of prepareAmplitudes.
*/
const map<vector<int>,CVector>& lastLargeNAmplitudes() const { return theLastLargeNAmplitudes; }
/**
* The leading N amplitude values which have been
* contributing to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector>& lastLargeNAmplitudes() { return theLastLargeNAmplitudes; }
/**
* True, if the the one-loop amplitudes need to be calculated
*/
bool calculateOneLoopAmplitudes() const { return theCalculateOneLoopAmplitudes; }
/**
* The one-loop amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
const map<vector<int>,CVector>& lastOneLoopAmplitudes() const { return theLastOneLoopAmplitudes; }
/**
* True, if the the one-loop amplitudes need to be calculated
*/
void haveOneLoopAmplitudes(bool f = true) { theCalculateOneLoopAmplitudes = !f; }
/**
* The one-loop amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector>& lastOneLoopAmplitudes() { return theLastOneLoopAmplitudes; }
/**
* True, if the tree-level matrix element squared needs to be
* calculated.
*/
bool calculateTreeME2() const { return theCalculateTreeME2; }
/**
* The last tree-level matrix element squared
*/
double lastTreeME2() const { return theLastTreeME2; }
/**
* The last tree-level matrix element squared
*/
void lastTreeME2(double v) {
theLastTreeME2 = v; theCalculateTreeME2 = false;
}
/**
* True, if the one-loop/tree-level interference.
* be calculated.
*/
bool calculateOneLoopInterference() const { return theCalculateOneLoopInterference; }
/**
* The last one-loop/tree-level interference.
*/
double lastOneLoopInterference() const { return theLastOneLoopInterference; }
/**
* The last one-loop/tree-level interference.
*/
void lastOneLoopInterference(double v) {
theLastOneLoopInterference = v; theCalculateOneLoopAmplitudes = false;
}
/**
* True, if the indexed colour correlated matrix element needs to be
* calculated.
*/
bool calculateColourCorrelator(pair<int,int> ij) const {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
map<pair<int,int>,bool>::const_iterator f =
theCalculateColourCorrelators.find(ij);
if ( f == theCalculateColourCorrelators.end() )
return true;
return f->second;
}
/**
* The colour correlated matrix element.
*/
double lastColourCorrelator(pair<int,int> ij) const {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
map<pair<int,int>,double>::const_iterator v =
theColourCorrelators.find(ij);
if ( v == theColourCorrelators.end() )
return 0.;
return v->second;
}
/**
* The colour correlated matrix element.
*/
void lastColourCorrelator(pair<int,int> ij, double v) {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
theColourCorrelators[ij] = v;
theCalculateColourCorrelators[ij] = false;
}
/**
* True, if the indexed large-N colour correlated matrix element needs to be
* calculated.
*/
bool calculateLargeNColourCorrelator(pair<int,int> ij) const {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
map<pair<int,int>,bool>::const_iterator f =
theCalculateLargeNColourCorrelators.find(ij);
if ( f == theCalculateLargeNColourCorrelators.end() )
return true;
return f->second;
}
/**
* The large-N colour correlated matrix element.
*/
double lastLargeNColourCorrelator(pair<int,int> ij) const {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
map<pair<int,int>,double>::const_iterator v =
theLargeNColourCorrelators.find(ij);
if ( v == theLargeNColourCorrelators.end() )
return 0.;
return v->second;
}
/**
* The large-N colour correlated matrix element.
*/
void lastLargeNColourCorrelator(pair<int,int> ij, double v) {
if ( ij.first > ij.second )
swap(ij.first,ij.second);
theLargeNColourCorrelators[ij] = v;
theCalculateLargeNColourCorrelators[ij] = false;
}
/**
* True, if the indexed colour/spin correlated matrix element needs to be
* calculated.
*/
bool calculateColourSpinCorrelator(const pair<int,int>& ij) const {
map<pair<int,int>,bool>::const_iterator f =
theCalculateColourSpinCorrelators.find(ij);
if ( f == theCalculateColourSpinCorrelators.end() )
return true;
return f->second;
}
/**
* The colour/spin correlated matrix element.
*/
Complex lastColourSpinCorrelator(const pair<int,int>& ij) const {
map<pair<int,int>,Complex>::const_iterator v =
theColourSpinCorrelators.find(ij);
if ( v == theColourSpinCorrelators.end() )
return 0.;
return v->second;
}
/**
* The colour/spin correlated matrix element.
*/
void lastColourSpinCorrelator(const pair<int,int>& ij, Complex v) {
theColourSpinCorrelators[ij] = v;
theCalculateColourSpinCorrelators[ij] = false;
}
/**
* Return the number of light flavours to be considered for this process.
*/
unsigned int nLight() const { return theNLight; }
/**
* Set the number of light flavours to be considered for this process.
*/
void nLight(unsigned int n) { theNLight = n; }
/**
* Get the dimensionality of the colour basis for this process.
*/
size_t colourBasisDim() const { return theColourBasisDim; }
/**
* Set the dimensionality of the colour basis for this process.
*/
void colourBasisDim(size_t d) { theColourBasisDim = d; }
/**
* Return the number of degrees of freedom required by the phase space generator
*/
int nDimPhasespace() const { return theNDimPhasespace; }
/**
* Set the number of degrees of freedom required by the phase space generator
*/
void nDimPhasespace(int d) { theNDimPhasespace = d; }
/**
* Return the number of degrees of freedom required by the amplitude
*/
int nDimAmplitude() const { return theNDimAmplitude; }
/**
* Set the number of degrees of freedom required by the amplitude
*/
void nDimAmplitude(int d) { theNDimAmplitude = d; }
/**
* Return the number of degrees of freedom required by the insertion operators
*/
int nDimInsertions() const { return theNDimInsertions; }
/**
* Set the number of degrees of freedom required by the insertion operators
*/
void nDimInsertions(int d) { theNDimInsertions = d; }
/**
* Get the additional random numbers required by the amplitude
*/
const vector<double>& amplitudeRandomNumbers() const { return theAmplitudeRandomNumbers; }
/**
* Access the additional random numbers required by the amplitude
*/
vector<double>& amplitudeRandomNumbers() { return theAmplitudeRandomNumbers; }
/**
* Get the additional random numbers required by the insertion operator
*/
const vector<double>& insertionRandomNumbers() const { return theInsertionRandomNumbers; }
/**
* Access the additional random numbers required by the insertion operator
*/
vector<double>& insertionRandomNumbers() { return theInsertionRandomNumbers; }
/**
* Return the diagram weights indexed by diagram id.
*/
const map<int,double>& diagramWeights() const { return theDiagramWeights; }
/**
* Access the diagram weights indexed by diagram id.
*/
map<int,double>& diagramWeights() { return theDiagramWeights; }
/**
* Return the singular limits
*/
const set<pair<size_t,size_t> >& singularLimits() const { return theSingularLimits; }
/**
* Access the singular limits
*/
set<pair<size_t,size_t> >& singularLimits() { return theSingularLimits; }
/**
* Return the last matched singular limit.
*/
const set<pair<size_t,size_t> >::const_iterator& lastSingularLimit() const { return theLastSingularLimit; }
/**
* Access the last matched singular limit.
*/
set<pair<size_t,size_t> >::const_iterator& lastSingularLimit() { return theLastSingularLimit; }
/**
* Set the Herwig++ StandardModel object
*/
void hwStandardModel(Ptr<StandardModel>::tcptr sm) { theStandardModel = sm; }
/**
* Get the Herwig++ StandardModel object
*/
Ptr<StandardModel>::tcptr hwStandardModel() const { return theStandardModel; }
/**
* Return the symmetry factor
*/
double symmetryFactor() const { return theSymmetryFactor; }
/**
* Set the symmetry factor
*/
void symmetryFactor(double f) { theSymmetryFactor = f; }
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
+
+ /**
+ * Put a CVector to the persistent ostream
+ */
+ static void putCVector(PersistentOStream&, const CVector&);
+
+ /**
+ * Get a CVector from the persistent istream
+ */
+ static void getCVector(PersistentIStream&, CVector&);
+
+ /**
+ * Put an amplitude map to the persistent ostream
+ */
+ static void putAmplitudeMap(PersistentOStream&, const map<vector<int>,CVector>&);
+
+ /**
+ * Get an amplitude map from the persistent istream
+ */
+ static void getAmplitudeMap(PersistentIStream&, map<vector<int>,CVector>&);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
MatchboxXCombData & operator=(const MatchboxXCombData &);
private:
/**
* The matrix element
*/
Ptr<MatchboxMEBase>::tptr theMatchboxME;
/**
* The dipole
*/
Ptr<SubtractionDipole>::tptr theSubtractionDipole;
/**
* The crossing information as filled by the last call to
* fillCrossingMap()
*/
vector<int> theCrossingMap;
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t> theAmplitudeToColourMap;
/**
* The colour crossing information as filled by the last call to
* fillCrossingMap()
*/
map<size_t,size_t> theColourToAmplitudeMap;
/**
* The crossing sign as filled by the last call to
* fillCrossingMap()
*/
double theCrossingSign;
/**
* The amplitude parton data.
*/
cPDVector theAmplitudePartonData;
/**
* The ccrossed momenta
*/
vector<Lorentz5Momentum> theAmplitudeMomenta;
/**
* True, if the the tree level amplitudes need to be calculated
*/
bool theCalculateTreeAmplitudes;
/**
* The amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector> theLastAmplitudes;
/**
* The leading N amplitude values which have been
* contributing to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector> theLastLargeNAmplitudes;
/**
* True, if the the one-loop amplitudes need to be calculated
*/
bool theCalculateOneLoopAmplitudes;
/**
* The one-loop amplitude values which have been contributing
* to the last call of prepareAmplitudes.
*/
map<vector<int>,CVector> theLastOneLoopAmplitudes;
/**
* True, if the tree-level matrix element squared needs to be
* calculated.
*/
bool theCalculateTreeME2;
/**
* The last tree-level matrix element squared
*/
double theLastTreeME2;
/**
* True, if the one-loop/tree-level interference.
* be calculated.
*/
bool theCalculateOneLoopInterference;
/**
* The last one-loop/tree-level interference.
*/
double theLastOneLoopInterference;
/**
* True, if the indexed colour correlated matrix element needs to be
* calculated.
*/
map<pair<int,int>,bool> theCalculateColourCorrelators;
/**
* The colour correlated matrix element.
*/
map<pair<int,int>,double> theColourCorrelators;
/**
* True, if the indexed large-N colour correlated matrix element needs to be
* calculated.
*/
map<pair<int,int>,bool> theCalculateLargeNColourCorrelators;
/**
* The large-N colour correlated matrix element.
*/
map<pair<int,int>,double> theLargeNColourCorrelators;
/**
* True, if the indexed colour/spin correlated matrix element needs to be
* calculated.
*/
map<pair<int,int>,bool> theCalculateColourSpinCorrelators;
/**
* The colour/spin correlated matrix element.
*/
map<pair<int,int>,Complex> theColourSpinCorrelators;
/**
* The number of light flavours to be considered for this process.
*/
unsigned int theNLight;
/**
* The dimensionality of the colour basis for this process.
*/
size_t theColourBasisDim;
/**
* The number of degrees of freedom required by the phase space generator
*/
int theNDimPhasespace;
/**
* The number of degrees of freedom required by the amplitude
*/
int theNDimAmplitude;
/**
* The number of degrees of freedom required by the insertion operators
*/
int theNDimInsertions;
/**
* Additional random numbers required by the amplitude
*/
vector<double> theAmplitudeRandomNumbers;
/**
* Additional random numbers required by the insertion operator
*/
vector<double> theInsertionRandomNumbers;
/**
* The diagram weights indexed by diagram id.
*/
map<int,double> theDiagramWeights;
/**
* If not empty, the entries here serve to limit phasespace
* generation to the corresponding collinear limits, or soft limits
* if both pair entries are equal.
*/
set<pair<size_t,size_t> > theSingularLimits;
/**
* The last matched singular limit.
*/
set<pair<size_t,size_t> >::const_iterator theLastSingularLimit;
/**
* The Herwig++ StandardModel object
*/
Ptr<StandardModel>::tcptr theStandardModel;
/**
* The symmetry factor
*/
double theSymmetryFactor;
};
- inline PersistentOStream& operator<<(PersistentOStream& os,
- const CVector& v) {
- size_t n = v.size();
- os << n;
- for ( size_t k = 0; k < n; k++ )
- os << v(k);
- return os;
- }
-
- inline PersistentIStream& operator>>(PersistentIStream& is,
- CVector& v) {
- size_t n; is >> n;
- v.resize(n);
- Complex value;
- for ( size_t k = 0; k < n; k++ ) {
- is >> value; v(k) = value;
- }
- return is;
- }
-
}
#endif /* Herwig_MatchboxXCombData_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:31 PM (1 d, 8 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805838
Default Alt Text
(288 KB)

Event Timeline