diff --git a/Config/Containers.cc b/Config/Containers.cc
deleted file mode 100644
--- a/Config/Containers.cc
+++ /dev/null
@@ -1,38 +0,0 @@
-// -*- C++ -*-
-//
-// Containers.cc is a part of ThePEG - Toolkit for HEP Event Generation
-// Copyright (C) 1999-2011 Leif Lonnblad
-//
-// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
-// Please respect the MCnet academic guidelines, see GUIDELINES for details.
-//
-
-// This file contains the implementations of the container
-// declarations in Containers.h.
-
-#include "ThePEG.h"
-#include "ThePEG/PDT/ParticleData.h"
-#include "ThePEG/PDT/MatcherBase.h"
-#include "ThePEG/PDT/DecayMode.h"
-#include "ThePEG/Interface/InterfacedBase.h"
-#include "ThePEG/Repository/EventGenerator.h"
-
-#ifdef ThePEG_TEMPLATES_IN_CC_FILE
-// #include "Containers.tcc"
-#endif
-
-using namespace ThePEG;
-
-ThePEG_IMPLEMENT_SET(PDPTr,ParticleDataSet)
-ThePEG_IMPLEMENT_SET(PMPtr,MatcherSet)
-ThePEG_IMPLEMENT_SET(DMPtr,DecayModeSet)
-ThePEG_IMPLEMENT_SET(tDMPtr,DecaySet)
-ThePEG_IMPLEMENT_SET(IBPtr,ObjectSet)
-ThePEG_IMPLEMENT_SET(IBPtr,DependencySet)
-ThePEG_IMPLEMENT_MAP(long,PDPtr,ParticleMap)
-ThePEG_IMPLEMENT_MAP(string,IBPtr,ObjectMap)
-ThePEG_IMPLEMENT_MAP(IBPtr,DependencySet,DependencyMap)
-ThePEG_IMPLEMENT_MAP(string,const InterfaceBase *,InterfaceMap)
-ThePEG_IMPLEMENT_MAP(string,EGPtr,GeneratorMap)
-ThePEG_IMPLEMENT_SET(string,StringSet)
-// typedef set<const InterfaceBase *> InterfaceSet;
diff --git a/Config/Makefile.am b/Config/Makefile.am
--- a/Config/Makefile.am
+++ b/Config/Makefile.am
@@ -1,18 +1,18 @@
-mySOURCES = ThePEG.cc Containers.cc
+mySOURCES = ThePEG.cc
DOCFILES = ThePEG.h Constants.h Containers.h Pointers.h \
Unitsystem.h algorithm.h \
std.h Complex.h TemplateTools.h \
PhysicalQty.h PhysicalQtyOps.h PhysicalQtyComplex.h
if HAVE_HEPMC
DOCFILES += HepMCHelper.h
endif
INCLUDEFILES = $(DOCFILES)
noinst_LTLIBRARIES = libThePEGConfig.la
libThePEGConfig_la_SOURCES = $(mySOURCES) $(INCLUDEFILES)
include $(top_srcdir)/Config/Makefile.aminclude
diff --git a/Config/std.h b/Config/std.h
--- a/Config/std.h
+++ b/Config/std.h
@@ -1,493 +1,180 @@
// -*- C++ -*-
//
// std.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_std_H
#define ThePEG_std_H
/** \file
* This file introduces a number of <code>std::</code> classes into
* the ThePEG namespace. Also introduces some useful functions for
* standard library classes.
*
* Do not make changes in this file. If you want to use alternatives
* to the <code>std::</code> classes in ThePEG, edit a copy of this
* file and include it in an alternative config file which can be
* included in the main ThePEG.h config file using the macro
* <code>ThePEG_ALTERNATE_CONFIG</code>.
*/
#include <map>
#include <set>
#include <string>
#include <vector>
#include <list>
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <stack>
#include <utility>
#include <typeinfo>
#include <stdexcept>
#include <cmath>
namespace std {
/** @cond TRAITSPECIALIZATIONS */
/**
* This specialization of the std::less class is needed in order to be
* able use put pointers to type_info objects as keys in maps and
* sets.
*/
template <>
struct less<const type_info *> :
public binary_function<const type_info *, const type_info *, bool>
{
/**
* This is the function called when comparing two pointers to
* type_info.
*/
bool operator()(const type_info * x, const type_info * y) const {
return x->before(*y); }
};
/** @endcond */
}
namespace ThePEG {
using std::deque;
using std::stack;
using std::vector;
using std::multiset;
using std::set;
using std::map;
using std::list;
using std::multimap;
using std::pair;
using std::make_pair;
using std::less;
using std::string;
using std::type_info;
using std::exception;
using std::range_error;
using std::ios;
using std::ostream;
using std::istream;
using std::ofstream;
using std::ifstream;
using std::ostringstream;
using std::istringstream;
using std::cin;
using std::cout;
using std::cerr;
using std::endl;
using std::flush;
using std::setprecision;
using std::setw;
using std::swap;
using std::min;
using std::max;
using std::mem_fun;
using std::sqrt;
//using std::pow;
using std::abs;
using std::atan2;
using std::isfinite;
/** Powers - standard or non-standard */
template <class ExponentT>
inline constexpr double pow(double x, ExponentT p) {
return std::pow(x,double(p));
}
/** Square root of an integer. */
inline double sqrt(int x) {
return std::sqrt(double(x));
}
/** Check if a given object is a part of a container. */
template <typename Container, typename Key>
inline bool member(const Container & c, const Key & k) {
return c.find(k) != c.end();
}
/** Check if a given object is a part of a vector. */
template <typename T, typename Key>
inline bool member(const vector<T> & v, const Key & k) {
for ( typename vector<T>::const_iterator i = v.begin(); i != v.end(); ++i )
if ( *i == k ) return true;
return false;
// return find(v.begin(), v.end(), k) != v.end();
}
/** Return an insert iterator for a given container. */
template <typename Cont>
inline std::insert_iterator<Cont> inserter(Cont & c) {
return std::insert_iterator<Cont>(c, c.end());
}
/** Return an insert iterator for a given vector. Overrides the
* general version. */
template <typename T, typename A>
inline std::back_insert_iterator< vector<T,A> > inserter(vector<T,A> & v) {
return back_inserter(v);
}
/** Return an insert iterator for a given vector. Overrides the
* general version. */
template <typename T, typename A>
inline std::back_insert_iterator< deque<T,A> > inserter(deque<T,A> & v) {
return back_inserter(v);
}
/** Stream manipulator setting an ostream to left-adjust its ouput. */
inline ostream& left(ostream& os) {
os.setf(ios::left, ios::adjustfield);
return os;
}
/** Stream manipulator setting an ostream to right-adjust its ouput. */
inline ostream& right(ostream& os) {
os.setf(ios::right, ios::adjustfield);
return os;
}
}
-#ifndef ThePEG_WRAP_STL_CONTAINERS
-
/** Macro for declaring a set. */
#define ThePEG_DECLARE_SET(VALTYPE,NAME) \
/** A set of VALTYPE. */ \
typedef set<VALTYPE, less<VALTYPE> > NAME
/** Macro for declaring a multiset. */
#define ThePEG_DECLARE_MULTISET(VALTYPE,NAME) \
/** A multiset of VALTYPE. */ \
typedef multiset<VALTYPE, less<VALTYPE> > NAME
/** Macro for declaring a map. */
#define ThePEG_DECLARE_MAP(KEYTYPE,VALTYPE,NAME) \
/** A map of VALTYPE indexed by KEYTYPE. */ \
typedef map<KEYTYPE, VALTYPE, less<KEYTYPE> > NAME
-/** Macro for implementing a set. */
-#define ThePEG_IMPLEMENT_SET(VALTYPE,NAME)
-
-/** Macro for implementing a multiset. */
-#define ThePEG_IMPLEMENT_MULTISET(VALTYPE,NAME)
-
-/** Macro for implementing a map. */
-#define ThePEG_IMPLEMENT_MAP(KEYTYPE,VALTYPE,NAME)
-
-#else
-
-/** Macro for declaring a set. */
-#define ThePEG_DECLARE_SET(VALTYPE,NAME) \
-class NAME : public set<VALTYPE, less<VALTYPE> > { \
-public: \
- typedef set<VALTYPE, less<VALTYPE> > SETTYPE; \
- NAME(); \
- explicit NAME(const key_compare & c, \
- const allocator_type & a = allocator_type()); \
- template <typename InputIterator> \
- NAME(InputIterator first, InputIterator last) \
- : SETTYPE(first, last) {} \
- template <typename InputIterator> \
- NAME(InputIterator first, InputIterator last, const key_compare & c, \
- const allocator_type & a = allocator_type()) \
- : SETTYPE(first, last, c, a) {} \
- NAME(const SETTYPE & s); \
- NAME(const NAME & s); \
- ~NAME(); \
- NAME & operator=(const NAME &); \
- NAME & operator=(const SETTYPE &); \
- pair<iterator,bool> insert(const value_type & x); \
- iterator insert(iterator position, const value_type & x); \
- template <typename InputIterator> \
- void insert(InputIterator first, InputIterator last) { \
- SETTYPE::insert(first, last); \
- } \
- void erase(iterator position); \
- size_type erase(const key_type & x); \
- void erase(iterator first, iterator last); \
- void clear(); \
- iterator find(const key_type & x) const; \
- size_type count(const key_type & x) const; \
- iterator lower_bound(const key_type & x) const; \
- iterator upper_bound(const key_type & x) const; \
- pair<iterator,iterator> equal_range(const key_type & x) const; \
-}
-
-
-/** Macro for declaring a multiset. */
-#define ThePEG_DECLARE_MULTISET(VALTYPE,NAME) \
-class NAME : \
- public multiset<VALTYPE, less<VALTYPE> > { \
-public: \
- typedef multiset<VALTYPE, less<VALTYPE> > SETTYPE;\
- NAME(); \
- explicit NAME(const key_compare & c, \
- const allocator_type & a = allocator_type()); \
- template <typename InputIterator> \
- NAME(InputIterator first, InputIterator last) \
- : SETTYPE(first, last) {} \
- template <typename InputIterator> \
- NAME(InputIterator first, InputIterator last, const key_compare & c, \
- const allocator_type & a = allocator_type()) \
- : SETTYPE(first, last, c, a) {} \
- NAME(const SETTYPE & s); \
- NAME(const NAME & s); \
- ~NAME(); \
- NAME & operator=(const NAME &); \
- NAME & operator=(const SETTYPE &); \
- iterator insert(const value_type & x); \
- iterator insert(iterator position, const value_type & x); \
- template <typename InputIterator> \
- void insert(InputIterator first, InputIterator last) { \
- SETTYPE::insert(first, last); \
- } \
- void erase(iterator position); \
- size_type erase(const key_type & x); \
- void erase(iterator first, iterator last); \
- void clear(); \
- iterator find(const key_type & x) const; \
- size_type count(const key_type & x) const; \
- iterator lower_bound(const key_type & x) const; \
- iterator upper_bound(const key_type & x) const; \
- pair<iterator,iterator> equal_range(const key_type & x) const; \
-}
-
-
-/** Macro for declaring a map. */
-#define ThePEG_DECLARE_MAP(KEYTYPE,VALTYPE,NAME) \
-class NAME : \
- public map<KEYTYPE, VALTYPE, less<KEYTYPE> > { \
-public: \
- typedef map<KEYTYPE, VALTYPE, less<KEYTYPE> > MAPTYPE; \
- NAME(); \
- explicit NAME(const key_compare & c, \
- const allocator_type & a = allocator_type()); \
- template <typename InputIterator> \
- NAME(InputIterator first, InputIterator last) \
- : MAPTYPE(first, last) {} \
- template <typename InputIterator> \
- NAME(InputIterator first, InputIterator last, const key_compare & c, \
- const allocator_type & a = allocator_type()) \
- : MAPTYPE(first, last, c, a) {} \
- NAME(const NAME & s); \
- NAME(const MAPTYPE & s); \
- ~NAME(); \
- NAME & operator=(const NAME &); \
- NAME & operator=(const MAPTYPE &); \
- data_type & operator[](const key_type & k); \
- pair<iterator,bool> insert(const value_type & x); \
- iterator insert(iterator position, const value_type & x); \
- template <typename InputIterator> \
- void insert(InputIterator first, InputIterator last) { \
- MAPTYPE::insert(first, last); \
- } \
- void erase(iterator position); \
- size_type erase(const key_type & x); \
- void erase(iterator first, iterator last); \
- void clear(); \
- iterator find(const key_type & x); \
- const_iterator find(const key_type & x) const; \
- size_type count(const key_type & x) const; \
- iterator lower_bound(const key_type & x); \
- const_iterator lower_bound(const key_type & x) const; \
- iterator upper_bound(const key_type & x); \
- const_iterator upper_bound(const key_type & x) const; \
- pair<iterator,iterator> equal_range(const key_type & x); \
- pair<const_iterator,const_iterator> \
- equal_range(const key_type & x) const; \
-}
-
-
-/** Macro for implementing a set. */
-#define ThePEG_IMPLEMENT_SET(VALTYPE,NAME) \
-NAME::NAME() {} \
-NAME::NAME(const key_compare & c, const allocator_type & a) \
- :SETTYPE(c, a) {} \
-NAME::NAME(const NAME & x) : SETTYPE(x) {} \
-NAME::NAME(const SETTYPE & x) : SETTYPE(x) {} \
-NAME::~NAME() {} \
-NAME & NAME::operator=(const NAME & x) { \
- SETTYPE::operator=(x); \
- return *this; \
-} \
-NAME & NAME::operator=(const SETTYPE & x) { \
- SETTYPE::operator=(x); \
- return *this; \
-} \
-pair<NAME::iterator,bool> NAME::insert(const value_type & x) { \
- return SETTYPE::insert(x); \
-} \
-NAME::iterator NAME::insert(iterator position, const value_type & x) { \
- return SETTYPE::insert(position, x); \
-} \
-void NAME::erase(iterator position) { \
- SETTYPE::erase(position); \
-} \
-NAME::size_type NAME::erase(const key_type & x) { \
- return SETTYPE::erase(x); \
-} \
-void NAME::erase(iterator first, iterator last) { \
- SETTYPE::erase(first, last); \
-} \
-void NAME::clear() { \
- SETTYPE::clear(); \
-} \
-NAME::iterator NAME::find(const key_type & x) const { \
- return SETTYPE::find(x); \
-} \
-NAME::size_type NAME::count(const key_type & x) const { \
- return SETTYPE::count(x); \
-} \
-NAME::iterator NAME::lower_bound(const key_type & x) const { \
- return SETTYPE::lower_bound(x); \
-} \
-NAME::iterator NAME::upper_bound(const key_type & x) const { \
- return SETTYPE::upper_bound(x); \
-} \
-pair<NAME::iterator,NAME::iterator> \
-NAME::equal_range(const key_type & x) const { \
- return SETTYPE::equal_range(x); \
-} \
-
-
-/** Macro for implementing a multiset. */
-#define ThePEG_IMPLEMENT_MULTISET(VALTYPE,NAME) \
-NAME::NAME() {} \
-NAME::NAME(const key_compare & c, const allocator_type & a) \
- :SETTYPE(c, a) {} \
-NAME::NAME(const NAME & x) : SETTYPE(x) {} \
-NAME::NAME(const SETTYPE & x) : SETTYPE(x) {} \
-NAME::~NAME() {} \
-NAME & NAME::operator=(const NAME & x) { \
- SETTYPE::operator=(x); \
- return *this; \
-} \
-NAME & NAME::operator=(const SETTYPE & x) { \
- SETTYPE::operator=(x); \
- return *this; \
-} \
-NAME::iterator NAME::insert(const value_type & x) { \
- return SETTYPE::insert(x); \
-} \
-NAME::iterator NAME::insert(iterator position, const value_type & x) { \
- return SETTYPE::insert(position, x); \
-} \
-void NAME::erase(iterator position) { \
- SETTYPE::erase(position); \
-} \
-NAME::size_type NAME::erase(const key_type & x) { \
- return SETTYPE::erase(x); \
-} \
-void NAME::erase(iterator first, iterator last) { \
- SETTYPE::erase(first, last); \
-} \
-void NAME::clear() { \
- SETTYPE::clear(); \
-} \
-NAME::iterator NAME::find(const key_type & x) const { \
- return SETTYPE::find(x); \
-} \
-NAME::size_type NAME::count(const key_type & x) const { \
- return SETTYPE::count(x); \
-} \
-NAME::iterator NAME::lower_bound(const key_type & x) const { \
- return SETTYPE::lower_bound(x); \
-} \
-NAME::iterator NAME::upper_bound(const key_type & x) const { \
- return SETTYPE::upper_bound(x); \
-} \
-pair<NAME::iterator,NAME::iterator> \
-NAME::equal_range(const key_type & x) const { \
- return SETTYPE::equal_range(x); \
-} \
-
-
-/** Macro for implementing a map. */
-#define ThePEG_IMPLEMENT_MAP(KEYTYPE,VALTYPE,NAME) \
-NAME::NAME() {} \
-NAME::NAME(const key_compare & c, const allocator_type & a) \
- :MAPTYPE(c, a) {} \
-NAME::NAME(const NAME & x) : MAPTYPE(x) {} \
-NAME::NAME(const MAPTYPE & x) : MAPTYPE(x) {} \
-NAME::~NAME() {} \
-NAME & NAME::operator=(const NAME & x) { \
- MAPTYPE::operator=(x); \
- return *this; \
-} \
-NAME & NAME::operator=(const MAPTYPE & x) { \
- MAPTYPE::operator=(x); \
- return *this; \
-} \
-pair<NAME::iterator,bool> NAME::insert(const value_type & x) { \
- return MAPTYPE::insert(x); \
-} \
-NAME::iterator NAME::insert(iterator position, const value_type & x) { \
- return MAPTYPE::insert(position, x); \
-} \
-void NAME::erase(iterator position) { \
- MAPTYPE::erase(position); \
-} \
-NAME::size_type NAME::erase(const key_type & x) { \
- return MAPTYPE::erase(x); \
-} \
-void NAME::erase(iterator first, iterator last) { \
- MAPTYPE::erase(first, last); \
-} \
-void NAME::clear() { \
- MAPTYPE::clear(); \
-} \
-NAME::iterator NAME::find(const key_type & x) { \
- return MAPTYPE::find(x); \
-} \
-NAME::const_iterator NAME::find(const key_type & x) const { \
- return MAPTYPE::find(x); \
-} \
-NAME::size_type NAME::count(const key_type & x) const { \
- return MAPTYPE::count(x); \
-} \
-NAME::iterator NAME::lower_bound(const key_type & x) { \
- return MAPTYPE::lower_bound(x); \
-} \
-NAME::const_iterator NAME::lower_bound(const key_type & x) const { \
- return MAPTYPE::lower_bound(x); \
-} \
-NAME::iterator NAME::upper_bound(const key_type & x) { \
- return MAPTYPE::upper_bound(x); \
-} \
-NAME::const_iterator NAME::upper_bound(const key_type & x) const { \
- return MAPTYPE::upper_bound(x); \
-} \
-pair<NAME::iterator,NAME::iterator> \
-NAME::equal_range(const key_type & x) { \
- return MAPTYPE::equal_range(x); \
-} \
-pair<NAME::const_iterator,NAME::const_iterator> \
-NAME::equal_range(const key_type & x) const { \
- return MAPTYPE::equal_range(x); \
-} \
-NAME::data_type & NAME::operator[](const key_type & k) { \
- return MAPTYPE::operator[](k); \
-} \
-
-#endif
-
-// #include "std.icc"
-#ifndef ThePEG_TEMPLATES_IN_CC_FILE
-// #include "std.tcc"
-#endif
-
#endif /* ThePEG_std_H */
diff --git a/EventRecord/Particle.cc b/EventRecord/Particle.cc
--- a/EventRecord/Particle.cc
+++ b/EventRecord/Particle.cc
@@ -1,524 +1,521 @@
// -*- C++ -*-
//
// Particle.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG 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 member functions of
// the Particle class.
//
#include "Particle.h"
#include "ThePEG/EventRecord/Step.h"
#include "ThePEG/EventRecord/Event.h"
#include "ThePEG/EventRecord/ColourLine.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Config/algorithm.h"
#include "ThePEG/EventRecord/ParticleTraits.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include <iostream>
#include <iomanip>
#include <cctype>
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
#include "Particle.tcc"
#endif
using namespace ThePEG;
Particle::ParticleRep::ParticleRep(const ParticleRep & p)
: theParents(p.theParents), theChildren(p.theChildren),
thePrevious(p.thePrevious), theNext(p.theNext),
theBirthStep(p.theBirthStep), theVertex(p.theVertex),
theLifeLength(p.theLifeLength), theScale(p.theScale),
theVetoScale(p.theVetoScale), theNumber(p.theNumber),
theExtraInfo(p.theExtraInfo.size()) {
if ( p.theColourInfo )
theColourInfo = dynamic_ptr_cast<CBPtr>(p.theColourInfo->clone());
if ( p.theSpinInfo )
theSpinInfo = dynamic_ptr_cast<SpinPtr>(p.theSpinInfo->clone());
for ( int i = 0, N = p.theExtraInfo.size(); i < N; ++i )
theExtraInfo[i] = p.theExtraInfo[i]->clone();
}
Particle::Particle(const Particle & p)
: Base(p), theData(p.theData), theMomentum(p.theMomentum), theRep(p.theRep) {
if ( p.theRep ) {
theRep = new ParticleRep(*p.theRep);
theRep->theParents.clear();
}
}
Particle::~Particle() {
if ( theRep ) {
if ( colourLine() ) colourLine()->removeColoured(this);
if ( antiColourLine() ) antiColourLine()->removeAntiColoured(this);
delete theRep;
}
theRep = 0;
theData = cEventPDPtr();
}
void Particle::initFull() {
if ( theRep ) return;
theRep = new ParticleRep;
Energy width = data().generateWidth(mass());
if ( width > ZERO ) {
Time lifetime = data().generateLifeTime(mass(), width);
theRep->theLifeLength.setTau(lifetime);
theRep->theLifeLength.
setVect((momentum().vect()*(lifetime /
max(mass(), Constants::epsilon*GeV))));
theRep->theLifeLength.rescaleEnergy();
}
}
PPtr Particle::clone() const {
return ptr_new<PPtr>(*this);
}
void Particle::rebind(const EventTranslationMap & trans) {
for ( ParticleVector::iterator pit = rep().theChildren.begin();
pit != rep().theChildren.end(); ++pit ) *pit = trans.translate(*pit);
for ( tParticleVector::iterator pit = rep().theParents.begin();
pit != rep().theParents.end(); ++pit ) *pit = trans.translate(*pit);
rep().thePrevious = trans.translate(rep().thePrevious);
rep().theNext = trans.translate(rep().theNext);
if ( hasColourInfo() ) colourInfo()->rebind(trans);
if ( spinInfo() ) spinInfo()->rebind(trans);
rep().theBirthStep = trans.translate(rep().theBirthStep);
for ( EIVector::const_iterator ie = rep().theExtraInfo.begin();
ie != rep().theExtraInfo.end(); ++ie ) (**ie).rebind(trans);
}
tParticleSet Particle::siblings() const {
tParticleSet theSiblings;
for ( tParticleVector::const_iterator pit = parents().begin();
pit != parents().end(); ++pit )
theSiblings.insert((*pit)->children().begin(), (*pit)->children().end());
theSiblings.erase(const_cast<Particle *>(this));
return theSiblings;
}
void Particle::colourNeighbour(tPPtr p, bool anti) {
tColinePtr line = colourLine(!anti);
if ( !line ) line = ColourLine::create(this, !anti);
line->addColoured(p, anti);
}
void Particle::outgoingColour(tPPtr p, bool anti) {
tColinePtr line = colourLine(anti);
if ( !line ) line = ColourLine::create(this, anti);
line->addColoured(p, anti);
}
tPPtr Particle::incomingColour(bool anti) const {
if ( !hasColourInfo() ) return tPPtr();
tColinePtr line = colourLine(anti);
if ( !line ) return tPPtr();
for ( int i = 0, N = parents().size(); i < N; ++i )
if ( parents()[i]->hasColourLine(line, anti) ) return parents()[i];
return tPPtr();
}
tPPtr Particle::outgoingColour(bool anti) const {
if ( !hasColourInfo() ) return tPPtr();
tColinePtr line = colourLine(anti);
if ( !line ) return tPPtr();
for ( int i = 0, N = children().size(); i < N; ++i )
if ( children()[i]->hasColourLine(line, anti) ) return children()[i];
return tPPtr();
}
LorentzPoint Particle::labVertex() const {
LorentzPoint r(rep().theBirthStep && rep().theBirthStep->collision()?
vertex() + rep().theBirthStep->collision()->vertex():
vertex());
return r;
}
void Particle::setLabVertex(const LorentzPoint & p) {
rep().theVertex = ( rep().theBirthStep && rep().theBirthStep->collision()?
p - rep().theBirthStep->collision()->vertex() : p );
}
void Particle::transform(const LorentzRotation & r) {
if ( hasRep() && spinInfo() ) spinInfo()->transform(momentum(), r);
theMomentum.transform(r);
if ( !hasRep() ) return;
rep().theVertex.transform(r);
rep().theLifeLength.transform(r);
}
void Particle::deepTransform(const LorentzRotation & r) {
transform(r);
if ( !theRep ) return;
for ( int i = 0, N = children().size(); i < N; ++i )
rep().theChildren[i]->deepTransform(r);
if ( rep().theNext ) rep().theNext->deepTransform(r);
}
void Particle::rotateX(double a) {
LorentzRotation r;
r.rotateX(a);
transform(r);
}
void Particle::deepRotateX(double a) {
LorentzRotation r;
r.rotateX(a);
deepTransform(r);
}
void Particle::rotateY(double a) {
LorentzRotation r;
r.rotateY(a);
transform(r);
}
void Particle::deepRotateY(double a) {
LorentzRotation r;
r.rotateY(a);
deepTransform(r);
}
void Particle::rotateZ(double a) {
LorentzRotation r;
r.rotateZ(a);
transform(r);
}
void Particle::deepRotateZ(double a) {
LorentzRotation r;
r.rotateZ(a);
deepTransform(r);
}
void Particle::rotate(double a, const Axis & axis) {
LorentzRotation r;
r.rotate(a, axis);
transform(r);
}
void Particle::deepRotate(double a, const Axis & axis) {
LorentzRotation r;
r.rotate(a, axis);
deepTransform(r);
}
string Particle::outputFormat =
"%n3%s10 %i7 %p[,]0 %c(,) %^^0%vv0 %>>0%<>0 %l{,}0\n"
" %x10.3%y10.3%z10.3%e10.3%m10.3\n";
int getNumber(string::const_iterator & pos, int def) {
if ( !isdigit(*pos) ) return def;
def = *pos++ - '0';
while ( isdigit(*pos) ) def = 10*def + *pos++ - '0';
return def;
}
void writePrecision(ostream & os, string::const_iterator & pos,
int defw, int defp, double x) {
defw = getNumber(pos, defw);
if ( *pos == '.' ) defp = getNumber(++pos, defp);
int oldp = os.precision();
os << setprecision(defp) << setw(defw) << x << setprecision(oldp);
}
void writeStringAdjusted(ostream & os, bool left, int w, string str) {
while ( !left && w-- > int(str.size()) ) os << ' ';
os << str;
while ( left && w-- > int(str.size()) ) os << ' ';
}
template <typename Container>
void writeParticleRanges(ostream & os, const Container & co, char sep, int w) {
set<int> cnum;
for ( typename Container::const_iterator it = co.begin();
it != co.end(); ++it) cnum.insert((**it).number());
bool elipsis = false;
int last = -10;
for ( set<int>::iterator it = cnum.begin(); it != cnum.end(); ++it) {
int n = *it;
int next = 0;
set<int>::iterator itn = it;
if ( ++itn != cnum.end() ) next = *itn;
bool writeit = true;
bool writesep = false;
if ( elipsis && ( n != last + 1 || n != next - 1 ) )
elipsis = false;
else if ( !elipsis && n == last + 1 && n == next -1 ) {
os << "..";
elipsis = true;
writeit = false;
}
else if ( elipsis && n == last + 1 && n == next -1 )
writeit = false;
else if ( it != cnum.begin() )
writesep = true;
if ( writeit ) {
if ( writesep ) os << sep;
os << setw(w) << n;
}
last = n;
}
}
ostream & ThePEG::operator<<(ostream & os, const Particle & p) {
return p.print(os, p.birthStep());
}
ostream & Particle::print(ostream & os, tcStepPtr step) const {
if ( !step ) step = birthStep();
tCollPtr coll = step? step->collision(): tCollPtr();
tEventPtr event = coll? coll->event(): tEventPtr();
string::const_iterator pos = Particle::outputFormat.begin();
ios::fmtflags saveflags = os.setf(ios::fixed, ios::floatfield);
while ( pos != Particle::outputFormat.end() ) {
if ( *pos == '%' && ++pos != Particle::outputFormat.end() ) {
bool left = false;
if ( *pos == '-' ) {
left = true;
os.setf(ios::left, ios::adjustfield);
++pos;
} else {
os.setf(ios::right, ios::adjustfield);
}
char mark;
char open;
char close;
char sep;
int w;
string str;
string fill;
if ( pos == Particle::outputFormat.end() ) break;
bool fullColour = false;
switch ( *pos ) {
case 'n':
os << setw(getNumber(++pos, 3)) << number();
break;
case 'i':
os << setw(getNumber(++pos, 8)) << id();
break;
case 's':
writeStringAdjusted(os, left, getNumber(++pos, 8), PDGName());
break;
case 'x':
writePrecision(os, ++pos, 10, 3, momentum().x()/GeV);
break;
case 'y':
writePrecision(os, ++pos, 10, 3, momentum().y()/GeV);
break;
case 'z':
writePrecision(os, ++pos, 10, 3, momentum().z()/GeV);
break;
case 'e':
writePrecision(os, ++pos, 10, 3, momentum().e()/GeV);
break;
case 'm':
writePrecision(os, ++pos, 10, 3, momentum().mass()/GeV);
break;
case 'P':
fullColour = true;
case 'p':
open = *++pos;
sep = *++pos;
close = *++pos;
w = getNumber(++pos, 0);
if ( parents().empty() ) break;
if ( open ) os << open;
writeParticleRanges(os, parents(), sep, w);
if ( fullColour && hasColourInfo() &&
( incomingColour() || incomingAntiColour() ) ) {
if ( close ) os << open;
if ( incomingColour() )
os << "+" << incomingColour()->number();
if ( incomingAntiColour() )
os << "-" << incomingAntiColour()->number();
if ( close ) os << close;
}
if ( close ) os << close;
break;
case 'l':
open = *++pos;
sep = *++pos;
close = *++pos;
w = getNumber(++pos, 0);
if ( hasColourInfo() &&
( colourLine() || antiColourLine() ) && event) {
if ( open ) os << open;
vector<tcColinePtr> clines = colourInfo()->colourLines();
for ( int i = 0, N = clines.size(); i < N; ++i ) {
if ( i > 0 && sep ) os << sep;
clines[i]->write(os, event, false);
}
vector<tcColinePtr> aclines = colourInfo()->antiColourLines();
for ( int i = 0, N = aclines.size(); i < N; ++i ) {
if ( ( i > 0 || clines.size() ) && sep ) os << sep;
aclines[i]->write(os, event, true);
}
if ( close ) os << close;
}
break;
case 'C':
fullColour = true;
case 'c':
open = *++pos;
sep = *++pos;
close = *++pos;
w = getNumber(++pos, 0);
if ( children().empty() ) break;
if ( open ) os << open;
writeParticleRanges(os, children(), sep, w);
if ( fullColour && hasColourInfo() &&
( outgoingColour() || outgoingAntiColour() ) ) {
if ( close ) os << open;
if ( outgoingColour() )
os << "+" << outgoingColour()->number();
if ( outgoingAntiColour() )
os << "-" << outgoingAntiColour()->number();
if ( close ) os << close;
}
if ( close ) os << close;
break;
case '>':
mark = *++pos;
w = getNumber(++pos, 0);
if ( hasColourInfo() && step && step->colourNeighbour(this) ) {
os << setw(w-1) << step->colourNeighbour(this)->number() << mark;
}
break;
case '<':
mark = *++pos;
w = getNumber(++pos, 0);
if ( hasColourInfo() && step && step->antiColourNeighbour(this) ) {
int n = step->antiColourNeighbour(this)->number();
ostringstream oss;
oss << mark << n;
writeStringAdjusted(os, left, w, oss.str());
}
break;
case 'v':
mark = *++pos;
w = getNumber(++pos, 0);
if ( next() ) {
if ( left && mark ) os << mark;
os << setw(w) << next()->number();
if ( !left && mark ) os << mark;
}
break;
case '^':
mark = *++pos;
w = getNumber(++pos, 0);
if ( previous() ) {
if ( left && mark ) os << mark;
os << setw(w) << previous()->number();
if ( !left && mark ) os << mark;
}
break;
case 'd':
switch ( *++pos ) {
case 'x':
writePrecision(os, ++pos, 10, 3, lifeLength().x()/mm);
break;
case 'y':
writePrecision(os, ++pos, 10, 3, lifeLength().y()/mm);
break;
case 'z':
writePrecision(os, ++pos, 10, 3, lifeLength().z()/mm);
break;
case 't':
writePrecision(os, ++pos, 10, 3, lifeLength().e()/mm);
break;
case 'T':
writePrecision(os, ++pos, 10, 3, lifeLength().tau()/mm);
break;
}
break;
case 'V':
switch ( *++pos ) {
case 'x':
writePrecision(os, ++pos, 10, 3, vertex().x()/mm);
break;
case 'y':
writePrecision(os, ++pos, 10, 3, vertex().y()/mm);
break;
case 'z':
writePrecision(os, ++pos, 10, 3, vertex().z()/mm);
break;
case 't':
writePrecision(os, ++pos, 10, 3, vertex().e()/mm);
break;
}
case 'L':
switch ( *++pos ) {
case 'x':
writePrecision(os, ++pos, 10, 3, labVertex().x()/mm);
break;
case 'y':
writePrecision(os, ++pos, 10, 3, labVertex().y()/mm);
break;
case 'z':
writePrecision(os, ++pos, 10, 3, labVertex().z()/mm);
break;
case 't':
writePrecision(os, ++pos, 10, 3, labVertex().e()/mm);
break;
}
break;
default:
os << *pos++;
}
} else {
if ( pos != Particle::outputFormat.end() ) os << *pos++;
}
}
os.flags(saveflags);
return os;
}
void Particle::debugme() const {
cerr << *this;
EventRecordBase::debugme();
}
void Particle::persistentOutput(PersistentOStream & os) const {
EventConfig::putParticleData(os, theData);
os << ounit(theMomentum, GeV) << bool( theRep != 0 );
if ( !theRep ) return;
os << rep().theParents << rep().theChildren
<< rep().thePrevious << rep().theNext << rep().theBirthStep
<< ounit(rep().theVertex, mm) << ounit(rep().theLifeLength, mm)
<< ounit(rep().theScale, GeV2) << ounit(rep().theVetoScale, GeV2)
<< rep().theNumber << rep().theDecayMode
<< rep().theColourInfo << rep().theSpinInfo << rep().theExtraInfo;
}
void Particle::persistentInput(PersistentIStream & is, int) {
bool readRep;
EventConfig::getParticleData(is, theData);
is >> iunit(theMomentum, GeV) >> readRep;
if ( !readRep ) return;
if ( !hasRep() ) theRep = new ParticleRep;
is >> rep().theParents >> rep().theChildren
>> rep().thePrevious >> rep().theNext >> rep().theBirthStep
>> iunit(rep().theVertex, mm) >> iunit(rep().theLifeLength, mm)
>> iunit(rep().theScale, GeV2) >> iunit(rep().theVetoScale, GeV2)
>> rep().theNumber >> rep().theDecayMode
>> rep().theColourInfo >> rep().theSpinInfo >> rep().theExtraInfo;
}
ClassDescription<Particle> Particle::initParticle;
void Particle::Init() {}
-
-ThePEG_IMPLEMENT_SET(PPtr,ParticleSet)
-
diff --git a/EventRecord/Step.cc b/EventRecord/Step.cc
--- a/EventRecord/Step.cc
+++ b/EventRecord/Step.cc
@@ -1,450 +1,448 @@
// -*- C++ -*-
//
// Step.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG 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 Step class.
//
#include "Step.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Config/algorithm.h"
#include <iostream>
#include <iomanip>
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
#include "Step.tcc"
#endif
using namespace ThePEG;
Step::Step(const Step & s)
: Base(s),
theParticles(s.theParticles), theIntermediates(s.theIntermediates),
allParticles(s.allParticles), theCollision(s.theCollision),
theHandler(s.theHandler) {}
Step::~Step() {
for ( ParticleSet::iterator it = allParticles.begin();
it != allParticles.end(); ++it )
if ( (**it).hasRep() && (**it).birthStep() == this )
(**it).rep().theBirthStep = tStepPtr();
theParticles.clear();
theIntermediates.clear();
theSubProcesses.clear();
allParticles.clear();
theCollision = tCollPtr();
theHandler = tcEventBasePtr();
}
StepPtr Step::clone() const {
return ptr_new<StepPtr>(*this);
}
void Step::rebind(const EventTranslationMap & trans) {
ParticleSet::const_iterator pit;
allParticles.clear();
for ( pit = theParticles.begin(); pit != theParticles.end(); ++pit )
allParticles.insert(trans.translate(*pit));
theParticles.swap(allParticles);
allParticles.clear();
for ( pit = theIntermediates.begin(); pit != theIntermediates.end(); ++pit )
allParticles.insert(trans.translate(*pit));
theIntermediates = allParticles;
allParticles.insert(theParticles.begin(), theParticles.end());
for ( pit = allParticles.begin(); pit != allParticles.end(); ++pit )
(**pit).rebind(trans);
}
void Step::addParticle(tPPtr p) {
if ( !p->birthStep() ) p->rep().theBirthStep = this;
theParticles.insert(p);
allParticles.insert(p);
if ( collision() ) collision()->addParticle(p);
}
void Step::addSubProcess(tSubProPtr sp) {
if (( !member(allParticles, sp->incoming().first) && collision() &&
sp->incoming().first != collision()->incoming().first) &&
std::find(sp->incoming().first->parents().begin(),
sp->incoming().first->parents().end(),
collision()->incoming().first)==
sp->incoming().first->parents().end()) {
collision()->incoming().first->
rep().theChildren.push_back(sp->incoming().first);
sp->incoming().first->rep().theParents.push_back(collision()->incoming().first);
}
if (( !member(allParticles, sp->incoming().second) && collision() &&
sp->incoming().second != collision()->incoming().second) &&
std::find(sp->incoming().second->parents().begin(),
sp->incoming().second->parents().end(),collision()->incoming().second)==
sp->incoming().second->parents().end()) {
collision()->incoming().second->
rep().theChildren.push_back(sp->incoming().second);
sp->incoming().second->rep().theParents.push_back(collision()->incoming().second);
}
if ( !sp->incoming().first->birthStep() )
sp->incoming().first->rep().theBirthStep = this;
if ( !sp->incoming().second->birthStep() )
sp->incoming().second->rep().theBirthStep = this;
addIntermediate(sp->incoming().first);
addIntermediate(sp->incoming().second);
addIntermediates(sp->intermediates().begin(), sp->intermediates().end());
addParticles(sp->outgoing().begin(), sp->outgoing().end());
theSubProcesses.push_back(sp);
if ( collision() ) collision()->addSubProcess(sp);
}
void Step::removeSubProcess(tSubProPtr sp) {
SubProcessVector::iterator sit = ThePEG::find(theSubProcesses, sp);
if ( sit == theSubProcesses.end() ) return;
for ( int i = 0, N = sp->outgoing().size(); i < N; ++i )
removeParticle(sp->outgoing()[i]);
for ( int i = 0, N = sp->intermediates().size(); i < N; ++i )
removeParticle(sp->intermediates()[i]);
removeParticle(sp->incoming().first);
removeParticle(sp->incoming().second);
theSubProcesses.erase(sit);
if ( collision() ) collision()->removeSubProcess(sp);
}
void Step::addIntermediate(tPPtr p) {
theIntermediates.insert(p);
ParticleSet::iterator pit = theParticles.find(p);
if ( pit != theParticles.end() ) theParticles.erase(pit);
else {
if ( !p->birthStep() ) p->rep().theBirthStep = this;
allParticles.insert(p);
if ( collision() ) collision()->addParticle(p);
}
}
void Step::
insertIntermediate(tPPtr p, tPPtr parent, tPPtr child) {
if ( !p->birthStep() ) p->rep().theBirthStep = this;
addIntermediate(p);
parent->removeChild(child);
child->removeParent(parent);
if ( parent ) {
parent->rep().theChildren.push_back(p);
p->rep().theParents.push_back(parent);
}
if ( child ) {
p->rep().theChildren.push_back(child);
child->rep().theParents.push_back(p);
}
}
void Step::removeEntry(tPPtr p) {
ParticleSet::iterator it = allParticles.find(p);
if ( it == allParticles.end() ) return;
allParticles.erase(it);
it = theParticles.find(p);
if ( it != theParticles.end() ) theParticles.erase(it);
if ( p->previous() ) {
it = theIntermediates.find(p->previous());
if ( it != theIntermediates.end() ) theIntermediates.erase(it);
theParticles.insert(p->previous());
allParticles.insert(p->previous());
}
while ( !p->parents().empty() ) {
PPtr parent = p->parents().back();
p->removeParent(parent);
parent->removeChild(p);
if ( !parent->children().empty() ) continue;
it = theIntermediates.find(parent);
if ( it != theIntermediates.end() ) theIntermediates.erase(it);
theParticles.insert(parent);
allParticles.insert(parent);
}
if ( p->hasColourInfo() ) {
if ( colourNeighbour(p) )
colourNeighbour(p)->antiColourNeighbour(antiColourNeighbour(p));
if ( antiColourNeighbour(p) )
antiColourNeighbour(p)->colourNeighbour(colourNeighbour(p));
if ( p->incomingColour() ) p->outgoingColour(tPPtr());
if ( p->incomingAntiColour() ) p->outgoingAntiColour(tPPtr());
}
it = theIntermediates.find(p);
if ( it != theIntermediates.end() ) theIntermediates.erase(it);
}
void Step::removeParticle(tPPtr p) {
if ( p->next() ) removeParticle(p->next());
while ( !p->children().empty() ) removeParticle(p->children().back());
removeEntry(p);
}
bool Step::nullStep() const {
for ( ParticleSet::const_iterator it = allParticles.begin();
it != allParticles.end(); ++it )
if ( (**it).birthStep() == this ) return false;
return true;
}
tPPtr Step::copyParticle(tcPPtr pin) {
PPtr cp;
tPPtr p = const_ptr_cast<tPPtr>(pin);
if ( !collision() ) return cp;
ParticleSet::iterator pit = theParticles.find(p);
if ( collision()->finalStep() != this || p->next() ||
! p->children().empty() || pit == theParticles.end() ) return cp;
cp = p->clone();
cp->rep().thePrevious = p;
p->rep().theNext = cp;
if ( p->hasColour() ) p->colourFlow(cp);
if ( p->hasAntiColour() ) p->antiColourFlow(cp);
cp->rep().theBirthStep = this;
theParticles.erase(pit);
if ( p->birthStep() == this ) theIntermediates.insert(p);
addParticle(cp);
return cp;
}
bool Step::setCopy(tcPPtr poldin, tPPtr pnew) {
if ( poldin->id() != pnew->id() ) return false;
tPPtr pold = const_ptr_cast<tPPtr>(poldin);
pold->rep().theNext = pnew;
pnew->rep().thePrevious = pold;
theParticles.erase(pold);
if ( pold->birthStep() == this ) theIntermediates.insert(pold);
pnew->rep().theBirthStep = this;
addParticle(pnew);
return true;
}
tPPtr Step::insertCopy(tcPPtr pin) {
PPtr cp;
tPPtr p = const_ptr_cast<tPPtr>(pin);
if ( !collision() ) return cp;
if ( collision()->all().find(p) == collision()->all().end() ) return cp;
cp = p->clone();
cp->rep().theNext = p;
cp->rep().theChildren.clear();
if ( p->previous() ) {
p->previous()->rep().theNext = cp;
cp->rep().thePrevious = p->previous();
} else {
for ( int i = 0, N = p->parents().size(); i < N; ++i ) {
tPPtr parent = p->parents()[i];
for ( int j = 0, M = parent->children().size(); j < M; ++j )
if ( parent->children()[j] == p ) parent->rep().theChildren[j] = cp;
}
}
p->rep().theParents.clear();
p->rep().thePrevious = cp;
if ( p->hasColour() ) cp->colourFlow(p);
if ( p->hasAntiColour() ) cp->antiColourFlow(p);
cp->rep().theBirthStep = this;
theIntermediates.insert(cp);
return cp;
}
bool Step::removeDecayProduct(tcPPtr par, tPPtr child) {
if ( !collision() ) return false;
tPPtr parent = const_ptr_cast<tPPtr>(par->final());
if ( collision()->all().find(parent) == collision()->all().end() )
return false;
if ( !par->hasRep() ) return false;
PVector::iterator it = ThePEG::find(parent->rep().theChildren, child);
if ( it == parent->rep().theChildren.end() ) return false;
parent->rep().theChildren.erase(it);
ParticleSet::iterator cit = theParticles.find(child);
if ( cit != theParticles.end() ) {
theParticles.erase(cit);
if ( child->birthStep() == this ) theIntermediates.insert(child);
}
return true;
}
bool Step::addDecayProduct(tcPPtr par, tPPtr child, bool fixColour) {
if ( !collision() ) return false;
tPPtr parent = const_ptr_cast<tPPtr>(par->final());
if ( collision()->finalStep() != this || parent->next() ) return false;
ParticleSet::iterator pit = theParticles.find(parent);
if ( pit != theParticles.end() ) {
theParticles.erase(pit);
if ( parent->birthStep() == this ) theIntermediates.insert(parent);
} else {
if ( parent != collision()->incoming().first &&
parent != collision()->incoming().second &&
parent->children().empty() ) return false;
}
parent->rep().theChildren.push_back(child);
child->rep().theParents.push_back(parent);
child->rep().theBirthStep = this;
addParticle(child);
if ( !fixColour || !parent->hasColourInfo() || !parent->coloured() ||
!child->coloured() ) return true;
if ( parent->hasColour() && child->hasColour() &&
!parent->outgoingColour() && !child->colourLine() )
parent->colourFlow(child);
if ( parent->hasAntiColour() && child->hasAntiColour() &&
!child->antiColourLine() ) {
if ( parent->outgoingAntiColour() )
parent->antiColourLine()->
removeAntiColoured(parent->outgoingAntiColour());
parent->antiColourFlow(child);
}
return true;
}
void Step::addDecayNoCheck(tPPtr parent, tPPtr child) {
ParticleSet::iterator pit = theParticles.find(parent);
if ( pit != theParticles.end() ) {
theParticles.erase(pit);
if ( parent->birthStep() == this ) theIntermediates.insert(parent);
}
child->rep().theBirthStep = this;
addParticle(child);
}
void Step::addDecayProduct(tPPtr child) {
for ( int i = 0, N = child->parents().size(); i < N; ++i ) {
ParticleSet::iterator pit = theParticles.find(child->parents()[i]);
if ( pit != theParticles.end() ) {
theParticles.erase(pit);
if ( child->parents()[i]->birthStep() == this )
theIntermediates.insert(child->parents()[i]);
}
}
child->rep().theBirthStep = this;
addParticle(child);
}
void Step::fixColourFlow() {
tParticleVector news;
for ( ParticleSet::iterator pi = theParticles.begin();
pi != theParticles.end(); ++pi )
if ( (**pi).birthStep() == this ) news.push_back(*pi);
for ( int i = 0, N = news.size(); i < N; ++i ) {
tPPtr p = news[i];
if ( p->hasColour() && !antiColourNeighbour(p) ) {
tPPtr ng = p;
while ( ( ng = ng->incomingColour() ) && !antiColourNeighbour(ng) ) {}
if ( ng ) {
ng = antiColourNeighbour(ng);
if ( !ng->outgoingColour() ) ng = copyParticle(ng);
while ( ng->outgoingColour() ) ng = ng->outgoingColour();
p->antiColourConnect(ng);
}
}
if ( p->hasAntiColour() && !colourNeighbour(p) ) {
tPPtr ng = p;
while ( ( ng = ng->incomingAntiColour() ) && !colourNeighbour(ng) ) {}
if ( ng ) {
ng = colourNeighbour(ng);
if ( !ng->outgoingAntiColour() ) ng = copyParticle(ng);
while ( ng->outgoingAntiColour() ) ng = ng->outgoingAntiColour();
p->colourConnect(ng);
}
}
}
}
tPPtr Step::antiColourNeighbour(tcPPtr p) const {
return colourNeighbour(p, true);
}
tPPtr Step::colourNeighbour(tcPPtr p, bool anti) const {
if ( !member(particles(), const_ptr_cast<tPPtr>(p)) ) return tPPtr();
tColinePtr line = p->colourLine(!anti);
if ( !line ) return tPPtr();
for ( ParticleSet::const_iterator it = particles().begin();
it != particles().end(); ++it )
if ( (**it).hasColourLine(line, anti) ) return *it;
return tPPtr();
}
vector<tPVector> Step::getSinglets(tParticleSet & left) {
vector<tPVector> ret;
while ( !left.empty() ) {
tPPtr first = *left.begin();
left.erase(first);
if ( !first->hasColourInfo() || !first->coloured() ) continue;
tPPtr last = first;
tPPtr test;
while ( ( test = last->antiColourNeighbour(left.begin(), left.end()) ) &&
test != first )
last = test;
while ( ( test = first->colourNeighbour(left.begin(), left.end()) ) &&
test != last )
first = test;
ret.push_back(tPVector());
for ( ; first != last;
first = first->antiColourNeighbour(left.begin(), left.end()) ) {
left.erase(first);
ret.back().push_back(first);
}
left.erase(first);
ret.back().push_back(first);
}
return ret;
}
ostream & ThePEG::operator<<(ostream & os, const Step & s) {
if ( !s.intermediates().empty() ) os << "--- intermediates:" << endl;
Particle::PrintParticles(os, s.intermediates(), &s);
os << "--- final:" << endl;
LorentzMomentum sum;
Energy2 sumx = Energy2();
Energy2 sumy = Energy2();
Energy2 sumz = Energy2();
Particle::PrintParticles(os, s.particles(), &s);
for ( ParticleSet::const_iterator it = s.particles().begin();
it != s.particles().end(); ++it ) {
sum += (**it).momentum();
sumx += sqr((**it).momentum().x());
sumy += sqr((**it).momentum().y());
sumz += sqr((**it).momentum().z());
}
os << string(78, '-') << endl << " Sum of momenta: ";
int oldprecision = os.precision();
Energy sumx1
= ( sqr(sum.x()) > Constants::epsilon*sumx ?
sum.x(): ZERO );
Energy sumy1
= ( sqr(sum.y()) > Constants::epsilon*sumy ?
sum.y(): ZERO );
Energy sumz1
= ( sqr(sum.z()) > Constants::epsilon*sumz ?
sum.z(): ZERO );
os << setprecision(3) << setw(10) << sumx1/GeV << setw(10) << sumy1/GeV
<< setw(10) << sumz1/GeV << setw(10) << sum.e()/GeV
<< setw(10) << sum.m()/GeV << endl << setprecision(oldprecision);
return os;
}
void Step::debugme() const {
cerr << *this;
EventRecordBase::debugme();
}
void Step::persistentOutput(PersistentOStream & os) const {
os << theParticles << theIntermediates << theSubProcesses << allParticles
<< theCollision;
EventConfig::putHandler(os, theHandler);
}
void Step::persistentInput(PersistentIStream & is, int) {
is >> theParticles >> theIntermediates >> theSubProcesses >> allParticles
>> theCollision;
EventConfig::getHandler(is, theHandler);
}
ClassDescription<Step> Step::initStep;
void Step::Init() {}
-
-ThePEG_IMPLEMENT_SET(StepPtr,StepSet)
diff --git a/EventRecord/SubProcess.cc b/EventRecord/SubProcess.cc
--- a/EventRecord/SubProcess.cc
+++ b/EventRecord/SubProcess.cc
@@ -1,152 +1,150 @@
// -*- C++ -*-
//
// SubProcess.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
// Copyright (C) 2009-2011 Simon Platzer
//
// ThePEG 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 SubProcess class.
//
#include "SubProcess.h"
#include "ThePEG/EventRecord/Collision.h"
#include "ThePEG/Config/algorithm.h"
#include "ThePEG/EventRecord/ParticleTraits.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include <iostream>
#ifdef ThePEG_TEMPLATES_IN_CC_FILE
#include "SubProcess.tcc"
#endif
using namespace ThePEG;
SubProcess::
SubProcess(const PPair & newIncoming,
tCollPtr newCollision, tcEventBasePtr newHandler,
tSubProPtr newHead, double newGroupWeight)
: theHandler(newHandler), theCollision(newCollision),
theIncoming(newIncoming), isDecayed(false),
theHead(newHead), theGroupWeight(newGroupWeight) {}
SubProcess::~SubProcess() {}
SubProPtr SubProcess::clone() const {
return ptr_new<SubProPtr>(*this);
}
void SubProcess::addIntermediate(tPPtr p, bool fixrelations) {
if ( fixrelations ) {
incoming().first->rep().theChildren.push_back(p);
incoming().second->rep().theChildren.push_back(p);
p->rep().theParents.push_back(incoming().first);
p->rep().theParents.push_back(incoming().second);
}
theIntermediates.push_back(p);
}
void SubProcess::addOutgoing(tPPtr p, bool fixrelations) {
if ( fixrelations ) {
if ( intermediates().empty() ) {
incoming().first->rep().theChildren.push_back(p);
incoming().second->rep().theChildren.push_back(p);
p->rep().theParents.push_back(incoming().first);
p->rep().theParents.push_back(incoming().second);
} else {
for (ParticleVector::iterator it = theIntermediates.begin();
it != theIntermediates.end(); ++it ) {
(**it).rep().theChildren.push_back(p);
p->rep().theParents.push_back(*it);
}
}
}
theOutgoing.push_back(p);
}
void SubProcess::changeIncoming(tPPtr pnew, tPPtr pold) {
if(pold==theIncoming.first) {
theIntermediates.push_back(pold);
theIncoming.first = pnew;
}
else if(pold==theIncoming.second) {
theIntermediates.push_back(pold);
theIncoming.second = pnew;
}
}
void SubProcess::rebind(const EventTranslationMap & trans) {
theIncoming.first = trans.translate(theIncoming.first);
theIncoming.second = trans.translate(theIncoming.second);
theCollision = trans.translate(theCollision);
for ( ParticleVector::iterator pit = theOutgoing.begin();
pit != theOutgoing.end(); ++pit )
*pit = trans.translate(*pit);
for ( ParticleVector::iterator pit = theIntermediates.begin();
pit != theIntermediates.end(); ++pit )
*pit = trans.translate(*pit);
}
void SubProcess::removeEntry(tPPtr p) {
if ( p == theIncoming.first ) theIncoming.first = PPtr();
if ( p == theIncoming.second ) theIncoming.second = PPtr();
ParticleVector::iterator pit = theOutgoing.begin();
while ( pit != theOutgoing.end() ) {
if ( *pit == p ) pit = theOutgoing.erase(pit);
else ++pit;
}
pit = theIntermediates.begin();
while ( pit != theIntermediates.end() ) {
if ( *pit == p ) pit = theIntermediates.erase(pit);
else ++pit;
}
}
void SubProcess::transform(const LorentzRotation & r) {
incoming().first->transform(r);
incoming().second->transform(r);
for_each(intermediates(), Transformer(r));
for_each(outgoing(), Transformer(r));
}
void SubProcess::printMe(ostream& os) const {
os << "--- incoming:" << endl
<< *incoming().first << *incoming().second;
if ( !intermediates().empty() ) os << "--- intermediates:" << endl;
Particle::PrintParticles(os, intermediates().begin(),
intermediates().end());
os << "--- outgoing:" << endl;
Particle::PrintParticles(os, outgoing().begin(), outgoing().end());
}
ostream & ThePEG::operator<<(ostream & os, const SubProcess & sp) {
sp.printMe(os);
return os;
}
void SubProcess::debugme() const {
cerr << *this;
EventRecordBase::debugme();
}
void SubProcess::persistentOutput(PersistentOStream & os) const {
EventConfig::putHandler(os, theHandler);
os << theCollision << theIncoming << theIntermediates << theOutgoing
<< isDecayed << theHead << theGroupWeight;
}
void SubProcess::persistentInput(PersistentIStream & is, int) {
EventConfig::getHandler(is, theHandler);
is >> theCollision >> theIncoming >> theIntermediates >> theOutgoing
>> isDecayed >> theHead >> theGroupWeight;
}
ClassDescription<SubProcess> SubProcess::initSubProcess;
void SubProcess::Init() {}
-
-ThePEG_IMPLEMENT_SET(SubProPtr,SubProcessSet)
diff --git a/PDT/DecayMode.cc b/PDT/DecayMode.cc
--- a/PDT/DecayMode.cc
+++ b/PDT/DecayMode.cc
@@ -1,685 +1,680 @@
// -*- C++ -*-
//
// DecayMode.cc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG 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 DecayMode class.
//
#include "DecayMode.h"
#include "DecayMode.xh"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/Rebinder.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Utilities/EnumIO.h"
using namespace ThePEG;
DecayMode::DecayMode()
: theBrat(0.0), isOn(false) {}
DecayMode::DecayMode(tPDPtr newParticle, double newBrat, bool newOn)
: theBrat(newBrat), isOn(newOn), theParent(newParticle) {}
DecayMode::DecayMode(const DecayMode & dm)
: Interfaced(dm), theTag(dm.theTag), theBrat(dm.theBrat), isOn(dm.isOn),
theParent(dm.theParent), theProducts(dm.theProducts),
theOrderedProducts(dm.theOrderedProducts),
theCascadeProducts(dm.theCascadeProducts), theMatchers(dm.theMatchers),
theWildMatcher(dm.theWildMatcher), theExcluded(dm.theExcluded),
theOverlap(dm.theOverlap), theDecayer(dm.theDecayer),
theAntiPartner(dm.theAntiPartner), theLinks(dm.theLinks) {}
DecayMode::~DecayMode() {}
DMPtr DecayMode::Create(tPDPtr newParent, double newBrat, bool newOn) {
DMPtr dm = new_ptr(DecayMode(newParent, newBrat, newOn));
Repository::Register(dm, newParent->fullName() + "/NEWMODE");
if ( !newParent->CC() ) return dm;
DMPtr adm = new_ptr(DecayMode(newParent->CC(), newBrat, newOn));
Repository::Register(adm, newParent->CC()->fullName() + "/NEWMODE");
dm->theAntiPartner = adm;
adm->theAntiPartner = dm;
return dm;
}
void DecayMode::readSetup(istream & is) {
string decnam;
is >> theBrat >> ienum(isOn) >> decnam;
if ( decnam.empty() ) return;
BaseRepository::DirectoryAppend(decnam);
setDecayer(BaseRepository::GetObject<tDecayerPtr>(decnam));
}
IBPtr DecayMode::clone() const {
return dmclone();
}
DMPtr DecayMode::dmclone() const {
return new_ptr(*this);
}
IBPtr DecayMode::fullclone() const {
DMPtr dm = dmclone();
Repository::Register(dm);
if ( !CC() ) return dm;
DMPtr adm = CC()->dmclone();
Repository::Register(adm);
dm->theAntiPartner = adm;
adm->theAntiPartner = dm;
return dm;
}
void DecayMode::doupdate() {
Interfaced::doupdate();
bool redo = touched();
UpdateChecker::check(decayer(), redo);
if ( !redo ) return;
if ( theBrat > 0.0 && isOn && !decayer()->accept(*this) )
throw DecModNoAccept(tag(), decayer()->name());
}
void DecayMode::rebind(const TranslationMap & trans) {
try {
theParent = trans.alwaysTranslate(theParent);
ParticleMSet newProducts;
trans.alwaysTranslate(inserter(newProducts),
products().begin(), products().end());
products().swap(newProducts);
tPDVector newOrdered;
trans.alwaysTranslate(inserter(newOrdered),
orderedProducts().begin(), orderedProducts().end());
theOrderedProducts.swap(newOrdered);
ModeMSet newCasc;
trans.alwaysTranslate(inserter(newCasc),
cascadeProducts().begin(), cascadeProducts().end());
cascadeProducts().swap(newCasc);
MatcherMSet newMatchers;
trans.alwaysTranslate(inserter(newMatchers),
productMatchers().begin(), productMatchers().end());
productMatchers().swap(newMatchers);
wildProductMatcher() = trans.alwaysTranslate(wildProductMatcher());
ParticleMSet newExclude;
trans.alwaysTranslate(inserter(newExclude),
excluded().begin(), excluded().end());
excluded().swap(newExclude);
theAntiPartner = trans.alwaysTranslate(CC());
for ( int i = 0, N = theLinks.size(); i < N; ++i ) {
theLinks[i].first = trans.alwaysTranslate(theLinks[i].first);
theLinks[i].second = trans.alwaysTranslate(theLinks[i].second);
}
}
catch (IBPtr ip) {
throw DecModRebind(name(), ip->name());
}
catch (...) {
throw DecModRebind(name(), "<unknown>");
}
}
IVector DecayMode::getReferences() {
IVector ret;
ret.push_back(theParent);
ret.insert(ret.end(), products().begin(), products().end());
ret.insert(ret.end(), cascadeProducts().begin(), cascadeProducts().end());
ret.insert(ret.end(), productMatchers().begin(), productMatchers().end());
if ( wildProductMatcher() ) ret.push_back(wildProductMatcher());
ret.insert(ret.end(), excluded().begin(), excluded().end());
if ( CC() ) ret.push_back(CC());
return ret;
}
bool DecayMode::addOverlap(tcDMPtr d) {
bool inc = includes(*d);
if ( !inc ) return false;
if ( find(theOverlap.begin(), theOverlap.end(), d) != theOverlap.end() )
return true;
theOverlap.push_back(d);
return true;
}
void DecayMode::resetOverlap() {
theOverlap.clear();
}
bool DecayMode::
compareId(const ParticleMSet & s1, const ParticleMSet & si2) const {
if ( generator() ) return s1 == si2;
ParticleMSet s2 = si2;
for ( ParticleMSet::const_iterator p1 = s1.begin(); p1 != s1.end(); ++p1 ) {
ParticleMSet::const_iterator p2 = s2.begin();
while ( p2 != s2.end() && (**p2).id() != (**p1).id() ) ++p2;
if ( p2 == s2.end() ) return false;
s2.erase(p2);
}
return s2.empty();
}
ParticleMSet::const_iterator
DecayMode::findId(const ParticleMSet & s, const ParticleData & p) const {
for ( ParticleMSet::const_iterator pit = s.begin(); pit != s.end(); ++pit )
if ( (**pit).id() == p.id() ) return pit;
return s.end();
}
struct IdCmp {
bool operator()(tcPDPtr p1, tcPDPtr p2) { return p1->id() == p2->id(); }
};
bool DecayMode::includes(const DecayMode & d) const {
// Fast check for ordinary decay modes.
if ( cascadeProducts().empty() && productMatchers().empty() &&
excluded().empty() && !wildProductMatcher() &&
d.cascadeProducts().empty() && d.productMatchers().empty() &&
d.excluded().empty() && !d.wildProductMatcher() ) {
if ( links().size() != d.links().size() ) return false;
if ( !compareId(products(), d.products()) ) return false;
LinkVector dlinks = d.links();
for ( int i = 0, N = links().size(); i < N; ++i ) {
for ( int j = 0, M = dlinks.size(); j < M; ++j ) {
if ( ( links()[i].first->id() == dlinks[j].first->id() &&
links()[i].second->id() == dlinks[j].second->id() ) ||
( links()[i].first->id() == dlinks[j].second->id() &&
links()[i].second->id() == dlinks[j].first->id() ) ) {
dlinks.erase(dlinks.begin() + j);
break;
}
}
return false;
}
return dlinks.empty();
}
// First check that none of the excluded products in this are
// present in the other.
ParticleMSet::const_iterator pit;
for ( pit = excluded().begin(); pit != excluded().end(); ++pit ) {
if ( findId(d.products(), **pit ) != d.products().end() ) return false;
for ( ModeMSet::const_iterator mit = d.cascadeProducts().begin();
mit != d.cascadeProducts().end(); ++mit )
if ( (**pit).id() == (**mit).parent()->id() ) return false;
}
// Check that all cascade decays in this overlaps with one in the
// other. Save the ones that are left
ModeMSet cascleft = d.cascadeProducts();
for ( ModeMSet::iterator mit = cascadeProducts().begin();
mit != cascadeProducts().end(); ++mit ) {
ModeMSet::iterator mit2 = cascleft.begin();
while ( mit2 != cascleft.end() && !(**mit).includes(**mit2) ) ++mit2;
if ( mit2 == cascleft.end() ) return false;
}
// Check that all cascade product parents in the other matches
// something in this. Otherwise expand the cascade product.
ParticleMSet partleft = d.products();
MatcherMSet matchleft = d.productMatchers();
ParticleMSet excludeleft = d.excluded();
MatcherMSet wildleft;
if ( d.wildProductMatcher() ) wildleft.insert(wildProductMatcher());
ParticleMSet part = products();
MatcherMSet match = productMatchers();
while ( cascleft.size() ) {
ModeMSet::iterator cdmit = cascleft.begin();
cDMPtr cdm = *cdmit;
ParticleMSet::iterator pit = findId(part, *(cdm->parent()));
if ( pit != part.end() ) {
cascleft.erase(cdmit);
part.erase(pit);
} else {
MatcherMSet::iterator mit = match.begin();
while ( mit != match.end() && !(**mit).matches(*(cdm->parent())) ) ++mit;
if ( mit != match.end() ) {
cascleft.erase(cdmit);
match.erase(mit);
} else {
if ( wildProductMatcher() &&
wildProductMatcher()->matches(*(cdm->parent())) ) {
cascleft.erase(cdmit);
} else {
cascleft.erase(cdmit);
partleft.insert(cdm->products().begin(), cdm->products().end());
matchleft.insert(cdm->productMatchers().begin(),
cdm->productMatchers().end());
if ( cdm->wildProductMatcher() )
wildleft.insert(cdm->wildProductMatcher());
excludeleft.insert(cdm->excluded().begin(), cdm->excluded().end());
cascleft.insert(cdm->cascadeProducts().begin(),
cdm->cascadeProducts().end());
}
}
}
}
// Check that all excluded left in the other are absent in this.
if ( find_first_of(excludeleft.begin(), excludeleft.end(), part.begin(),
part.end(), IdCmp()) != excludeleft.end() ) return false;
// Now all particles and matches left in this must match something
// in the other.
pit = part.begin();
while ( pit != part.end() ) {
ParticleMSet::iterator pit2 = findId(partleft, **pit++);
if ( pit2 == partleft.end() ) return false;
partleft.erase(pit2);
}
MatcherMSet::const_iterator pmit = match.begin();
while ( pmit != match.end() ) {
ParticleMSet::iterator pit2 = partleft.begin();
while ( pit2 != partleft.end() && ! (**pmit).matches(**pit2) ) ++pit2;
if ( pit2 != partleft.end() ) {
partleft.erase(pit2);
} else {
MatcherMSet::iterator pmit2 = matchleft.begin();
while ( pmit2 != matchleft.end() && ! (**pmit).matches(**pmit2) ) ++pmit2;
if ( pmit2 != matchleft.end() ) {
matchleft.erase(pmit2);
} else
return false;
}
}
// Now all particles and matchers left in the other must be matched
// by the wild match in this.
if ( wildProductMatcher() ) {
pit = partleft.begin();
while ( pit != partleft.end() )
if ( !(wildProductMatcher()->matches(**pit++)) ) return false;
pmit = matchleft.begin();
while (pmit != matchleft.end() )
if ( !(wildProductMatcher()->matches(**pmit++)) ) return false;
pmit = wildleft.begin();
while (pmit != wildleft.end() )
if ( !(wildProductMatcher()->matches(**pmit++)) ) return false;
} else
return partleft.empty() && matchleft.empty() && wildleft.empty();
return true;
}
DMPtr DecayMode::clone(tPDPtr pd) const {
DMPtr dm = dmclone();
dm->theParent = pd;
Repository::Register(dm, pd->fullName() + "/" + dm->name());
if ( !theDecayer || !theDecayer->accept(*dm) ) dm->isOn = false;
if ( pd->CC() ) {
DMPtr adm = CC()? CC()->dmclone(): dmclone();
adm->theParent = pd->CC();
Repository::Register(adm, pd->CC()->fullName() + "/" + adm->name());
dm->theAntiPartner = adm;
adm->theAntiPartner = dm;
if ( !adm->theDecayer->accept(*adm) ) adm->isOn = false;
} else
dm->theAntiPartner = DMPtr();
return dm;
}
void DecayMode::synchronize() {
if ( !CC() ) return;
theBrat = CC()->theBrat;
isOn = CC()->isOn;
theDecayer = CC()->theDecayer;
}
struct ParticleOrdering {
bool operator()(tcPDPtr p1, tcPDPtr p2) {
return abs(p1->id()) > abs(p2->id()) ||
( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
}
};
struct ModeOrdering {
bool operator()(tcDMPtr d1, tcDMPtr d2) {
ParticleOrdering ord;
return ord(d1->parent(), d2->parent()) ||
( !ord(d2->parent(), d1->parent()) &&
( d1->tag() < d2->tag() ||
( d1->tag() == d2->tag() && d1->fullName() < d2->fullName() ) ) );
}
};
struct MatcherOrdering {
bool operator()(tcPMPtr m1, tcPMPtr m2) {
return m1->name() < m2->name() ||
( m1->name() == m2->name() && m1->fullName() < m2->fullName() );
}
};
PVector DecayMode::produceProducts() const {
PVector ret;
for ( int i = 0, N = orderedProducts().size(); i < N; ++i )
ret.push_back(orderedProducts()[i]->produceParticle());
return ret;
}
string DecayMode::makeTag() const {
string ret;
typedef multiset<tcPDPtr,ParticleOrdering> OrderedParticles;
typedef multiset<tcPMPtr,MatcherOrdering> OrderedMatchers;
typedef multiset<tcDMPtr,ModeOrdering> OrderedModes;
LinkVector dlinks = links();
ret = theParent->PDGName() + "->";
if ( dlinks.empty() ) {
OrderedParticles prod(products().begin(), products().end());
for (OrderedParticles::iterator pit = prod.begin();
pit != prod.end(); ++pit )
ret += (**pit).PDGName() + ",";
} else {
unsigned int dl = 0;
for ( int i = 0, N = orderedProducts().size(); i < N; ++i ) {
if ( dl < dlinks.size() && orderedProducts()[i] == dlinks[dl].first ) {
ret += orderedProducts()[i]->PDGName() + "=";
++dl;
} else
ret += orderedProducts()[i]->PDGName() + ",";
}
}
OrderedModes casc(cascadeProducts().begin(), cascadeProducts().end());
for ( OrderedModes::iterator dmit = casc.begin();dmit != casc.end(); ++dmit )
ret += "[" + (**dmit).tag() + "],";
OrderedMatchers match(productMatchers().begin(), productMatchers().end());
for ( OrderedMatchers::iterator mit = match.begin();
mit != match.end(); ++mit ) ret += "?" +(**mit).name() + ",";
if ( theWildMatcher )
ret += "*" + theWildMatcher->name() + ",";
OrderedParticles ex(excluded().begin(), excluded().end());
for ( OrderedParticles::iterator pit = ex.begin(); pit != ex.end(); ++pit )
ret += "!" + (**pit).PDGName() + ",";
ret[ret.size()-1] = ';';
return ret;
}
void DecayMode::brat(double newBrat) {
theBrat = newBrat;
if ( theBrat <= 0.0 ) switchOff();
if ( CC() && parent()->synchronized() ) CC()->theBrat = newBrat;
}
double DecayMode::brat() const {
return isOn? theDecayer->brat(*this, *theParent, theBrat): 0.0;
}
double DecayMode::brat(const Particle & p) const {
return isOn && p.dataPtr() == parent()?
theDecayer->brat(*this, p, theBrat): 0.0;
}
void DecayMode::switchOn() {
isOn = true;
if ( CC() && parent()->synchronized() ) CC()->isOn = true;
}
void DecayMode::switchOff() {
isOn = false;
if ( CC() && parent()->synchronized() ) CC()->isOn = false;
}
void DecayMode::addProduct(tPDPtr pd) {
products().insert(pd);
theOrderedProducts.push_back(pd);
if ( CC() ) {
CC()->products().insert(pd->CC()? pd->CC(): pd);
CC()->theOrderedProducts.push_back(pd->CC()? pd->CC(): pd);
}
resetTag();
}
void DecayMode::addLink(tPDPtr a, tPDPtr b) {
theLinks.push_back(make_pair(a, b));
if ( CC() ) CC()->theLinks.push_back(make_pair(a->CC()? a->CC(): a,
b->CC()? b->CC(): b));
resetTag();
}
void DecayMode::addCascadeProduct(tDMPtr dm) {
cascadeProducts().insert(dm);
if ( CC() ) CC()->cascadeProducts().insert(dm->CC()? dm->CC(): dm);
resetTag();
}
void DecayMode::addProductMatcher( tPMPtr pm) {
productMatchers().insert(pm);
if ( CC() ) CC()->productMatchers().insert(pm->CC()? pm->CC(): pm);
resetTag();
}
void DecayMode::setWildMatcher(tPMPtr pm) {
wildProductMatcher() = pm;
if ( CC() ) CC()->wildProductMatcher() = pm->CC()? pm->CC(): pm;
resetTag();
}
void DecayMode::addExcluded(tPDPtr pd) {
excluded().insert(pd);
if ( CC() ) CC()->excluded().insert(pd->CC()? pd->CC(): pd);
resetTag();
}
void DecayMode::decayer(tDecayerPtr dec) {
if ( !dec || !dec->accept(*this) )
throw DecModSetupNoAccept(tag(), dec->name());
if ( CC() && parent()->synchronized() ) {
if ( !dec->accept(*CC()) )
throw DecModSetupNoAccept(CC()->tag(), dec->name());
CC()->theDecayer = dec;
}
theDecayer = dec;
}
DMPtr DecayMode::constructDecayMode(string & tag, vector<DMPtr> * save) {
DMPtr rdm;
DMPtr adm;
int level = 0;
string::size_type end = 0;
while ( end < tag.size() && ( tag[end] != ']' || level ) ) {
switch ( tag[end++] ) {
case '[':
++level;
break;
case ']':
--level;
break;
}
}
string::size_type next = tag.find("->");
if ( next == string::npos ) return rdm;
if ( tag.find(';') == string::npos ) return rdm;
tPDPtr pd = Repository::findParticle(tag.substr(0,next));
if ( !pd ) return rdm;
rdm = ptr_new<DMPtr>();
rdm->parent(pd);
if ( pd->CC() ) {
adm = ptr_new<DMPtr>();
adm->parent(pd->CC());
rdm->theAntiPartner = adm;
adm->theAntiPartner = rdm;
}
bool error = false;
tag = tag.substr(next+2);
tPDPtr lastprod;
bool dolink = false;
do {
switch ( tag[0] ) {
case '[':
{
tag = tag.substr(1);
DMPtr cdm = constructDecayMode(tag, save);
if ( save ) save->push_back(cdm);
if ( cdm ) rdm->addCascadeProduct(cdm);
else error = true;
} break;
case '=':
dolink = true;
case ',':
case ']':
tag = tag.substr(1);
break;
case '?':
{
next = min(tag.find(','), tag.find(';'));
tPMPtr pm = Repository::findMatcher(tag.substr(1,next-1));
if ( pm ) rdm->addProductMatcher(pm);
else error = true;
tag = tag.substr(next);
} break;
case '!':
{
next = min(tag.find(','), tag.find(';'));
tPDPtr pd = Repository::findParticle(tag.substr(1,next-1));
if ( pd ) rdm->addExcluded(pd);
else error = true;
tag = tag.substr(next);
} break;
case '*':
{
next = min(tag.find(','), tag.find(';'));
tPMPtr pm = Repository::findMatcher(tag.substr(1,next-1));
if ( pm ) rdm->setWildMatcher(pm);
else error = true;
tag = tag.substr(next);
} break;
default:
{
next = min(tag.find('='), min(tag.find(','), tag.find(';')));
tPDPtr pdp = Repository::findParticle(tag.substr(0,next));
if ( pdp ) rdm->addProduct(pdp);
else error = true;
tag = tag.substr(next);
if ( dolink && lastprod ) {
rdm->addLink(lastprod, pdp);
dolink = false;
}
lastprod = pdp;
} break;
}
} while ( tag[0] != ';' && tag.size() );
if ( tag[0] != ';' || error ) {
return DMPtr();
}
tag = tag.substr(1);
for ( DecaySet::const_iterator dit = pd->decayModes().begin();
dit != pd->decayModes().end(); ++dit )
if ( (**dit).tag() == rdm->tag() ) return *dit;
if ( save ) {
save->push_back(rdm);
save->push_back(adm);
} else {
pd->addDecayMode(rdm);
Repository::Register(rdm, pd->fullName() + "/" + rdm->tag());
if ( adm )
Repository::Register(adm, pd->CC()->fullName() + "/" + adm->tag());
}
return rdm;
}
void DecayMode::persistentOutput(PersistentOStream & os) const {
multiset<tcPDPtr,ParticleOrdering> prod(products().begin(), products().end());
multiset<tcDMPtr,ModeOrdering>
casc(cascadeProducts().begin(), cascadeProducts().end());
multiset<tcPMPtr,MatcherOrdering>
match(productMatchers().begin(), productMatchers().end());
multiset<tcPDPtr,ParticleOrdering> ex(excluded().begin(), excluded().end());
multiset<tcDMPtr,ModeOrdering> ovlap(overlap().begin(), overlap().end());
os << theTag << theBrat << isOn << theParent << prod << theOrderedProducts
<< casc << match << theWildMatcher << ex << ovlap << theDecayer
<< theAntiPartner << theLinks;
}
void DecayMode::persistentInput(PersistentIStream & is, int) {
is >> theTag >> theBrat >> isOn >> theParent >> theProducts
>> theOrderedProducts >> theCascadeProducts >> theMatchers
>> theWildMatcher >> theExcluded >> theOverlap >> theDecayer
>> theAntiPartner >> theLinks;
}
ClassDescription<DecayMode> DecayMode::initDecayMode;
void DecayMode::setOn(long i) {
isOn = i;
}
long DecayMode::getOn() const {
return isOn;
}
void DecayMode::setDecayer(DecayerPtr dp) {
decayer(dp);
}
void DecayMode::Init() {
static ClassDocumentation<DecayMode> documentation
("Represents a specific decay channel of a particle.");
static Parameter<DecayMode,double> interfaceBrat
("BranchingRatio",
"The branching fraction for this decay mode. Note that if the sum of "
"branching ratios for one particle is always renormalized to 1. Also, "
"the decaying particle may change this branching ratio if it has a "
"ThePEG::WidthGenerator object assigned to it. ",
&DecayMode::theBrat, 0.0, 0.0, 1.0, false, false, true);
interfaceBrat.setHasDefault(false);
static Switch<DecayMode> interfaceOn
("OnOff",
"Indicates if the decay mode is switched on or off.",
0, 0, false, false, &DecayMode::setOn, &DecayMode::getOn);
static SwitchOption interfaceOnYes
(interfaceOn, "On", "The decay channel is switched on.", 1);
static SwitchOption interfaceOnNo
(interfaceOn, "Off", "The decay channel is switched off.", 0);
interfaceOn.setHasDefault(false);
static Reference<DecayMode,Decayer> interfaceDecayer
("Decayer",
"The ThePEG::Decayer object responsible for performing this decay.",
&DecayMode::theDecayer, false, false, true, false,
&DecayMode::setDecayer);
interfaceBrat.rank(10);
interfaceDecayer.rank(9);
interfaceOn.rank(8);
}
DecModNoAccept::DecModNoAccept(string tag, string dec) {
theMessage << "The Decayer '" << dec << "' is not capable to "
<< "perform the decay in the DecayMode '" << tag << "'.";
severity(warning);
}
DecModSetupNoAccept::DecModSetupNoAccept(string tag, string dec) {
theMessage << "The Decayer '" << dec << "' is not capable to "
<< "perform the decay in the DecayMode '" << tag << "'.";
severity(warning);
}
DecModRebind::DecModRebind(string tag, string obj) {
theMessage << "'Rebind' of DecayMode '" << tag << "' failed because "
<< "the object '" << obj << "' refered to lacked a translation.";
severity(abortnow);
}
-
-ThePEG_IMPLEMENT_MULTISET(tPDPtr,ParticleMSet)
-ThePEG_IMPLEMENT_MULTISET(tPMPtr,MatcherMSet)
-ThePEG_IMPLEMENT_MULTISET(tDMPtr,ModeMSet)
-