diff --git a/Analysis/FactoryBase.cc b/Analysis/FactoryBase.cc --- a/Analysis/FactoryBase.cc +++ b/Analysis/FactoryBase.cc @@ -1,195 +1,195 @@ // -*- C++ -*- // // FactoryBase.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 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 FactoryBase class. // #include "FactoryBase.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Config/algorithm.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "AIDA_helper.h" using namespace ThePEG; FactoryBase::FactoryBase() : theFilename(""), theSuffix("aida"), theStoreType("xml"), theAnalysisFactory(0), theTree(0), theHistogramFactory(0), theDataSetFactory(0) {} FactoryBase::FactoryBase(const FactoryBase & x) : Interfaced(x), theFilename(x.theFilename), theSuffix(x.theSuffix), theStoreType(x.theStoreType), theAnalysisFactory(0), theTree(0), theHistogramFactory(0), theDataSetFactory(0) {} FactoryBase::~FactoryBase() {} FactoryBase::DataFiller::~DataFiller() { int N = v.size()/(3*dset->dimension()); for ( int i = 0; i < N; ++i ) { AIDA::IDataPoint * p = dset->addPoint(); for ( int j = 0; j < p->dimension(); ++j ) { p->coordinate(j)->setValue(v.front()); v.pop_front(); p->coordinate(j)->setErrorPlus(v.front()); v.pop_front(); p->coordinate(j)->setErrorMinus(v.front()); v.pop_front(); } } } void FactoryBase::clear() { if ( theTree ) delete theTree; if ( theAnalysisFactory ) delete theAnalysisFactory; theHistogramFactory = 0; theDataSetFactory = 0; theTree = 0; theAnalysisFactory = 0; } void FactoryBase::dofinish() { Interfaced::dofinish(); - for_each(clients, mem_fun(&InterfacedBase::finish)); + for_each(clients, mem_fn(&InterfacedBase::finish)); tree().commit(); clear(); } void FactoryBase::doinitrun() { Interfaced::doinitrun(); string file = filename(); if ( file == "" ) file = generator()->filename(); else if ( file[0] != '/' ) file = generator()->path() + "/" + file; file += "." + suffix(); theTree = analysisFactory().createTreeFactory()->create (file, storeType(), false, true); theTree->setOverwrite(false); theHistogramFactory = analysisFactory().createHistogramFactory(tree()); theDataSetFactory = analysisFactory().createDataPointSetFactory(tree()); } void FactoryBase::persistentOutput(PersistentOStream & os) const { os << theFilename << theSuffix << theStoreType; } void FactoryBase::persistentInput(PersistentIStream & is, int) { clear(); is >> theFilename >> theSuffix >> theStoreType; } AbstractClassDescription FactoryBase::initFactoryBase; // Definition of the static class description member. void FactoryBase::Init() { static ClassDocumentation documentation ("There is no documentation for the FactoryBase class"); static Parameter interfaceFilename ("Filename", "Together with Suffix, the name of the file " "where the resulting histograms will be stored. If empty, the run-name " "provided by the current EventGenerator will be used instead.", &FactoryBase::theFilename, "", true, false); static Parameter interfaceSuffix ("Suffix", "Together with Filename, the name of the file " "where the resulting histograms will be stored.", &FactoryBase::theSuffix, "aida", true, false); static Parameter interfaceStoreType ("StoreType", "The format in which the histograms are stored in the output file. " "The allowed values depend on the actual AIDA implementation used.", &FactoryBase::theStoreType, "xml", true, false); } AIDA::ITree & FactoryBase::tree() const { return *theTree; } AIDA::IHistogramFactory & FactoryBase::histogramFactory() const { return *theHistogramFactory; } AIDA::IDataPointSetFactory & FactoryBase::dataSetFactory() const { return *theDataSetFactory; } void FactoryBase::mkdir(const string & path) { tree().mkdir(path); } void FactoryBase::mkdirs(const string & path) { tree().mkdirs(path); } void FactoryBase::cd(const string & path) { tree().cd(path); } FactoryBase::tH1DPtr FactoryBase::createHistogram1D(const string & path, int nb, double lo, double up) { return histogramFactory().createHistogram1D(path, nb, lo, up); } FactoryBase::tH1DPtr FactoryBase::createHistogram1D(const string & path, const string & title, int nb, double lo, double up) { return histogramFactory().createHistogram1D(path, title, nb, lo, up); } FactoryBase::tH1DPtr FactoryBase::createHistogram1D(const string & path, const string & title, const std::vector & edges) { return histogramFactory().createHistogram1D(path, title, edges); } FactoryBase::tH2DPtr FactoryBase::createHistogram2D(const string & path, int nbx, double xlo, double xup, int nby, double ylo, double yup) { return histogramFactory().createHistogram2D(path, nbx, xlo, xup, nby, ylo, yup); } FactoryBase::tH2DPtr FactoryBase::createHistogram2D(const string & path, const string & title, int nbx, double xlo, double xup, int nby, double ylo, double yup) { return histogramFactory().createHistogram2D(path, title, nbx, xlo, xup, nby, ylo, yup); } FactoryBase::tH2DPtr FactoryBase::createHistogram2D(const string & path, const string & title, const std::vector & xedges, const std::vector & yedges) { return histogramFactory().createHistogram2D(path, title, xedges, yedges); } FactoryBase::DataFiller FactoryBase::createDataSet(const string & path, const string & title, int dim) { return DataFiller(dataSetFactory().create(path, title, dim)); } void FactoryBase::registerClient(tIPtr client) { initrun(); clients.insert(client); } diff --git a/Analysis/LWH/DataPointSet.h b/Analysis/LWH/DataPointSet.h --- a/Analysis/LWH/DataPointSet.h +++ b/Analysis/LWH/DataPointSet.h @@ -1,372 +1,370 @@ // -*- C++ -*- // // DataPointSet.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef LWH_DataPointSet_H #define LWH_DataPointSet_H // // This is the declaration of the DataPointSet class representing // #include #include #include #include #include "AIDataPointSet.h" #include "ManagedObject.h" #include "DataPoint.h" namespace LWH { using namespace AIDA; /** * An DataPointSet represents a binned histogram axis. A 1D Histogram would have * one DataPointSet representing the X axis, while a 2D Histogram would have two * axes representing the X and Y DataPointSet. */ class DataPointSet: public IDataPointSet, public ManagedObject { public: /** * Standard constructor takes the dimension, \a D, of the data * points as argument. */ DataPointSet(int D): dim(D) {} /** * Destructor. */ virtual ~DataPointSet() {} /** * Not implemented in LWH. will throw an exception. */ IAnnotation & annotation() { throw std::runtime_error("LWH cannot handle annotations"); } /** * Not implemented in LWH. will throw an exception. */ const IAnnotation & annotation() const { throw std::runtime_error("LWH cannot handle annotations"); } /** * Get the data set's title. * @return The data set's title. */ std::string title() const { return theTitle; } /** * Get the data set's title. * @return The data set's title. */ std::string name() const { return theTitle; } /** * Set the data set's title. * @param title The title. * @return false If title cannot be changed. */ bool setTitle(const std::string & title) { theTitle = title; return true; } /** * Get the dimension of the IDataPoints that can be stored in the set. * @return The dimension of the IDataPoints storable in the set. * */ int dimension() const { return dim; } /** * Remove all the IDataPoints in the set. * After this the IDataPointSet is as just created. */ void clear() { dset.clear(); } /** * Get the current size of the IDataPointSet, i.e. the number * of IDataPoints contained in the set. * @return The size of the IDataPointSet. */ int size() const { return dset.size(); } /** * Get the IDataPoint at a give index in the set. * @param index The IDataPoint index. * @return The corresponding IDataPoint. */ IDataPoint * point(int index) { return &(dset[index]); } /** * Set the values and errors of a given coordinate all at once. If * this method is called on an empty IDataPointSet, a number of * points equal to the size of the arrays provided is created; if * the IDataPointSet is not empty the dimension of the array must * match with the size of the IDataPointSet. * @param coord The coordinate's index * @param val The array of the values for the given coordinate * @param err The array with the symmetric errors. * @return false if an illegal coordinate is provided or if there is * a mismatch between the size of the array and the size of the * IDataPointSet. */ bool setCoordinate(int coord, const std::vector & val, const std::vector & err) { return setCoordinate(coord, val, err, err); } /** * Set the values and errors of a given coordinate all at once. If * this method is called on an empty IDataPointSet, a number of * points equal to the size of the arrays provided is created; if * the IDataPointSet is not empty the dimension of the array must * match with the size of the IDataPointSet. * @param coord The coordinate's index * @param val The array of the values for the given coordinate * @param errp The array with the plus errors. * @param errm The array with the minus errors. * @return false if an illegal coordinate is provided or if there is * a mismatch between the size of the array and the size of the * IDataPointSet. * */ bool setCoordinate(int coord, const std::vector & val, const std::vector & errp, const std::vector & errm) { if ( coord < 0 || coord >= dimension() ) return false; if ( val.size() != dset.size() || errp.size() != dset.size() || errm.size() != dset.size() ) return false; for ( int i = 0, N = val.size(); i < N; ++i ) { dset[i].coordinate(coord)->setValue(val[i]); dset[i].coordinate(coord)->setErrorPlus(errp[i]); dset[i].coordinate(coord)->setErrorMinus(errm[i]); } return true; } /** * Return the data point at the given index. * @return 0 if index is out of range. */ const IDataPoint * point(int index) const { if ( index < 0 || unsigned(index) >= dset.size() ) return 0; return &(dset[index]); } /** * Add a new empty IDataPoint at the end of the set. * @return The newly added point. */ IDataPoint * addPoint() { dset.push_back(DataPoint(dimension())); return &(dset.back()); } /** * Add a copy of an IDataPoint at the end of the set. * @param point The IDataPoint to be added. * @return false If the point has the wrong dimension or * if the point cannot be added. */ bool addPoint(const IDataPoint & point) { if ( dimension() && dimension() != point.dimension() ) return false; dset.push_back(DataPoint(point)); return true; } /** * Remove the IDataPoint at a given index. * @param index The index of the IDataPoint to be removed. * @return false If the index is < 0 or >= size(). */ bool removePoint(int index) { if ( index < 0 || unsigned(index) >= dset.size() ) return false; dset.erase(dset.begin() + index); return true; } /** * Get the lower value for a give axis. * @param coord The coordinate of the axis. * @return The lower edge of the corresponding axis. * If coord < 0 or coord >= dimension(), or if the * set is empty NaN is returned. */ double lowerExtent(int coord) const { if ( dset.empty() ) return std::numeric_limits::quiet_NaN(); if ( coord < 0 || coord >= dimension() ) return std::numeric_limits::quiet_NaN(); double low = dset[0].coordinate(coord)->value(); for ( int i = 1, N = dset.size(); i < N; ++i ) low = std::min(low, dset[i].coordinate(coord)->value()); return low; } /** * Get the upper value for a give axis. * @param coord The coordinate of the axis. * @return The upper edge of the corresponding axis. * If coord < 0 or coord >= dimension(), or if the set * is empty NaN is returned. */ double upperExtent(int coord) const { if ( dset.empty() ) return std::numeric_limits::quiet_NaN(); if ( coord < 0 || coord >= dimension() ) return std::numeric_limits::quiet_NaN(); double upp = dset[0].coordinate(coord)->value(); for ( int i = 1, N = dset.size(); i < N; ++i ) upp = std::max(upp, dset[i].coordinate(coord)->value()); return upp; } /** * Scales the values and the errors of all the measurements * of each point by a given factor. * @param scale The scale factor. * @return false If an illegal scaleFactor is provided. */ bool scale(double scale) { for ( int i = 0, N = dset.size(); i < N; ++i ) for ( int j = 0, M = dset[i].dimension(); j < M; ++j ) { IMeasurement * m = dset[i].coordinate(j); m->setValue(m->value()*scale); m->setErrorPlus(m->errorPlus()*scale); m->setErrorMinus(m->errorPlus()*scale); } return true; } /** * Scales the values of all the measurements * of each point by a given factor. * @param scale The scale factor. * @return false If an illegal scaleFactor is provided. */ bool scaleValues(double scale) { for ( int i = 0, N = dset.size(); i < N; ++i ) for ( int j = 0, M = dset[i].dimension(); j < M; ++j ) { IMeasurement * m = dset[i].coordinate(j); m->setValue(m->value()*scale); } return true; } /** * Scales the errors of all the measurements * of each point by a given factor. * @param scale The scale factor. * @return false If an illegal scaleFactor is provided. */ bool scaleErrors(double scale) { for ( int i = 0, N = dset.size(); i < N; ++i ) for ( int j = 0, M = dset[i].dimension(); j < M; ++j ) { IMeasurement * m = dset[i].coordinate(j); m->setErrorPlus(m->errorPlus()*scale); m->setErrorMinus(m->errorPlus()*scale); } return true; } /** * Not implemented in LWH. * @return null pointer always. */ void * cast(const std::string &) const { return 0; } /** * Write out the data set in the AIDA xml format. */ bool writeXML(std::ostream & os, std::string path, std::string name) { os << " \n"; for ( int d = 0; d < dimension(); ++d ) os << " \n"; for ( int i = 0, N = size(); i < N; ++i ) { os << " \n"; for ( int j = 0, M = dimension(); j < M; ++j ) os << " coordinate(j)->value() << "\" errorPlus=\"" << point(i)->coordinate(j)->errorPlus() << "\" errorMinus=\"" << point(i)->coordinate(j)->errorMinus() << "\"/>\n"; os << " \n"; } os << " " << std::endl; return true; } /** * Write out the data set in a flat text file suitable for * eg. gnuplot to read. The coloums are layed out as 'x1 x2 ... xn * dx1+ dx2+ ... dxn+ dx1- dx2- ... dxn-'. */ bool writeFLAT(std::ostream & os, std::string path, std::string name) { os << "# " << path << "/" << name << " " << size() << " \"" << title() << " \" dimension " << dimension() << std::endl; for ( int i = 0, N = size(); i < N; ++i ) { for ( int j = 0, M = dimension(); j < M; ++j ) os << point(i)->coordinate(j)->value() << " "; for ( int j = 0, M = dimension(); j < M; ++j ) os << point(i)->coordinate(j)->errorPlus() << " "; for ( int j = 0, M = dimension(); j < M; ++j ) os << point(i)->coordinate(j)->errorMinus() << " "; os << std::endl; } os << std::endl; return true; } private: /** The title */ std::string theTitle; /** * The included data points. */ std::vector dset; /** * The dimension of the points in this set. */ unsigned int dim; - /** dummy pointer to non-existen annotation. */ - IAnnotation * anno; }; } #endif /* LWH_DataPointSet_H */ diff --git a/Analysis/LWH/Histogram1D.h b/Analysis/LWH/Histogram1D.h --- a/Analysis/LWH/Histogram1D.h +++ b/Analysis/LWH/Histogram1D.h @@ -1,523 +1,520 @@ // -*- C++ -*- // // Histogram1D.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef LWH_Histogram1D_H #define LWH_Histogram1D_H // // This is the declaration of the Histogram1D class. // #include "AIHistogram1D.h" #include "ManagedObject.h" #include "Axis.h" #include "VariAxis.h" #include #include namespace LWH { using namespace AIDA; /** * User level interface to 1D Histogram. */ class Histogram1D: public IHistogram1D, public ManagedObject { public: /** HistFactory is a friend. */ friend class HistogramFactory; public: /** * Standard constructor. */ Histogram1D(int n, double lo, double up) : fax(new Axis(n, lo, up)), vax(0), sum(n + 2), sumw(n + 2), sumw2(n + 2), sumxw(n + 2), sumx2w(n + 2) { ax = fax; } /** * Standard constructor for variable bin width. */ Histogram1D(const std::vector & edges) : fax(0), vax(new VariAxis(edges)), sum(edges.size() + 1), sumw(edges.size() + 1), sumw2(edges.size() + 1), sumxw(edges.size() + 1), sumx2w(edges.size() + 1) { ax = vax; } /** * Copy constructor. */ Histogram1D(const Histogram1D & h) : IBaseHistogram(h), IHistogram(h), IHistogram1D(h), ManagedObject(h), fax(0), vax(0), sum(h.sum), sumw(h.sumw), sumw2(h.sumw2), sumxw(h.sumxw), sumx2w(h.sumx2w) { const VariAxis * hvax = dynamic_cast(h.ax); if ( vax ) ax = vax = new VariAxis(*hvax); else ax = fax = new Axis(dynamic_cast(*h.ax)); } /// Destructor. virtual ~Histogram1D() { delete ax; } /** * Get the Histogram's title. * @return The Histogram's title. */ std::string title() const { return theTitle; } /** * Get the Histogram's title. * @return The Histogram's title. */ std::string name() const { return theTitle; } /** * Set the histogram title. * @param title The title. * @return false If title cannot be changed. */ bool setTitle(const std::string & title) { theTitle = title; return true; } /** * Not implemented in LWH. will throw an exception. */ IAnnotation & annotation() { throw std::runtime_error("LWH cannot handle annotations"); } /** * Not implemented in LWH. will throw an exception. */ const IAnnotation & annotation() const { throw std::runtime_error("LWH cannot handle annotations"); } /** * Get the Histogram's dimension. * @return The Histogram's dimension. */ int dimension() const { return 1; } /** * Reset the Histogram; as if just created. * @return false If something goes wrong. */ bool reset() { sum = std::vector(ax->bins() + 2); sumw = std::vector(ax->bins() + 2); sumxw = std::vector(ax->bins() + 2); sumx2w = std::vector(ax->bins() + 2); sumw2 = std::vector(ax->bins() + 2); return true; } /** * Get the number of in-range entries in the Histogram. * @return The number of in-range entries. * */ int entries() const { int si = 0; for ( int i = 2; i < ax->bins() + 2; ++i ) si += sum[i]; return si; } /** * Sum of the entries in all the IHistogram's bins, * i.e in-range bins, UNDERFLOW and OVERFLOW. * This is equivalent to the number of times the * method fill was invoked. * @return The sum of all the entries. */ int allEntries() const { return entries() + extraEntries(); } /** * Number of entries in the UNDERFLOW and OVERFLOW bins. * @return The number of entries outside the range of the IHistogram. */ int extraEntries() const { return sum[0] + sum[1]; } /** * Number of equivalent entries, * i.e. SUM[ weight ] ^ 2 / SUM[ weight^2 ] * @return The number of equivalent entries. */ double equivalentBinEntries() const { double sw = 0.0; double sw2 = 0.0; for ( int i = 2; i < ax->bins() + 2; ++i ) { sw += sumw[i]; sw2 += sumw2[i]; } return sw2/(sw*sw); } /** * Sum of in-range bin heights in the IHistogram, * UNDERFLOW and OVERFLOW bins are excluded. * @return The sum of the in-range bins heights. * */ double sumBinHeights() const { double sw = 0.0; for ( int i = 2; i < ax->bins() + 2; ++i ) sw += sumw[i]; return sw; } /** * Sum of the heights of all the IHistogram's bins, * i.e in-range bins, UNDERFLOW and OVERFLOW. * @return The sum of all the bins heights. */ double sumAllBinHeights() const { return sumBinHeights() + sumExtraBinHeights(); } /** * Sum of heights in the UNDERFLOW and OVERFLOW bins. * @return The sum of the heights of the out-of-range bins. */ double sumExtraBinHeights() const { return sumw[0] + sumw[1]; } /** * Minimum height of the in-range bins, * i.e. not considering the UNDERFLOW and OVERFLOW bins. * @return The minimum height among the in-range bins. */ double minBinHeight() const { double minw = sumw[2]; for ( int i = 3; i < ax->bins() + 2; ++i ) minw = std::min(minw, sumw[i]); return minw; } /** * Maximum height of the in-range bins, * i.e. not considering the UNDERFLOW and OVERFLOW bins. * @return The maximum height among the in-range bins. */ double maxBinHeight() const{ double maxw = sumw[2]; for ( int i = 3; i < ax->bins() + 2; ++i ) maxw = std::max(maxw, sumw[i]); return maxw; } /** * Fill the IHistogram1D with a value and the * corresponding weight. * @param x The value to be filled in. * @param weight The corresponding weight (by default 1). * @return false If the weight is <0 or >1 (?). */ bool fill(double x, double weight = 1.) { int i = ax->coordToIndex(x) + 2; ++sum[i]; sumw[i] += weight; sumxw[i] += x*weight; sumx2w[i] += x*x*weight; sumw2[i] += weight*weight; return weight >= 0 && weight <= 1; } /** * The weighted mean of a bin. * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW. * @return The mean of the corresponding bin. */ double binMean(int index) const { int i = index + 2; return sumw[i] != 0.0? sumxw[i]/sumw[i]: ( vax? vax->binMidPoint(index): fax->binMidPoint(index) ); }; /** * The weighted RMS of a bin. * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW. * @return The RMS of the corresponding bin. */ double binRms(int index) const { int i = index + 2; return sumw[i] == 0.0 || sum[i] < 2? ax->binWidth(index): std::sqrt(std::max(sumw[i]*sumx2w[i] - sumxw[i]*sumxw[i], 0.0))/sumw[i]; }; /** * Number of entries in the corresponding bin (ie the number of * times fill was called for this bin). * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW. * @return The number of entries in the corresponding bin. */ int binEntries(int index) const { return sum[index + 2]; } /** * Total height of the corresponding bin (ie the sum of the weights * in this bin). * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW. * @return The height of the corresponding bin. */ double binHeight(int index) const { return sumw[index + 2]; } /** * The error of a given bin. * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW. * @return The error on the corresponding bin. * */ double binError(int index) const { return std::sqrt(sumw2[index + 2]); } /** * The mean of the whole IHistogram1D. * @return The mean of the IHistogram1D. */ double mean() const { double s = 0.0; double sx = 0.0; for ( int i = 2; i < ax->bins() + 2; ++i ) { s += sumw[i]; sx += sumxw[i]; } return s != 0.0? sx/s: 0.0; } /** * The RMS of the whole IHistogram1D. * @return The RMS if the IHistogram1D. */ double rms() const { double s = 0.0; double sx = 0.0; double sx2 = 0.0; for ( int i = 2; i < ax->bins() + 2; ++i ) { s += sumw[i]; sx += sumxw[i]; sx2 += sumx2w[i]; } return s != 0.0? std::sqrt(std::max(s*sx2 - sx*sx, 0.0))/s: ax->upperEdge() - ax->lowerEdge(); } /** * Get the x axis of the IHistogram1D. * @return The x coordinate IAxis. */ const IAxis & axis() const { return *ax; } /** * Get the bin number corresponding to a given coordinate along the * x axis. This is a convenience method, equivalent to * axis().coordToIndex(coord). * @param coord The coordinalte along the x axis. * @return The corresponding bin number. */ int coordToIndex(double coord) const { return ax->coordToIndex(coord); } /** * Add to this Histogram1D the contents of another IHistogram1D. * @param h The Histogram1D to be added to this IHistogram1D. * @return false If the IHistogram1Ds binnings are incompatible. */ bool add(const Histogram1D & h) { if ( ax->upperEdge() != h.ax->upperEdge() || ax->lowerEdge() != h.ax->lowerEdge() || ax->bins() != h.ax->bins() ) return false; for ( int i = 0; i < ax->bins() + 2; ++i ) { sum[i] += h.sum[i]; sumw[i] += h.sumw[i]; sumxw[i] += h.sumxw[i]; sumx2w[i] += h.sumx2w[i]; sumw2[i] += h.sumw2[i]; } return true; } /** * Add to this IHistogram1D the contents of another IHistogram1D. * @param hist The IHistogram1D to be added to this IHistogram1D. * @return false If the IHistogram1Ds binnings are incompatible. */ bool add(const IHistogram1D & hist) { return add(dynamic_cast(hist)); } /** * Scale the contents of this histogram with the given factor. * @param s the scaling factor to use. */ bool scale(double s) { for ( int i = 0; i < ax->bins() + 2; ++i ) { sumw[i] *= s; sumxw[i] *= s; sumx2w[i] *= s; sumw2[i] *= s*s; } return true; } /** * Scale the given histogram so that the integral over all bins * (including overflow) gives \a intg. This function also corrects * for the bin-widths, which means that it should only be run once * for each histogram. Further rescaling must be done with the * scale(double) function. */ void normalize(double intg) { double oldintg = sumAllBinHeights(); if ( oldintg == 0.0 ) return; for ( int i = 0; i < ax->bins() + 2; ++i ) { double fac = intg/oldintg; if ( i >= 2 ) fac /= (ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2)); sumw[i] *= fac; sumxw[i] *= fac; sumx2w[i] *= fac; sumw2[i] *= fac*fac; } } /** * Return the integral over the histogram bins assuming it has been * normalize()d. */ double integral() const { double intg = sumw[0] + sumw[1]; for ( int i = 2; i < ax->bins() + 2; ++i ) intg += sumw[i]*(ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2)); return intg; } /** * Not implemented in LWH. * @return null pointer always. */ void * cast(const std::string &) const { return 0; } /** * Write out the histogram in the AIDA xml format. */ bool writeXML(std::ostream & os, std::string path, std::string name) { os << " \n upperEdge() << "\" numberOfBins=\"" << ax->bins() << "\" min=\"" << ax->lowerEdge() << "\" direction=\"x\""; if ( vax ) { os << ">\n"; for ( int i = 0, N = ax->bins() - 1; i < N; ++i ) os << " binUpperEdge(i) << "\"/>\n"; os << " \n"; } else { os << "/>\n"; } os << " \n \n \n \n"; for ( int i = 0; i < ax->bins() + 2; ++i ) if ( sum[i] ) { os << " \n"; } os << " \n " << std::endl; return true; } /** * Write out the histogram in a flat text file suitable for * eg. gnuplot to read. The coloums are layed out as 'x w w2 n'. */ bool writeFLAT(std::ostream & os, std::string path, std::string name) { os << "# " << path << "/" << name << " " << ax->lowerEdge() << " " << ax->bins() << " " << ax->upperEdge() << " \"" << title() << " \"" << std::endl; for ( int i = 2; i < ax->bins() + 2; ++i ) os << 0.5*(ax->binLowerEdge(i - 2) + ax->binUpperEdge(i - 2)) << " " << sumw[i] << " " << sqrt(sumw2[i]) << " " << sum[i] << std::endl; os << std::endl; return true; } private: /** The title */ std::string theTitle; /** The axis. */ IAxis * ax; /** Pointer (possibly null) to a axis with fixed bin width. */ Axis * fax; /** Pointer (possibly null) to a axis with fixed bin width. */ VariAxis * vax; /** The counts. */ std::vector sum; /** The weights. */ std::vector sumw; /** The squared weights. */ std::vector sumw2; /** The weighted x-values. */ std::vector sumxw; /** The weighted x-square-values. */ std::vector sumx2w; - /** dummy pointer to non-existen annotation. */ - IAnnotation * anno; - }; } #endif /* LWH_Histogram1D_H */ diff --git a/Config/std.h b/Config/std.h --- a/Config/std.h +++ b/Config/std.h @@ -1,187 +1,188 @@ // -*- C++ -*- // // std.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 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 std:: 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 std:: 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 * ThePEG_ALTERNATE_CONFIG. */ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include 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 : public binary_function { /** * 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::array; 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_fn; 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 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)); } /** factorial */ inline constexpr long double factorial(unsigned int n) { return (n < 2) ? 1.0 : n * factorial(n - 1); } /** Check if a given object is a part of a container. */ template 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 inline bool member(const vector & v, const Key & k) { for ( typename vector::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 inline std::insert_iterator inserter(Cont & c) { return std::insert_iterator(c, c.end()); } /** Return an insert iterator for a given vector. Overrides the * general version. */ template inline std::back_insert_iterator< vector > inserter(vector & v) { return back_inserter(v); } /** Return an insert iterator for a given vector. Overrides the * general version. */ template inline std::back_insert_iterator< deque > inserter(deque & 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; } } /** Macro for declaring a set. */ #define ThePEG_DECLARE_SET(VALTYPE,NAME) \ /** A set of VALTYPE. */ \ typedef set > NAME /** Macro for declaring a multiset. */ #define ThePEG_DECLARE_MULTISET(VALTYPE,NAME) \ /** A multiset of VALTYPE. */ \ typedef multiset > NAME /** Macro for declaring a map. */ #define ThePEG_DECLARE_MAP(KEYTYPE,VALTYPE,NAME) \ /** A map of VALTYPE indexed by KEYTYPE. */ \ typedef map > NAME #endif /* ThePEG_std_H */ diff --git a/Cuts/Cuts.cc b/Cuts/Cuts.cc --- a/Cuts/Cuts.cc +++ b/Cuts/Cuts.cc @@ -1,663 +1,663 @@ // -*- C++ -*- // // Cuts.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 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 Cuts class. // #include "Cuts.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/EventRecord/SubProcess.h" #include "ThePEG/EventRecord/Collision.h" #include "ThePEG/EventRecord/TmpTransform.h" #include "ThePEG/Utilities/UtilityBase.h" #include "ThePEG/Utilities/HoldFlag.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/Repository/CurrentGenerator.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Config/algorithm.h" using namespace ThePEG; Cuts::Cuts(Energy MhatMin) : theSMax(ZERO), theY(0), theCurrentSHat(-1.0*GeV2), theCurrentYHat(0), theMHatMin(MhatMin), theMHatMax(Constants::MaxEnergy), theYHatMin(-Constants::MaxRapidity), theYHatMax(Constants::MaxRapidity), theX1Min(0.0), theX1Max(1.0), theX2Min(0.0), theX2Max(1.0), theScaleMin(ZERO), theScaleMax(Constants::MaxEnergy2), theSubMirror(false), theCutWeight(1.0), theLastCutWeight(1.0) {} Cuts::~Cuts() {} IBPtr Cuts::clone() const { return new_ptr(*this); } IBPtr Cuts::fullclone() const { return new_ptr(*this); } void Cuts::doinitrun() { Interfaced::doinitrun(); if ( Debug::level ) { describe(); - for_each(theOneCuts, mem_fun(&OneCutBase::describe)); - for_each(theTwoCuts, mem_fun(&TwoCutBase::describe)); - for_each(theMultiCuts, mem_fun(&MultiCutBase::describe)); + for_each(theOneCuts, mem_fn(&OneCutBase::describe)); + for_each(theTwoCuts, mem_fn(&TwoCutBase::describe)); + for_each(theMultiCuts, mem_fn(&MultiCutBase::describe)); } } void Cuts::describe() const { CurrentGenerator::log() << fullName() << ":\n" << "MHat = " << theMHatMin/GeV << " .. " << theMHatMax/GeV << " GeV\n" << "Scale = " << theScaleMin/GeV2 << " .. " << theScaleMax/GeV2 << " GeV2\n" << "YHat = " << theYHatMin << " .. " << theYHatMax << '\n' << "X1 = " << theX1Min << " .. " << theX1Max << '\n' << "X2 = " << theX2Min << " .. " << theX2Max << "\n\n"; if ( theJetFinder ) theJetFinder->describe(); } void Cuts::initialize(Energy2 smax, double Y) { theSMax = smax; theMHatMax = min(theMHatMax, sqrt(smax)); theY = Y; theSubMirror = false; } void Cuts::initEvent() { theCurrentSHat = -1.0*GeV2; theCurrentYHat = 0.0; theSubMirror = false; } bool Cuts::initSubProcess(Energy2 shat, double yhat, bool mirror) const { theCutWeight = 1.0; theLastCutWeight = 1.0; theSubMirror = mirror; theCurrentSHat = shat; theCurrentYHat = yhat; if ( shat <= sHatMin() || shat > sHatMax()*(1.0 + 1000.0*Constants::epsilon) ) { theCutWeight = 0.0; return false; } if ( yhat <= yHatMin() || yhat >= yHatMax() ) { theCutWeight = 0.0; return false; } double x1 = min(1.0, sqrt(shat/SMax())*exp(yhat)); if ( x1 <= x1Min() || x1 > x1Max() ) { theCutWeight = 0.0; return false; } double x2 = min(1.0, sqrt(shat/SMax())*exp(-yhat)); if ( x2 <= x2Min() || x2 > x2Max() ) { theCutWeight = 0.0; return false; } return true; } bool Cuts::passCuts(const tcPDVector & ptype, const vector & p, tcPDPtr t1, tcPDPtr t2) const { if ( subMirror() ) { vector pmir = p; for ( int i = 0, N = pmir.size(); i < N; ++i ) pmir[i].setZ(-pmir[i].z()); swap(t1,t2); HoldFlag<> nomir(theSubMirror, false); return passCuts(ptype, pmir, t1, t2); } bool pass = true; theCutWeight = 1.0; theLastCutWeight = 1.0; if ( jetFinder() ) { if ( ptype.size() > jetFinder()->minOutgoing() ) { vector jets; tcPDVector jettype; if ( !jetFinder()->restrictConsitutents() ) { jets = p; jettype = ptype; } else { tcPDVector::const_iterator pd = ptype.begin(); vector::const_iterator pm = p.begin(); for ( ; pd != ptype.end(); ++pd, ++pm ) { if ( pm->rapidity() > jetFinder()->constituentRapidityRange().first && pm->rapidity() < jetFinder()->constituentRapidityRange().second ) { jets.push_back(*pm); jettype.push_back(*pd); } } } if ( jetFinder()->cluster(jettype,jets,this,t1,t2) ){ return passCuts(jettype,jets,t1,t2); } } } for ( int i = 0, N = p.size(); i < N; ++i ) for ( int j = 0, M = theOneCuts.size(); j < M; ++j ) { pass &= theOneCuts[j]->passCuts(this, ptype[i], p[i]); theCutWeight *= theLastCutWeight; theLastCutWeight = 1.0; if ( !pass ) { theCutWeight = 0.0; return false; } } for ( int i1 = 0, N1 = p.size() - 1; i1 < N1; ++i1 ) for ( int i2 = i1 + 1, N2 = p.size(); i2 < N2; ++i2 ) for ( int j = 0, M = theTwoCuts.size(); j < M; ++j ) { pass &= theTwoCuts[j]->passCuts(this, ptype[i1], ptype[i2], p[i1], p[i2]); theCutWeight *= theLastCutWeight; theLastCutWeight = 1.0; if ( !pass ) { theCutWeight = 0.0; return false; } } for ( int j = 0, M = theMultiCuts.size(); j < M; ++j ) { pass &= theMultiCuts[j]->passCuts(this, ptype, p); theCutWeight *= theLastCutWeight; theLastCutWeight = 1.0; if ( !pass ) { theCutWeight = 0.0; return false; } } if ( t1 ) { LorentzMomentum p1(ZERO, ZERO, 0.5*sqrt(currentSHat()), 0.5*sqrt(currentSHat())); for ( int i = 0, N = p.size(); i < N; ++i ) for ( int j = 0, M = theTwoCuts.size(); j < M; ++j ) { pass &= theTwoCuts[j]->passCuts(this, t1, ptype[i], p1, p[i], true, false); theCutWeight *= theLastCutWeight; theLastCutWeight = 1.0; if ( !pass ) { theCutWeight = 0.0; return false; } } } if ( t2 ) { LorentzMomentum p2(ZERO, ZERO, -0.5*sqrt(currentSHat()), 0.5*sqrt(currentSHat())); for ( int i = 0, N = p.size(); i < N; ++i ) for ( int j = 0, M = theTwoCuts.size(); j < M; ++j ) { pass &= theTwoCuts[j]->passCuts(this, ptype[i], t2, p[i], p2, false, true); theCutWeight *= theLastCutWeight; theLastCutWeight = 1.0; if ( !pass ) { theCutWeight = 0.0; return false; } } } return pass; } bool Cuts::passCuts(const tcPVector & p, tcPDPtr t1, tcPDPtr t2) const { tcPDVector ptype(p.size()); vector mom(p.size()); for ( int i = 0, N = p.size(); i < N; ++i ) { ptype[i] = p[i]->dataPtr(); mom[i] = p[i]->momentum(); } return passCuts(ptype, mom, t1, t2); } bool Cuts::passCuts(const SubProcess & sub) const { if ( !passCuts(tcPVector(sub.outgoing().begin(), sub.outgoing().end()), sub.incoming().first->dataPtr(), sub.incoming().second->dataPtr()) ) return false; return true; } bool Cuts::passCuts(const Collision & coll) const { tSubProPtr sub = coll.primarySubProcess(); LorentzMomentum phat = sub->incoming().first->momentum() + sub->incoming().second->momentum(); if ( !initSubProcess(phat.m2(), phat.rapidity()) ) return false; TmpTransform tmp(sub, Utilities::getBoostToCM(sub->incoming())); if ( !passCuts(*sub) ) return false; return true; } Energy2 Cuts::minS(const tcPDVector & pv) const { Energy2 mins = ZERO; for ( int i = 0, N = theMultiCuts.size(); i < N; ++i ) mins = max(mins, theMultiCuts[i]->minS(pv)); return mins; } Energy2 Cuts::maxS(const tcPDVector & pv) const { Energy2 maxs = SMax(); for ( int i = 0, N = theMultiCuts.size(); i < N; ++i ) maxs = min(maxs, theMultiCuts[i]->maxS(pv)); return maxs; } Energy2 Cuts::minSij(tcPDPtr pi, tcPDPtr pj) const { Energy2 mins = ZERO; for ( int i = 0, N = theTwoCuts.size(); i < N; ++i ) mins = max(mins, theTwoCuts[i]->minSij(pi, pj)); if ( mins > ZERO ) return mins; mins = sqr(pi->massMin() + pj->massMin()); mins = max(mins, sqr(minKTClus(pi, pj))/4.0); mins = max(mins, minDurham(pi, pj)*currentSHat()/2.0); mins = max(mins, minKT(pi)*minKT(pj)*minDeltaR(pi, pj)/4.0); return mins; } Energy2 Cuts::minTij(tcPDPtr pi, tcPDPtr po) const { Energy2 mint = ZERO; for ( int i = 0, N = theTwoCuts.size(); i < N; ++i ) mint = max(mint, theTwoCuts[i]->minTij(pi, po)); if ( mint > ZERO ) return mint; mint = max(mint, sqr(minKT(po))); return mint; } double Cuts::minDeltaR(tcPDPtr pi, tcPDPtr pj) const { double mindr = 0.0; for ( int i = 0, N = theTwoCuts.size(); i < N; ++i ) mindr = max(mindr, theTwoCuts[i]->minDeltaR(pi, pj)); return mindr; } Energy Cuts::minKTClus(tcPDPtr pi, tcPDPtr pj) const { Energy minkt = ZERO; for ( int i = 0, N = theTwoCuts.size(); i < N; ++i ) minkt = max(minkt, theTwoCuts[i]->minKTClus(pi, pj)); return minkt; } double Cuts::minDurham(tcPDPtr pi, tcPDPtr pj) const { double y = 0.0; for ( int i = 0, N = theTwoCuts.size(); i < N; ++i ) y = max(y, theTwoCuts[i]->minDurham(pi, pj)); return y; } Energy Cuts::minKT(tcPDPtr p) const { Energy minkt = ZERO; for ( int i = 0, N = theOneCuts.size(); i < N; ++i ) minkt = max(minkt, theOneCuts[i]->minKT(p)); if ( minkt > ZERO ) return minkt; minkt = minKTClus(p, tcPDPtr()); return minkt; } double Cuts::minEta(tcPDPtr p) const { double mineta = -Constants::MaxRapidity; for ( int i = 0, N = theOneCuts.size(); i < N; ++i ) mineta = max(mineta, theOneCuts[i]->minEta(p)); return mineta; } double Cuts::maxEta(tcPDPtr p) const { double maxeta = Constants::MaxRapidity; for ( int i = 0, N = theOneCuts.size(); i < N; ++i ) maxeta = min(maxeta, theOneCuts[i]->maxEta(p)); return maxeta; } double Cuts::minYStar(tcPDPtr p) const { if ( currentSHat() < ZERO ) return -Constants::MaxRapidity; if ( subMirror() ) { HoldFlag<> nomir(theSubMirror, false); return -maxYStar(p); } double etamin = minEta(p); double ytot = Y() + currentYHat(); if ( etamin > 0.0 ) { Energy minkt = minKT(p); Energy maxm = p->massMax(); return asinh(minkt*sinh(etamin)/sqrt(sqr(minkt) + sqr(maxm))) - ytot; } else { return etamin - ytot; } } double Cuts::maxYStar(tcPDPtr p) const { if ( currentSHat() < ZERO ) return Constants::MaxRapidity; if ( subMirror() ) { HoldFlag<> nomir(theSubMirror, false); return -minYStar(p); } double etamax = maxEta(p); double ytot = Y() + currentYHat(); if ( etamax > 0.0 ) { return etamax - ytot; } else { Energy minkt = minKT(p); Energy maxm = p->massMax(); return asinh(minkt*sinh(etamax)/sqrt(sqr(minkt) + sqr(maxm))) - ytot; } } double Cuts::minRapidityMax(tcPDPtr p) const { double minRapidityMax = -Constants::MaxRapidity; for ( int i = 0, N = theOneCuts.size(); i < N; ++i ) minRapidityMax = max(minRapidityMax, theOneCuts[i]->minRapidityMax(p)); return minRapidityMax; } double Cuts::maxRapidityMin(tcPDPtr p) const { double maxRapidityMin = Constants::MaxRapidity; for ( int i = 0, N = theOneCuts.size(); i < N; ++i ) maxRapidityMin = min(maxRapidityMin, theOneCuts[i]->maxRapidityMin(p)); return maxRapidityMin; } void Cuts::persistentOutput(PersistentOStream & os) const { os << ounit(theSMax, GeV2) << theY << ounit(theCurrentSHat, GeV2) << theCurrentYHat << ounit(theMHatMin, GeV) << ounit(theMHatMax, GeV) << theYHatMin << theYHatMax << theX1Min << theX1Max << theX2Min << theX2Max << ounit(theScaleMin, GeV2) << ounit(theScaleMax, GeV2) << theOneCuts << theTwoCuts << theMultiCuts << theJetFinder << theSubMirror << theCutWeight << theLastCutWeight << theFuzzyTheta; } void Cuts::persistentInput(PersistentIStream & is, int) { is >> iunit(theSMax, GeV2) >> theY >> iunit(theCurrentSHat, GeV2) >> theCurrentYHat >> iunit(theMHatMin, GeV) >> iunit(theMHatMax, GeV) >> theYHatMin >> theYHatMax >> theX1Min >> theX1Max >> theX2Min >> theX2Max >> iunit(theScaleMin, GeV2) >> iunit(theScaleMax, GeV2) >> theOneCuts >> theTwoCuts >> theMultiCuts >> theJetFinder >> theSubMirror >> theCutWeight >> theLastCutWeight >> theFuzzyTheta; } ClassDescription Cuts::initCuts; // Definition of the static class description member. Energy Cuts::maxMHatMin() const { return theMHatMax; } Energy Cuts::minMHatMax() const { return theMHatMin; } Energy2 Cuts::maxScaleMin() const { return theScaleMax; } Energy2 Cuts::minScaleMax() const { return theScaleMin; } double Cuts::maxYHatMin() const { return theYHatMax; } double Cuts::minYHatMax() const { return theYHatMin; } double Cuts::maxX1Min() const { return theX1Max; } double Cuts::minX1Max() const { return theX1Min; } double Cuts::maxX2Min() const { return theX2Max; } double Cuts::minX2Max() const { return theX2Min; } void Cuts::Init() { typedef double (ThePEG::Cuts::*IGFN)() const; typedef void (ThePEG::Cuts::*ISFN)(double); static ClassDocumentation documentation ("Cuts is a class for implementing kinematical cuts in ThePEG. The " "class itself only implements cuts on the total momentum of the hard " "sub-process, implemented as minimum and maximum values of \\f$x_1\\f$ " "and \\f$x_2\\f$ (or \\f$\\hat{s}\\f$ and \\f$\\hat{y}\\f$. Further cuts " "can be implemented either by inheriting from this base class, in which " "the virtual cut() function should be overridden, or by assigning " "objects of class OneCutBase, TwoCutBase and MultiCutBase defining " "cuts on single particles, pairs of particles and groups of " "particles respectively."); static Parameter interfaceMHatMin ("MHatMin", "The minimum allowed value of \\f$\\sqrt{\\hat{s}}\\f$.", &Cuts::theMHatMin, GeV, 2.0*GeV, ZERO, Constants::MaxEnergy, true, false, Interface::limited, 0, 0, 0, &Cuts::maxMHatMin, 0); interfaceMHatMin.setHasDefault(false); static Parameter interfaceMHatMax ("MHatMax", "The maximum allowed value of \\f$\\sqrt{\\hat{s}}\\f$.", &Cuts::theMHatMax, GeV, 100.0*GeV, ZERO, ZERO, true, false, Interface::lowerlim, 0, 0, &Cuts::minMHatMax, 0, 0); interfaceMHatMax.setHasDefault(false); static Parameter interfaceScaleMin ("ScaleMin", "The minimum allowed value of the scale to be used in PDFs and " "coupling constants.", &Cuts::theScaleMin, GeV2, ZERO, ZERO, Constants::MaxEnergy2, true, false, Interface::limited, 0, 0, 0, &Cuts::maxScaleMin, 0); interfaceScaleMin.setHasDefault(false); static Parameter interfaceScaleMax ("ScaleMax", "The maximum allowed value of the scale to be used in PDFs and " "coupling constants.", &Cuts::theScaleMax, GeV2, 10000.0*GeV2, ZERO, ZERO, true, false, Interface::lowerlim, 0, 0, &Cuts::minScaleMax, 0, 0); interfaceScaleMax.setHasDefault(false); static Parameter interfaceYHatMin ("YHatMin", "The minimum value of the rapidity of the hard sub-process " "(wrt. the rest system of the colliding particles).", &Cuts::theYHatMin, -10.0, 0.0, Constants::MaxRapidity, true, false, Interface::upperlim, (ISFN)0, (IGFN)0, (IGFN)0, &Cuts::maxYHatMin, (IGFN)0); interfaceYHatMin.setHasDefault(false); static Parameter interfaceYHatMax ("YHatMax", "The maximum value of the rapidity of the hard sub-process " "(wrt. the rest system of the colliding particles).", &Cuts::theYHatMax, 10.0, -Constants::MaxRapidity, 0.0, true, false, Interface::lowerlim, (ISFN)0, (IGFN)0, &Cuts::minYHatMax, (IGFN)0, (IGFN)0); interfaceYHatMax.setHasDefault(false); static Parameter interfaceX1Min ("X1Min", "The minimum value of the positive light-cone fraction of the hard " "sub-process.", &Cuts::theX1Min, 0.0, 0.0, 1.0, true, false, Interface::limited, (ISFN)0, (IGFN)0, (IGFN)0, &Cuts::maxX1Min, (IGFN)0); interfaceX1Min.setHasDefault(false); static Parameter interfaceX1Max ("X1Max", "The maximum value of the positive light-cone fraction of the hard " "sub-process.", &Cuts::theX1Max, 0.0, 0.0, 1.0, true, false, Interface::limited, (ISFN)0, (IGFN)0, &Cuts::minX1Max, (IGFN)0, (IGFN)0); interfaceX1Max.setHasDefault(false); static Parameter interfaceX2Min ("X2Min", "The minimum value of the negative light-cone fraction of the hard " "sub-process.", &Cuts::theX2Min, 0.0, 0.0, 1.0, true, false, Interface::limited, (ISFN)0, (IGFN)0, (IGFN)0, &Cuts::maxX2Min, (IGFN)0); interfaceX2Min.setHasDefault(false); static Parameter interfaceX2Max ("X2Max", "The maximum value of the negative light-cone fraction of the hard " "sub-process.", &Cuts::theX2Max, 0.0, 0.0, 1.0, true, false, Interface::limited, (ISFN)0, (IGFN)0, &Cuts::minX2Max, (IGFN)0, (IGFN)0); interfaceX2Max.setHasDefault(false); static RefVector interfaceOneCuts ("OneCuts", "The objects defining cuts on single outgoing partons from the " "hard sub-process.", &Cuts::theOneCuts, -1, true, false, true, false, false); static RefVector interfaceTwoCuts ("TwoCuts", "The objects defining cuts on pairs of particles in the " "hard sub-process.", &Cuts::theTwoCuts, -1, true, false, true, false, false); static RefVector interfaceMultiCuts ("MultiCuts", "The objects defining cuts on sets of outgoing particles from the " "hard sub-process.", &Cuts::theMultiCuts, -1, true, false, true, false, false); static Reference interfaceJetFinder ("JetFinder", "Set a JetFinder object used to define cuts on the" "level of reconstructed jets as needed for higher order corrections.", &Cuts::theJetFinder, false, false, true, true, false); static Reference interfaceFuzzy ("Fuzzy", "The fuzziness to be applied to cuts (may not be supported by all cut objects).", &Cuts::theFuzzyTheta, false, false, true, true, false); interfaceX1Min.rank(10); interfaceX1Max.rank(9); interfaceX2Min.rank(8); interfaceX2Max.rank(7); interfaceMHatMin.rank(6); interfaceMHatMax.rank(5); interfaceYHatMin.rank(4); interfaceYHatMax.rank(3); interfaceOneCuts.rank(2); interfaceTwoCuts.rank(1); } double Cuts::yHatMin() const { return theX1Min > 0.0 && theX2Max > 0.0? max(theYHatMin, 0.5*log(theX1Min/theX2Max)): theYHatMin; } double Cuts::yHatMax() const { return theX1Max > 0.0 && theX2Min > 0.0? min(theYHatMax, 0.5*log(theX1Max/theX2Min)): theYHatMax; } bool Cuts::yHat(double y) const { return y > yHatMin() && y < yHatMax(); } double Cuts::x1Min() const { return max(theX1Min, (theMHatMin/sqrt(SMax()))*exp(theYHatMin)); } double Cuts::x1Max() const { return min(theX1Max, (theMHatMax/sqrt(SMax()))*exp(theYHatMax)); } bool Cuts::x1(double x) const { return x > x1Min() && x <= x1Max(); } double Cuts::x2Min() const { return max(theX2Min, (theMHatMin/sqrt(SMax()))/exp(theYHatMax)); } double Cuts::x2Max() const { return min(theX2Max, (theMHatMax/sqrt(SMax()))/exp(theYHatMin)); } bool Cuts::x2(double x) const { return x > x2Min() && x <= x2Max(); } template vector::transient_const_pointer> Cuts::oneCutObjects() const { typedef typename Ptr::transient_const_pointer tcPtr; vector ret; for ( int i = 0, N = theOneCuts.size(); i < N; ++i ) if ( dynamic_ptr_cast(theOneCuts[i]) ) ret.push_back(dynamic_ptr_cast(theOneCuts[i])); return ret; } template vector::transient_const_pointer> Cuts::twoCutObjects() const { typedef typename Ptr::transient_const_pointer tcPtr; vector ret; for ( int i = 0, N = theTwoCuts.size(); i < N; ++i ) if ( dynamic_ptr_cast(theTwoCuts[i]) ) ret.push_back(dynamic_ptr_cast(theTwoCuts[i])); return ret; } template vector::transient_const_pointer> Cuts::multiCutObjects() const { typedef typename Ptr::transient_const_pointer tcPtr; vector ret; for ( int i = 0, N = theMultiCuts.size(); i < N; ++i ) if ( dynamic_ptr_cast(theMultiCuts[i]) ) ret.push_back(dynamic_ptr_cast(theMultiCuts[i])); return ret; } diff --git a/Helicity/Vertex/AbstractSSSSVertex.cc b/Helicity/Vertex/AbstractSSSSVertex.cc --- a/Helicity/Vertex/AbstractSSSSVertex.cc +++ b/Helicity/Vertex/AbstractSSSSVertex.cc @@ -1,32 +1,26 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the AbstractSSSSVertex class. // #include "AbstractSSSSVertex.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" using namespace ThePEG; using namespace Helicity; // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeThePEGAbstractSSSSVertex("ThePEG::AbstractSSSSVertex", "libThePEG.so"); void AbstractSSSSVertex::Init() { static ClassDocumentation documentation ("The AbstractSSSSVertex class is the base class for all " "scalar-scalar-scalar-scalar interactions"); } -ScalarWaveFunction AbstractSSSSVertex::evaluate(Energy2,int, tcPDPtr, - const ScalarWaveFunction & , - const ScalarWaveFunction & , - const ScalarWaveFunction & ) { - assert(false); -} diff --git a/Helicity/Vertex/AbstractSSSSVertex.h b/Helicity/Vertex/AbstractSSSSVertex.h --- a/Helicity/Vertex/AbstractSSSSVertex.h +++ b/Helicity/Vertex/AbstractSSSSVertex.h @@ -1,88 +1,88 @@ // -*- C++ -*- #ifndef HELICITY_AbstractSSSSVertex_H #define HELICITY_AbstractSSSSVertex_H // // This is the declaration of the AbstractSSSSVertex class. // #include "VertexBase.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "AbstractSSSSVertex.fh" namespace ThePEG { namespace Helicity { /** * The AbstractSSSSVertex class is the base class for all scalar-scalar-scalar * interactions in ThePEG */ class AbstractSSSSVertex: public VertexBase { public: /** * Default constructor */ AbstractSSSSVertex() : VertexBase(VertexType::SSSS) {} /** * Members to calculate the helicity amplitude expressions for vertices * and off-shell particles. */ //@{ /** * Evaluate the vertex. * @param q2 The scale \f$q^2\f$ for the coupling at the vertex. * @param sca1 The wavefunction for the first scalar. * @param sca2 The wavefunction for the second scalar. * @param sca3 The wavefunction for the third scalar. * @param sca4 The wavefunction for the fourth scalar. */ virtual Complex evaluate(Energy2 q2,const ScalarWaveFunction & sca1, const ScalarWaveFunction & sca2, const ScalarWaveFunction & sca3, const ScalarWaveFunction & sca4) = 0; /** * Evaluate the off-shell scalar coming from the vertex. * @param q2 The scale \f$q^2\f$ for the coupling at the vertex. * @param iopt Option of the shape of the Breit-Wigner for the off-shell scalar. * @param out The ParticleData pointer for the off-shell scalar. * @param sca1 The wavefunction for the first scalar. * @param sca2 The wavefunction for the second scalar. * @param sca3 The wavefunction for the third scalar. */ virtual ScalarWaveFunction evaluate(Energy2 q2,int iopt, tcPDPtr out, const ScalarWaveFunction & sca1, const ScalarWaveFunction & sca2, - const ScalarWaveFunction & sca3); + const ScalarWaveFunction & sca3)=0; //@} public: /** * 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. */ AbstractSSSSVertex & operator=(const AbstractSSSSVertex &) = delete; }; } } namespace ThePEG { } #endif /* HELICITY_AbstractSSSSVertex_H */ diff --git a/PDT/QuarksToHadronsDecayer.cc b/PDT/QuarksToHadronsDecayer.cc --- a/PDT/QuarksToHadronsDecayer.cc +++ b/PDT/QuarksToHadronsDecayer.cc @@ -1,235 +1,236 @@ // -*- C++ -*- // // QuarksToHadronsDecayer.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 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 QuarksToHadronsDecayer class. // #include "QuarksToHadronsDecayer.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/PDT/StandardMatchers.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/SimplePhaseSpace.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace ThePEG; QuarksToHadronsDecayer::~QuarksToHadronsDecayer() {} IBPtr QuarksToHadronsDecayer::clone() const { return new_ptr(*this); } IBPtr QuarksToHadronsDecayer::fullclone() const { return new_ptr(*this); } bool QuarksToHadronsDecayer::accept(const DecayMode & dm) const { int col = 0; int acol = 0; if ( !dm.productMatchers().empty() ) { for ( MatcherMSet::const_iterator it = dm.productMatchers().begin(); it != dm.productMatchers().end(); ++it ) { - if ( typeid(**it) == typeid(MatchLightQuark) ) ++col; - else if ( typeid(**it) == typeid(MatchLightAntiQuark) ) ++acol; + const auto & tmp=**it; + if ( typeid(tmp) == typeid(MatchLightQuark) ) ++col; + else if ( typeid(tmp) == typeid(MatchLightAntiQuark) ) ++acol; else return false; } if ( col != 1 || col != acol ) return false; } if ( dm.orderedProducts().size() + col + acol < 2 || !dm.cascadeProducts().empty() || dm.wildProductMatcher() ) return false; for ( int i = 0, N = dm.orderedProducts().size(); i < N; ++i ) { if ( DiquarkMatcher::Check(*dm.orderedProducts()[i]) ) { if ( i + 1 != N ) return false; if ( dm.orderedProducts()[i]->id() < 0 ) ++col; else ++acol; } if ( QuarkMatcher::Check(*dm.orderedProducts()[i]) ) { if ( dm.orderedProducts()[i]->id() > 0 ) ++col; else ++acol; } } if ( acol != col || col < 1 || col > 2 ) return false; return true; } PVector QuarksToHadronsDecayer::decay(const DecayMode & dm, const Particle & parent) const { PVector children; tcPDVector quarks; if ( !dm.productMatchers().empty() ) { tcPDPtr pd = getParticleData(flavourGenerator()->selectQuark()); quarks.push_back(pd); quarks.push_back(pd->CC()); } Energy summq = ZERO; Energy summp = ZERO; tPDVector prods = dm.orderedProducts(); for ( int i = 0, N = prods.size(); i < N; ++i ) if ( QuarkMatcher::Check(*prods[i]) || DiquarkMatcher::Check(*prods[i])) { quarks.push_back(prods[i]); summq += quarks.back()->mass(); } else { children.push_back(prods[i]->produceParticle()); summp += children.back()->mass(); } Energy summh = ZERO; PVector hadrons; if ( !quarks.empty() ) do { hadrons = getHadrons(getN(parent.mass(), summq, quarks.size()), quarks); summh = ZERO; for ( int i = 0, N = hadrons.size(); i < N; ++i ) summh += hadrons[i]->mass(); } while ( hadrons.empty() || summp + summh >= parent.mass() ); children.insert(children.end(), hadrons.begin(), hadrons.end()); distribute(parent, children); finalBoost(parent, children); setScales(parent, children); return children; } int QuarksToHadronsDecayer::getN(Energy m0, Energy summq, int Nq) const { int Nh = fixedN(); if ( Nh >= 2 ) return Nh; double c = c1()*log((m0 - summq)/c2()) + c3(); if ( c < 0.0 ) return minN(); while ( true ) { using namespace Constants; Nh = int(0.5 + double(Nq)/4.0 + c + sqrt(-2.0*c*log(max(1.0e-10, rnd())))*sin(2.0*pi*rnd())); if ( Nh >= minN() ) return Nh; } } PVector QuarksToHadronsDecayer:: getHadrons(int Nh, tcPDVector quarks) const { PVector hadrons; Nh -= quarks.size()/2; while ( Nh-- > 0 ) { int i = irnd(quarks.size() - 1); tcPDPair hq = flavourGenerator()->alwaysGenerateHadron(quarks[i]); hadrons.push_back(hq.first->produceParticle()); quarks[i] = hq.second; } if ( DiquarkMatcher::Check(*quarks[0]) && DiquarkMatcher::Check(*quarks[1]) ) return PVector(); tcPDPtr h = flavourGenerator()->alwaysGetHadron(quarks[0], quarks[1]); hadrons.push_back(h->produceParticle()); if ( quarks.size() <= 2 ) return hadrons; if ( DiquarkMatcher::Check(*quarks[2]) && DiquarkMatcher::Check(*quarks[3]) ) return PVector(); h = flavourGenerator()->alwaysGetHadron(quarks[2], quarks[3]); hadrons.push_back(h->produceParticle()); return hadrons; } void QuarksToHadronsDecayer:: distribute(const Particle & parent, PVector & children) const { do { try { SimplePhaseSpace::CMSn(children, parent.mass()); } catch ( ImpossibleKinematics & e) { children.clear(); return; } } while ( reweight(parent, children) < rnd() ); } double QuarksToHadronsDecayer:: reweight(const Particle &, const PVector &) const { return 1.0; } void QuarksToHadronsDecayer::persistentOutput(PersistentOStream & os) const { os << theFixedN << theMinN << theC1 << ounit(theC2,GeV) << theC3 << theFlavourGenerator; } void QuarksToHadronsDecayer::persistentInput(PersistentIStream & is, int) { is >> theFixedN >> theMinN >> theC1 >> iunit(theC2,GeV) >> theC3 >> theFlavourGenerator; } ClassDescription QuarksToHadronsDecayer::initQuarksToHadronsDecayer; // Definition of the static class description member. void QuarksToHadronsDecayer::Init() { static ClassDocumentation documentation ("This class decays particles to nq (2 or 4) quarks which then are " "decayes to hadrons according to phase space. The number of final " "hadrons can either be given by a fixed number or as a Gaussian " "multiplicity distribution centered around c+nq/4+c3 and a width " "sqrt(c), where c = c1 log((m - summ)/c2), m is the mass of the " "decaying particle, summ the sum of the quark masses and ci real " "parameters."); static Parameter interfaceFixedN ("FixedN", "The fixed number of hadrons to be produced. If less than 2, the " "number is instead given by a gaussian multiplicity distribution.", &QuarksToHadronsDecayer::theFixedN, 0, 0, 10, true, false, true); static Parameter interfaceMinN ("MinN", "The minimum hadrons to be produced.", &QuarksToHadronsDecayer::theMinN, 2, 2, 10, true, false, true); static Parameter interfaceC1 ("C1", "The c1 parameter of the gaussian multiplicity distribution centered " "around c1 log((m - summ)/c2) +c3.", &QuarksToHadronsDecayer::theC1, 4.5, 0.0, 10.0, true, false, true); static Parameter interfaceC2 ("C2", "The c2 parameter of the gaussian multiplicity distribution centered " "around c1 log((m - summ)/c2) +c3.", &QuarksToHadronsDecayer::theC2, GeV, 0.7*GeV, ZERO, 10.0*GeV, true, false, true); static Parameter interfaceC3 ("C3", "The c3 parameter of the gaussian multiplicity distribution centered " "around c1 log((m - summ)/c2) +c3.", &QuarksToHadronsDecayer::theC3, 0.0, 0.0, 10.0, true, false, true); static Reference interfaceFlavourGenerator ("FlavourGenerator", "The object in charge of generating hadrons spieces from given quark " "flavours.", &QuarksToHadronsDecayer::theFlavourGenerator, true, false, true, false, true); interfaceFixedN.rank(10); interfaceMinN.rank(9); interfaceFlavourGenerator.rank(8); interfaceMinN.setHasDefault(false);; } diff --git a/Persistency/PersistentOStream.cc b/Persistency/PersistentOStream.cc --- a/Persistency/PersistentOStream.cc +++ b/Persistency/PersistentOStream.cc @@ -1,169 +1,170 @@ // -*- C++ -*- // // PersistentOStream.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 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 PersistentOStream class. // #include "PersistentOStream.h" #include "ThePEG/Utilities/DynamicLoader.h" #include namespace ThePEG { PersistentOStream::PersistentOStream(ostream & os, const vector & libs) : theOStream(&os), badState(false), allocStream(false) { init(libs); } PersistentOStream::PersistentOStream(string file, const vector & libs) : badState(false), allocStream(true) { // if ( file[0] == '|' ) // theOStream = new opfstream(file.substr(1).c_str()); // else if ( file.substr(file.length()-3, file.length()) == ".gz" ) // theOStream = new opfstream(string("gzip > " + file).c_str()); // else theOStream = new ofstream(file.c_str()); if ( theOStream ) init(libs); else setBadState(); } void PersistentOStream::init(const vector & libs) { operator<<(string("ThePEG version 1 Database")); operator<<(version); operator<<(subVersion); *this << DynamicLoader::appendedPaths(); *this << DynamicLoader::prependedPaths(); vector libraries; for ( int i = 0, N = libs.size(); i < N; ++i ) libraries.push_back(DynamicLoader::dlnameversion(libs[i])); *this << libraries; } PersistentOStream::~PersistentOStream() { if ( allocStream ) delete theOStream; } void PersistentOStream:: putObjectPart(tcBPtr obj, const ClassDescriptionBase * db) { ClassDescriptionBase::DescriptionVector::const_iterator bit = db->descriptions().begin(); while ( bit != db->descriptions().end() ) { putObjectPart(obj, *bit++); endBase(); } db->output(obj, *this); } PersistentOStream & PersistentOStream::outputPointer(tcBPtr obj) { if ( !good() ) return *this; if ( !obj ) return operator<<(0); // It it's the null pointer, just print a zero. int oid = 0; const ClassDescriptionBase * desc = 0; try { // Check if the object has been written before. In that case just write // out it's number ObjectMap::const_iterator oit = writtenObjects.find(obj); if ( oit != writtenObjects.end() ) { *this << oit->second; return *this; } // This object hasn't been written before so we write it out, beginning // with a number, then the class information, and finally let it write // itself on the stream. beginObject(); oid = writtenObjects.size()+1; writtenObjects[obj] = oid; *this << oid; desc = writeClassId(obj); *this << obj->uniqueId; putObjectPart(obj, desc); endObject(); } catch (Exception & e) { e.handle(); string classname = ""; if ( desc ) classname = desc->name(); throw WriteError() << "While writing object number " << oid << " of class " << classname << ":\n" << e.message() << Exception::runerror; setBadState(); } catch (...) { setBadState(); } checkState(); return *this; } const ClassDescriptionBase * PersistentOStream::writeClassId(tcBPtr obj) { - const ClassDescriptionBase * db = DescriptionList::find(typeid(*obj)); + const auto & tmp=*obj; + const ClassDescriptionBase * db = DescriptionList::find(typeid(tmp)); if ( !db ) { throw MissingClass() << "PersistentOStream could not find the ClassDescription object " - << "corresponding to the class " << typeid(*obj).name() + << "corresponding to the class " << typeid(tmp).name() << ". Please check that the class has a properly instantiated " << "ClassDescription object." << Exception::runerror; } writeClassDescription(db); return db; } void PersistentOStream:: writeClassDescription(const ClassDescriptionBase * db) { // If objects of this class has been written out before, just write // the corresponding number ClassMap::iterator cit = writtenClasses.find(db); if ( cit != writtenClasses.end() ) { operator<<(cit->second); return; } // This class hasn't been written before, so append it to the list of // written classes and assign a number to it, before writing the string // containing the information int cid = writtenClasses.size(); writtenClasses[db] = cid; operator<<(cid); operator<<(db->name()); operator<<(db->version()); operator<<(DynamicLoader::dlnameversion(db->library())); // Now write its base classes or a zero if the base class is PersistentBase. operator<<(db->descriptions().size()); DescriptionVector::const_iterator bit = db->descriptions().begin(); while ( bit != db->descriptions().end() ) writeClassDescription(*bit++); } PersistentOStream & PersistentOStream::flush() { #ifdef _LIBCPP_VERSION typedef ObjectMap::const_iterator Iterator; #else typedef ObjectMap::iterator Iterator; #endif Iterator it = writtenObjects.begin(); while ( it != writtenObjects.end() ) { Iterator it2 = it++; if ( (*it2).second > lastSavedObject.top() ) writtenObjects.erase(it2); } os().flush(); return *this; } } diff --git a/Repository/BaseRepository.cc b/Repository/BaseRepository.cc --- a/Repository/BaseRepository.cc +++ b/Repository/BaseRepository.cc @@ -1,977 +1,986 @@ // -*- C++ -*- // // BaseRepository.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 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 BaseRepository class. // // macro is passed in from -D compile flag #ifndef THEPEG_PKGDATADIR #error Makefile.am needs to define THEPEG_PKGDATADIR #endif #include "BaseRepository.h" #include "ThePEG/Config/algorithm.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/InterfaceBase.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/Utilities/ClassDescription.h" #include "ThePEG/Utilities/DescriptionList.h" #include "ThePEG/Utilities/HoldFlag.h" #include "ThePEG/Utilities/TypeInfo.h" #include "ThePEG/Utilities/DynamicLoader.h" #include "ThePEG/Utilities/StringUtils.h" #include "ThePEG/Utilities/Throw.h" #include "ThePEG/PDT/DecayMode.h" #ifdef ThePEG_TEMPLATES_IN_CC_FILE #include "BaseRepository.tcc" #endif using namespace ThePEG; ostream *& BaseRepository::coutp() { static ostream * theCout = &std::cout; return theCout; } ostream *& BaseRepository::cerrp() { static ostream * theCerr = &std::cerr; return theCerr; } ostream *& BaseRepository::clogp() { static ostream * theClog = &std::clog; return theClog; } bool & BaseRepository::updating() { static bool theBool = false; return theBool; } ObjectMap & BaseRepository::objects() { static ObjectMap theObjectMap; return theObjectMap; } ObjectSet & BaseRepository::allObjects() { static ObjectSet theObjectSet; return theObjectSet; } BaseRepository::TypeInterfaceMap & BaseRepository::interfaces() { static TypeInterfaceMap theInterfaceMap; return theInterfaceMap; } BaseRepository::TypeDocumentationMap & BaseRepository::documentations() { static TypeDocumentationMap theDocumentationMap; return theDocumentationMap; } BaseRepository::DirectorySet & BaseRepository::directories() { static DirectorySet theDirectories = {"/"}; return theDirectories; } vector & BaseRepository::globalLibraries() { static vector theGlobalLibraries; return theGlobalLibraries; } stack & BaseRepository::currentReadDirStack() { static stack theCurrentReadDirStack; if ( theCurrentReadDirStack.empty() ) theCurrentReadDirStack.push(""); return theCurrentReadDirStack; } vector & BaseRepository::readDirs() { // macro is passed in from -D compile flag static vector theReadDirs(1, THEPEG_PKGDATADIR); return theReadDirs; } const vector & BaseRepository::listReadDirs() { return BaseRepository::readDirs(); } void BaseRepository::prependReadDir(string dir) { readDirs().insert(readDirs().begin(), dir); } void BaseRepository::prependReadDir(const std::vector& dirs) { readDirs().insert(readDirs().begin(), dirs.begin(), dirs.end()); } void BaseRepository::appendReadDir(string dir) { readDirs().push_back(dir); } void BaseRepository::appendReadDir(const std::vector& dirs) { readDirs().insert(readDirs().end(), dirs.begin(), dirs.end()); } BaseRepository::StringVector & BaseRepository::directoryStack() { static StringVector theDirectoryStack(1, "/"); return theDirectoryStack; } void BaseRepository::Register(const InterfaceBase & ib, const type_info & i) { const ClassDescriptionBase * db = DescriptionList::find(i); if ( db ) interfaces()[db].insert(&ib); } void BaseRepository:: Register(const ClassDocumentationBase & cd, const type_info & i) { const ClassDescriptionBase * db = DescriptionList::find(i); if ( db ) documentations()[db] = &cd; } void BaseRepository::Register(IBPtr ip, string newName) { DirectoryAppend(newName); ip->name(newName); Register(ip); } void BaseRepository::Register(IBPtr ip) { if ( !ip || member(allObjects(), ip) ) return; if ( member(objects(), ip->fullName()) ) throw RepoNameExistsException(ip->fullName()); objects()[ip->fullName()] = ip; allObjects().insert(ip); ip->clear(); ip->update(); ip->touch(); } void BaseRepository::DirectoryAppend(string & name) { if ( name == "." ) name = directoryStack().back(); if ( name[0] != '/' ) name = directoryStack().back() + name; } void BaseRepository::CreateDirectory(string name) { DirectoryAppend(name); if ( name[name.size()-1] != '/' ) name += "/"; if ( member(directories(), name) ) return; directories().insert(name); name = name.substr(0, name.size() - 1); name = name.substr(0, name.rfind('/')); if ( name.size() ) CreateDirectory(name); } void BaseRepository::CheckObjectDirectory(string name) { if ( name[name.size() - 1] != '/' ) name = name.substr(0, name.rfind('/') + 1); CheckDirectory(name); } void BaseRepository::CheckDirectory(string name) { DirectoryAppend(name); if ( name[name.size()-1] != '/' ) name += "/"; if ( member(directories(), name) ) return; throw RepositoryNoDirectory(name); } void BaseRepository::ChangeDirectory(string name) { DirectoryAppend(name); if ( name[name.size()-1] != '/' ) name += "/"; if ( member(directories(), name) ) { directoryStack().back() = name; return; } throw RepositoryNoDirectory(name); } void BaseRepository::PushDirectory(string name) { DirectoryAppend(name); if ( name[name.size()-1] != '/' ) name += "/"; if ( member(directories(), name) ) { directoryStack().push_back(name); return; } throw RepositoryNoDirectory(name); } void BaseRepository::PopDirectory() { if ( directoryStack().size() > 1 ) directoryStack().pop_back(); } IBPtr BaseRepository::GetPointer(string name) { ObjectMap::iterator it = objects().find(name); return it == objects().end()? IBPtr(): it->second; } IVector BaseRepository::SearchDirectory(string name, string className) { IVector ret; DirectoryAppend(name); const ClassDescriptionBase * cdb = 0; if ( className.size() ) { cdb = DescriptionList::find(className); if ( !cdb ) return ret; } if ( name[name.size()-1] != '/' ) name += "/"; string::size_type size = name.size(); for ( ObjectMap::const_iterator i = objects().begin(); i != objects().end(); ++i ) { - if ( cdb && !DescriptionList::find(typeid(*(i->second)))->isA(*cdb) ) + const auto & tmp=(*(i->second)); + if ( cdb && !DescriptionList::find(typeid(tmp))->isA(*cdb) ) continue; if ( i->first.substr(0, size) == name ) ret.push_back(i->second); } return ret; } IVector BaseRepository::GetObjectsReferringTo(IBPtr obj) { IVector ret; for ( ObjectMap::const_iterator i = objects().begin(); i != objects().end(); ++i ) { if ( obj == i->second ) continue; IVector ov = DirectReferences(i->second); if ( member(ov, obj) ) ret.push_back(i->second); } return ret; } IVector BaseRepository::DirectReferences(IBPtr obj) { IVector ov = obj->getReferences(); - InterfaceMap interfaceMap = getInterfaces(typeid(*obj)); + const auto & tmp=*obj; + InterfaceMap interfaceMap = getInterfaces(typeid(tmp)); for ( InterfaceMap::iterator iit = interfaceMap.begin(); iit != interfaceMap.end(); ++iit ) { IVector ovi = iit->second->getReferences(*obj); ov.insert(ov.end(), ovi.begin(), ovi.end()); } return ov; } void BaseRepository:: addReferences(tIBPtr obj, ObjectSet & refs) { if ( !obj ) return; refs.insert(obj); IVector ov = obj->getReferences(); for ( IVector::const_iterator it = ov.begin(); it != ov.end(); ++it ) if ( !member(refs, *it) ) addReferences(*it, refs); - InterfaceMap interfaceMap = getInterfaces(typeid(*obj)); + const auto & tmp=*obj; + InterfaceMap interfaceMap = getInterfaces(typeid(tmp)); for ( InterfaceMap::iterator iit = interfaceMap.begin(); iit != interfaceMap.end(); ++iit ) { IVector ov = iit->second->getReferences(*obj); for ( IVector::const_iterator it = ov.begin(); it != ov.end(); ++it ) if ( !member(refs, *it) ) addReferences(*it, refs); } } void BaseRepository:: addInterfaces(const ClassDescriptionBase & db, InterfaceMap & interfaceMap, bool all) { for ( ClassDescriptionBase::DescriptionVector::const_iterator it = db.descriptions().begin(); it != db.descriptions().end(); ++it ) if ( *it ) addInterfaces(**it, interfaceMap, all); TypeInterfaceMap::const_iterator cit = interfaces().find(&db); if ( cit == interfaces().end() ) return; for ( InterfaceSet::const_iterator iit = (cit->second).begin(); iit != (cit->second).end(); ++iit ) { string n = (**iit).name(); while ( all && member(interfaceMap, n) ) n = "+" + n; interfaceMap[n] = *iit; } } InterfaceMap BaseRepository::getInterfaces(const type_info & ti, bool all) { InterfaceMap interfaceMap; const ClassDescriptionBase * db = DescriptionList::find(ti); if ( !db ) return interfaceMap; addInterfaces(*db, interfaceMap, all); return interfaceMap; } void BaseRepository:: rebind(InterfacedBase & i, const TranslationMap & trans, const IVector & defaults) { InterfaceMap interfaceMap = getInterfaces(typeid(i), true); for ( InterfaceMap::iterator iit = interfaceMap.begin(); iit != interfaceMap.end(); ++iit ) iit->second->rebind(i, trans, defaults); i.rebind(trans); } void BaseRepository::update() { - for_each(allObjects(), mem_fun(&InterfacedBase::update)); + for_each(allObjects(), std::mem_fn(&InterfacedBase::update)); clearAll(allObjects()); } template bool overlap(const Set1 & s1, const Set2 & s2) { typename Set1::const_iterator i1 = s1.begin(); typename Set2::const_iterator i2 = s2.begin(); while ( i1 != s1.end() && i2 != s2.end() ) { if ( *i1 == *i2 ) return true; if ( *i1 < *i2 ) { i1 = s1.lower_bound(*i2); if ( *i1 == *i2 ) return true; ++i1; } else { i2 = s2.lower_bound(*i1); if ( *i1 == *i2 ) return true; ++i2; } } return false; } void BaseRepository::remove(tIBPtr ip) { ObjectMap::iterator it = objects().find(ip->fullName()); if ( it == objects().end() || ip != it->second ) return; objects().erase(it); allObjects().erase(ip); } string BaseRepository::remove(const ObjectSet & rmset) { ObjectSet refset; for ( ObjectSet::const_iterator oi = rmset.begin(); oi != rmset.end(); ++oi ) { IVector ov = GetObjectsReferringTo(*oi); refset.insert(ov.begin(), ov.end()); } for ( ObjectSet::iterator oi = rmset.begin(); oi != rmset.end(); ++oi ) refset.erase(*oi); if ( refset.empty() ) { for ( ObjectSet::iterator oi = rmset.begin(); oi != rmset.end(); ++oi ) remove(*oi); return ""; } string ret = "Error: cannot remove the objects because the following " "objects refers to some of them:\n"; for ( ObjectSet::iterator oi = refset.begin(); oi != refset.end(); ++oi ) ret += (**oi).fullName() + "\n"; return ret; } void BaseRepository::rename(tIBPtr ip, string newName) { ObjectSet::iterator it = allObjects().find(ip); if ( it == allObjects().end() ) { Register(ip, newName); return; } ObjectMap::iterator mit = objects().find(ip->fullName()); if ( mit == objects().end() || mit->second != ip ) throw RepoNameException(ip->fullName()); objects().erase(mit); ip->name(newName); while ( member(objects(), ip->fullName()) ) ip->name(ip->fullName() + "#"); objects()[ip->fullName()] = ip; } const InterfaceBase * BaseRepository::FindInterface(IBPtr ip, string name) { - InterfaceMap imap = getInterfaces(typeid(*ip), false); + const auto & tmp=*ip; + InterfaceMap imap = getInterfaces(typeid(tmp), false); InterfaceMap::iterator it = imap.find(name); return it == imap.end()? 0: it->second; } const ClassDocumentationBase * BaseRepository::getDocumentation(tcIBPtr ip) { + const auto & tmp=*ip; TypeDocumentationMap::const_iterator cdoc = - documentations().find(DescriptionList::find(typeid(*ip))); + documentations().find(DescriptionList::find(typeid(tmp))); return cdoc != documentations().end()? cdoc->second: 0; } string BaseRepository::getModelDescription(tcIBPtr ip) { const ClassDocumentationBase * cd = getDocumentation(ip); return cd? cd->modelDescription(): string(""); } string BaseRepository::getModelReferences(tcIBPtr ip) { const ClassDocumentationBase * cd = getDocumentation(ip); return cd? cd->modelReferences(): string(""); } IBPtr BaseRepository::TraceObject(string path) { DirectoryAppend(path); string::size_type colon = path.find(':'); IBPtr ip = GetPointer(path.substr(0, colon)); if ( !ip ) { // Do special check if this is a decay mode. string name = path.substr(0, colon); string::size_type slash = name.rfind('/'); if ( slash != string::npos ) name = name.substr(slash + 1); if ( name.find("->") != string::npos && name[name.length() - 1] == ';' ) { vector save; DMPtr dm = DecayMode::constructDecayMode(name, &save); if ( dm ) ip = dynamic_ptr_cast(GetPointer(path.substr(0, slash + 1) + dm->tag())); if ( ip ) Throw() << "Warning: rewriting DecayMode name '" << path.substr(0, colon).substr(slash + 1) << "' to '" << ip->name() << Exception::warning; } } while ( colon != string::npos ) { if ( !ip ) throw RepositoryNotFound(path); path = path.substr(colon+1); colon = path.find(':'); string::size_type bra = path.find('['); const InterfaceBase * ifb = FindInterface(ip, path.substr(0, min(colon, bra))); const ReferenceBase * rb = dynamic_cast(ifb); if ( rb ) { ip = rb->get(*ip); continue; } const RefVectorBase * rvb = dynamic_cast(ifb); if ( rvb ) { unsigned int place = 0; if ( bra < colon ) { string::size_type ket = path.find(']'); place = atoi(path.substr(bra + 1,ket - bra - 1).c_str()); } IVector iv = rvb->get(*ip); if ( place >= iv.size() ) throw RepositoryNotFound(path); ip = iv[place]; continue; } const CommandBase * cb = dynamic_cast(ifb); if ( cb ) { string::size_type ket = path.find(']'); string newobj = cb->cmd(*ip, path.substr(bra + 1,ket - bra - 1)); ip = GetPointer(newobj); continue; } throw RepositoryNotFound(path); } if ( !ip ) throw RepositoryNotFound(path); return ip; } IBPtr BaseRepository::getObjectFromNoun(string noun) { string::size_type colon = noun.rfind(':'); return TraceObject(noun.substr(0, colon)); } string BaseRepository::getInterfaceFromNoun(string noun) { string::size_type colon = noun.rfind(':'); string interface = noun.substr(colon+1); string::size_type bra = interface.find('['); if ( bra != string::npos ) return interface.substr(0, bra); else return interface; } string BaseRepository::getPosArgFromNoun(string noun) { string::size_type colon = noun.rfind(':'); string interface = noun.substr(colon+1); string::size_type bra = interface.find('['); if ( bra != string::npos ) { string::size_type ket = interface.find(']'); return interface.substr(bra + 1,ket - bra - 1); } return ""; } string BaseRepository:: GetInterfacedBaseClasses(const ClassDescriptionBase * cdb) { if ( !cdb || cdb->name() == "ThePEG::Interfaced" || cdb->name() == "ThePEG::InterfacedBase" ) return ""; string ret = cdb->name() + "\n"; for ( int i = 0, N = cdb->descriptions().size(); i < N; ++i ) ret += GetInterfacedBaseClasses(cdb->descriptions()[i]); return ret; } struct InterfaceOrder { bool operator()(const InterfaceBase * x, const InterfaceBase * y) const { return x->rank() > y->rank() || ( x->rank() == y->rank() && x->name() < y->name() ); } }; void BaseRepository::readSetup(tIBPtr ip, istream & is) { ip->setup(is); } string BaseRepository::exec(string command, ostream &) { string verb = StringUtils::car(command); command = StringUtils::cdr(command); if ( verb.empty() || verb[0] == '#' ) return ""; try { if ( verb == "DISABLEREADONLY" ) { InterfaceBase::NoReadOnly = true; return ""; } if ( verb == "ENABLEREADONLY" ) { InterfaceBase::NoReadOnly = false; return ""; } if ( verb == "cd" || verb == "pushd" || verb == "mkdir") { string dir = StringUtils::car(command); if ( verb == "cd" ) ChangeDirectory(dir); else if ( verb == "pushd" ) PushDirectory(dir); else CreateDirectory(dir); return ""; } if ( verb == "popd" ) { PopDirectory(); return ""; } if ( verb == "pwd" ) return directoryStack().back(); if ( verb == "dirs" ) { string ret; for ( StringVector::reverse_iterator it = directoryStack().rbegin(); it != directoryStack().rend(); ++it ) ret += *it; return ret; } if ( verb == "cp" || verb == "mv" ) { string oldname = StringUtils::car(command); DirectoryAppend(oldname); IBPtr obj = GetPointer(oldname); if ( !obj ) return "Error: No object named '" + oldname + "' available."; command = StringUtils::cdr(command); string newname = StringUtils::car(command); DirectoryAppend(newname); if ( newname[newname.size() - 1] == '/' ) newname += obj->name(); if ( verb == "cp" ) obj = obj->fullclone(); rename(obj, newname); return ""; } if ( verb == "check" ) { string name = StringUtils::car(command); if ( directories().find(name) != directories().end() ) return name; if ( objects().find(name) != objects().end() ) return name; return "Not found"; } if ( verb == "ls" ) { string className; string dir = StringUtils::car(command); if ( dir.size() ) { PushDirectory(dir); command = StringUtils::cdr(command); className = StringUtils::car(command); } string ret; string thisdir = directoryStack().back(); for ( DirectorySet::iterator it = directories().begin(); it != directories().end(); ++it ) { string d = *it; if ( d.size() <= thisdir.size() ) continue; string d0 = d.substr(0, thisdir.size()); string d1 = d.substr(thisdir.size()); if ( d0 == thisdir && d1.find('/') == d1.size() - 1 ) { if ( className.size() && SearchDirectory(d, className).empty() ) continue; ret += (dir.size()? d: d1) + "\n"; } } for ( ObjectMap::iterator it = objects().begin(); it != objects().end(); ++it ) { if ( className.size() ) { const ClassDescriptionBase * cdb = DescriptionList::find(className); + const auto & tmp=*(it->second); if ( cdb && - !DescriptionList::find(typeid(*(it->second)))->isA(*cdb) ) + !DescriptionList::find(typeid(tmp))->isA(*cdb) ) continue; } if ( thisdir + it->second->name() == it->first ) ret += (dir.size()? it->first: it->second->name()) + '\n'; } if ( dir.size() ) PopDirectory(); return ret; } if ( verb == "library" ) { string library = StringUtils::car(command); if ( library.empty() ) return "Error: No library specified."; if ( !DynamicLoader::load(library) ) return "Error: Could not load library " + library + "\n - " + DynamicLoader::lastErrorMessage; return ""; } if ( verb == "globallibrary" ) { string library = StringUtils::car(command); if ( library.empty() ) return "Error: No library specified."; if ( !DynamicLoader::load(library) ) return "Error: Could not load library " + library + "\n - " + DynamicLoader::lastErrorMessage; globalLibraries().push_back(library); return ""; } if ( verb == "rmgloballibrary" ) { string library = StringUtils::car(command); if ( library.empty() ) return "Error: No library specified."; vector::iterator it; while ( (it = find(globalLibraries(), library)) != globalLibraries().end() ) globalLibraries().erase(it); return ""; } if ( verb == "appendpath" ) { string path = StringUtils::car(command); if ( !path.empty() ) DynamicLoader::appendPath(path); return ""; } if ( verb == "lspaths" ) { string paths; for ( int i = 0, N = DynamicLoader::allPaths().size(); i < N; ++i ) paths += DynamicLoader::allPaths()[i] + "\n"; return paths; } if ( verb == "prependpath" ) { string path = StringUtils::car(command); if ( !path.empty() ) DynamicLoader::prependPath(path); return ""; } if ( verb == "create" ) { string className = StringUtils::car(command); command = StringUtils::cdr(command); string name = StringUtils::car(command); const ClassDescriptionBase * db = DescriptionList::find(className); command = StringUtils::cdr(command); while ( !db && command.length() ) { string library = StringUtils::car(command); command = StringUtils::cdr(command); DynamicLoader::load(library); db = DescriptionList::find(className); } if ( !db ) { string msg = "Error: " + className + ": No such class found."; if ( !DynamicLoader::lastErrorMessage.empty() ) msg += "\nerror message from dynamic loader:\n" + DynamicLoader::lastErrorMessage; return msg; } IBPtr obj = dynamic_ptr_cast(db->create()); if ( !obj ) return "Error: Could not create object of class "+className; if ( name.empty() ) return "Error: No name specified."; Register(obj, name); return ""; } if ( verb == "setup" ) { string name = StringUtils::car(command); DirectoryAppend(name); IBPtr obj = GetPointer(name); if ( !obj ) return "Error: Could not find object named " + name; istringstream is(StringUtils::cdr(command)); readSetup(obj, is); return ""; } if ( verb == "rm" ) { ObjectSet rmset; while ( !command.empty() ) { string name = StringUtils::car(command); DirectoryAppend(name); IBPtr obj = GetPointer(name); if ( !obj ) return "Error: Could not find object named " + name; rmset.insert(obj); command = StringUtils::cdr(command); } return remove(rmset); } if ( verb == "rmdir" || verb == "rrmdir" ) { string dir = StringUtils::car(command); DirectoryAppend(dir); if ( dir[dir.size() - 1] != '/' ) dir += '/'; if ( !member(directories(), dir) ) return verb == "rmdir"? "Error: No such directory.": ""; IVector ov = SearchDirectory(dir); if ( ov.size() && verb == "rmdir" ) return "Error: Cannot remove a non-empty directory. " "(Use rrmdir do remove all object and subdirectories.)"; ObjectSet rmset(ov.begin(), ov.end()); string ret = remove(rmset); if ( !ret.empty() ) return ret; StringVector dirs(directories().begin(), directories().end()); for ( int i = 0, N = dirs.size(); i < N; ++ i ) if ( dirs[i].substr(0, dir.size()) == dir ) directories().erase(dirs[i]); for ( int i = 0, N = directoryStack().size(); i < N; ++i ) if ( directoryStack()[i].substr(0, dir.size()) == dir ) directoryStack()[i] = '/'; return ""; } if ( verb == "rcp" ) { string name = StringUtils::car(command); DirectoryAppend(name); string newName = StringUtils::car(StringUtils::cdr(command)); if ( newName.empty() ) return "Error: No destination directory specified."; DirectoryAppend(newName); CreateDirectory(newName); if ( newName[newName.size() - 1] != '/' ) newName += '/'; IBPtr obj = GetPointer(name); if ( name[name.size() - 1] != '/' ) name += '/'; IVector ov = SearchDirectory(name); ov.push_back(obj); if ( ov.empty() ) return "Error: No such object or directory."; ObjectSet toclone; for ( IVector::iterator i = ov.begin(); i != ov.end(); ++i ) { toclone.insert(*i); addReferences(*i, toclone); } for ( ObjectSet::iterator i = toclone.begin(); i != toclone.end(); ++i ) Register((**i).clone(), newName + (**i).name()); return ""; } if ( verb == "rebind" ) { // For all objects in the repository, replace any references to // the first object given with references to the second // object. The two objects will change names IBPtr ip1 = TraceObject(StringUtils::car(command)); string newname = StringUtils::car(StringUtils::cdr(command)); DirectoryAppend(newname); IBPtr ip2 = GetPointer(newname); if ( !ip2 ) { ip2 = ip1->fullclone(); rename(ip2, newname); } TranslationMap trans; trans[ip1] = ip2; IVector objs = GetObjectsReferringTo(ip1); for ( int i = 0, N = objs.size(); i < N; ++i ) rebind(*objs[i], trans, IVector()); } if ( verb == "doxygendump" ) { string spacename = StringUtils::car(command); command = StringUtils::cdr(command); string filename = StringUtils::car(command); ofstream os(filename.c_str()); for ( TypeDocumentationMap::const_iterator it = documentations().begin(); it != documentations().end(); ++it ) { const ClassDescriptionBase & db = *(it->first); string classname = db.name(); if ( classname.substr(0, spacename.length()) != spacename ) continue; string briefname = classname.substr(spacename.length()); os << "/** \\page " << briefname << "Interfaces " << "Interfaces defined for the " << classname << " class.\n\n" << "\\par Brief class description:\n"; string doc = it->second->documentation(); if ( doc.substr(0,25) == "There is no documentation" ) os << "See " << classname << "\n\n"; else os << doc << "
See also " << classname << "\n\n"; TypeInterfaceMap::const_iterator isit = interfaces().find(it->first); if ( isit == interfaces().end() || isit->second.empty() ) { os << "There are no interfaces declared for this class.\n\n"; } else { const InterfaceSet & ints = isit->second; for ( InterfaceSet::const_iterator iit = ints.begin(); iit != ints.end(); ++iit ) (**iit).doxygenDescription(os); } string baserefs = ""; int nbases = 0; for ( int ib = 0, N = db.descriptions().size(); ib < N; ++ib ) { if ( documentations().find(db.descriptions()[ib]) == documentations().end() ) continue; const ClassDescriptionBase & bdb = *db.descriptions()[ib]; if ( nbases ) baserefs += " and "; string briefname = bdb.name().substr(bdb.name().rfind("::") + 2); baserefs += "\\ref " + briefname + "Interfaces \"" + bdb.name() + "\""; ++nbases; } if ( nbases == 1 ) os << "
There may be interfaces inherited from the " << baserefs << " class."; else if ( nbases > 1 ) os << "
There may be interfaces inherited from the " << "following classes: " << baserefs << "."; os << "\n\n*/\n\n"; } return ""; } if ( verb == "mset" || verb == "msetdef" || verb == "minsert" || verb == "mdo" || verb == "mget" || verb == "mdef" || verb == "mmin" || verb == "mmax" || verb == "merase" || verb == "msend" ) { if ( verb == "msend" ) verb = "mdo"; string dir = StringUtils::car(command); command = StringUtils::cdr(command); string className = StringUtils::car(command); command = StringUtils::cdr(command); string interface = StringUtils::car(command); string arguments = StringUtils::cdr(command); string::size_type bra = interface.find('['); if ( bra != string::npos ) { string::size_type ket = interface.find(']'); arguments = interface.substr(bra + 1,ket - bra - 1) + " " + arguments; interface = interface.substr(0, bra); } IVector ov = SearchDirectory(dir, className); if ( ov.empty() ) return "Error: no matching objects found."; string ret; verb = verb.substr(1); for ( IVector::size_type i = 0; i < ov.size(); ++i ) { const InterfaceBase * ifb = FindInterface(ov[i], interface); if ( !ifb ) continue; string mess = ifb->exec(*ov[i], verb, arguments); if ( !mess.empty() ) ret += ov[i]->fullName() + ": " + mess + "\n"; } return ret.substr(0, ret.size() - 1); } if ( verb == "set" || verb == "setdef" || verb == "insert" || verb == "do" || verb == "get" || verb == "def" || verb == "min" || verb == "max" || verb == "describe" || verb == "fulldescribe" || verb == "erase" || verb == "clear" || verb == "send" || verb == "newdef" ) { if ( verb == "send" ) verb = "do"; if ( verb == "newdef" && !InterfaceBase::NoReadOnly ) return "Error: The default value of an interface is a read-only " "entity. Use the command 'DISABLEREADONLY' to override."; string noun = StringUtils::car(command); string arguments = getPosArgFromNoun(noun) + " " + StringUtils::cdr(command); IBPtr ip = getObjectFromNoun(noun); const InterfaceBase * ifb = FindInterface(ip, getInterfaceFromNoun(noun)); if ( !ifb && verb != "describe" && verb != "fulldescribe" ) { string ret = "Error: The interface '" + noun + "' was not found.\n"; ret += "Valid interfaces:\n"; - InterfaceMap imap = getInterfaces(typeid(*ip)); + const auto & tmp=*ip; + InterfaceMap imap = getInterfaces(typeid(tmp)); for ( InterfaceMap::iterator it = imap.begin(); it != imap.end(); ++it ) ret += "* " + it->second->name() + "\n"; return ret; } if ( verb == "describe" ) { if ( ifb ) return ifb->description(); - const ClassDescriptionBase * cd = DescriptionList::find(typeid(*ip)); + const auto & tmp=*ip; + const ClassDescriptionBase * cd = DescriptionList::find(typeid(tmp)); string ret = "Object '" + ip->name() + "' of class '" + cd->name() + "':\n"; TypeDocumentationMap::const_iterator cdoc = documentations().find(cd); if ( cdoc != documentations().end() ) ret += cdoc->second->documentation() + "\n"; ret +="Interfaces:\n"; - InterfaceMap imap = getInterfaces(typeid(*ip)); + InterfaceMap imap = getInterfaces(typeid(tmp)); for ( InterfaceMap::iterator it = imap.begin(); it != imap.end(); ++it ) ret += "* " + it->second->name() + "\n"; return ret; } else if ( verb == "fulldescribe" ) { if ( ifb ) return ifb->fullDescription(*ip); ostringstream ret; - const ClassDescriptionBase * cd = DescriptionList::find(typeid(*ip)); + const auto & tmp=*ip; + const ClassDescriptionBase * cd = DescriptionList::find(typeid(tmp)); TypeDocumentationMap::const_iterator cdoc = documentations().find(cd); ret << ip->fullName() << endl << cd->name() << endl; if ( cdoc != documentations().end() ) ret << cdoc->second->documentation() << endl; ret << "Interfaces:" << endl; - InterfaceMap imap = getInterfaces(typeid(*ip)); + InterfaceMap imap = getInterfaces(typeid(tmp)); typedef set InterfaceSet; InterfaceSet iset; for ( InterfaceMap::iterator it = imap.begin(); it != imap.end(); ++it ) iset.insert(it->second); double rank = 1.0; for ( InterfaceSet::iterator it = iset.begin(); it != iset.end(); ++it ) { if ( rank >= 0.0 && (**it).rank() < 0.0 ) ret << "0" << endl; rank = (**it).rank(); ret << (**it).type() << " " << (**it).name() << endl; } return ret.str(); } else return ifb->exec(*ip, verb, arguments); } if ( verb == "baseclasses" ) { string className = StringUtils::car(command); const ClassDescriptionBase * cdb = 0; if ( className.size() ) { cdb = DescriptionList::find(className); if ( !cdb ) return "Error: no class '" + className + "' found."; } return GetInterfacedBaseClasses(cdb); } if ( verb == "describeclass" ) { string className = StringUtils::car(command); const ClassDescriptionBase * cdb = 0; if ( className.size() ) { cdb = DescriptionList::find(className); if ( !cdb ) return "Error: no class '" + className + "' found."; } TypeDocumentationMap::const_iterator cdoc = documentations().find(cdb); if ( cdoc != documentations().end() ) return cdoc->second->documentation() + "\n"; else return ""; } if ( verb == "lsclass" ) { string className = StringUtils::car(command); const ClassDescriptionBase * cdb = 0; if ( className.size() ) { cdb = DescriptionList::find(className); if ( !cdb ) return "Error: no class '" + className + "' found."; } vector classes; if ( cdb && !cdb->abstract() ) classes.push_back(cdb); for ( DescriptionList::DescriptionMap::const_iterator it = DescriptionList::all().begin(); it != DescriptionList::all().end(); ++it ) { if ( it->second == cdb || it->second->abstract() ) continue; if ( cdb && !it->second->isA(*cdb) ) continue; classes.push_back(it->second); } if ( classes.empty() ) return "Error: no classes found."; string ret; for ( int i = 0, N = classes.size(); i < N; ++i ) ret += classes[i]->name() + "\n"; return ret; } } catch (const Exception & e) { e.handle(); return "Error: " + e.message(); } return "Error: Unrecognized command '" + verb + "'."; } BadClassClone::BadClassClone(const InterfacedBase & o) { theMessage << "Could not clone the object '" << o.name() << "' of class '" << TypeInfo::name(o) << "' because the class does not" << " implement a working 'clone' method."; severity(abortnow); } BadClone::BadClone(const InterfacedBase & o) { theMessage << "Could not clone the object '" << o.name() << "' of class '" << TypeInfo::name(o) << "' because the clone method threw an unknown exception."; severity(abortnow); } RepoNameException::RepoNameException(string name) { theMessage << "The object '" << name << "' is present in the Repository but " << "under a different name. This means that the name of the " << "object has been illegally changed outside of the Repository."; severity(abortnow); } RepoNameExistsException::RepoNameExistsException(string name) { theMessage << "The object '" << name << "' was not created as another object with that name already exists."; severity(warning); } RepositoryNoDirectory::RepositoryNoDirectory(string name) { theMessage << "The directory '" << name << "' does not exist."; severity(warning); } RepositoryNotFound::RepositoryNotFound(string name) { theMessage << "There was no object named '" << name << "' in the repository."; severity(warning); } RepositoryClassMisMatch:: RepositoryClassMisMatch(const InterfacedBase & o, string name) { theMessage << "The requested object '" << o.fullName() << "' was not of the " << "specified type (" << name << ")."; severity(warning); } diff --git a/Repository/BaseRepository.h b/Repository/BaseRepository.h --- a/Repository/BaseRepository.h +++ b/Repository/BaseRepository.h @@ -1,583 +1,583 @@ // -*- C++ -*- // // BaseRepository.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_BaseRepository_H #define ThePEG_BaseRepository_H // This is the declaration of the BaseRepository class. #include "ThePEG/Config/ThePEG.h" #include "BaseRepository.xh" #include "ThePEG/Interface/InterfaceBase.fh" #include "ThePEG/Interface/ClassDocumentation.fh" #include "ThePEG/Interface/InterfacedBase.h" #include "ThePEG/Utilities/ClassDescription.fh" namespace ThePEG { /** * BaseRepository is a purely static class which keeps a set of * InterfacedBase objects indexed by their name. The objects and their * names are divided up in a tree-like structure inspired by the Unix * file system. * * The InterfacedBase objects may be manipulated using InterfaceBase * objects. This may be done directly or via a simple command * interface using the exec() method. * * RepositoryBase is closely related to the Repository sub-class. The * division may seem unnecessary, but the idea is that BaseRepository * is a general repository for administrating and manipulating a set * of InterfacedBase objects, while the Repository adds on utilites * which are special to ThePEG where the objects are Interfaced (a * sub-class of InterfacedBase). * * @see Repository * @see InterfacedBase * @see InterfaceBase * @see Interfaced * */ class BaseRepository { public: /** A set of strings. */ typedef StringSet DirectorySet; /** A vector of character strings. */ typedef vector StringVector; /** A set of pointers to InterfaceBase objects. */ typedef set InterfaceSet; /** A map of sets of IterfaceBase objects indexed by pointers to ClassDescriptionBase objects. */ typedef map TypeInterfaceMap; /** A map of ClassDocumentationBase objects indexed by pointers to ClassDescriptionBase objects. */ typedef map TypeDocumentationMap; public: /** * Interpret the command in \a cmd and return possible * messages. This is the main function for the command-line * interface. The syntax is described elsewhere. The ostream * argument is currently unused. */ static string exec(string cmd, ostream &); /** @name Functions for adding and deleting objects and interfaces. */ //@{ /** * Register an interface. This is called automatically in the * InterfaceBase constructor and should never be called explicitly. */ static void Register(const InterfaceBase &, const type_info &); /** * Register a class documentation. This is called automatically in * the ClassDocumentationBase constructor and should never be called * explicitly. */ static void Register(const ClassDocumentationBase &, const type_info &); /** * Register a new object using the its current name. If the object * is already in the repository, nothing happens. If another object * already exists with the same name, the new object will have * #'s appended to its name to make it unique. */ static void Register(IBPtr); /** * Register a new object giving it a new \a name. If the object is * already in the repository, nothing happens. If another object * already exists with the same name, the new object will have * #'s appended to its name to make it unique. */ static void Register(IBPtr, string name); /** * Remove the given object from the repository. If the object was * not present nothing will happen. */ static void remove(tIBPtr); /** * Remove objects. Remove the objects in \a rmset if there are no * other objects in the repository referring to them, otherwise * return an error message and the names of the objects refering to * them separated by new-line characters. */ static string remove(const ObjectSet & rmset); /** * Rename a given \a object. Syntacticly the same as * remove(object); Register(object, newName);. */ static void rename(tIBPtr object, string newName); //@} /** @name Access the directory stack. */ //@{ /** * Create a new directory with the given name. If the given name * starts with a / the name is assumed to be an absolute * path, otherwise it is assumed to be a path relative to the * current directory. */ static void CreateDirectory(string); /** * Check if directory exixts. Check if the name given as argument * corresponds to an existing directory. If the argument string does * not end in a / it is assumed to be the name of an * object in a directory, and only the directory part of the name is * checked. If the given name starts with a / the name * is assumed to be an absolute path, otherwise it is assumed to be * a path relative to the current directory. * * @throws RepositoryNoDirectory if the correspinding directory is * non-existent. */ static void CheckObjectDirectory(string); /** * Check if directory exixts. Check if the name given as argument * corresponds to an existing directory. If the given name starts * with a / the name is assumed to be an absolute path, * otherwise it is assumed to be a path relative to the current * directory. * * @throws RepositoryNoDirectory if the correspinding directory is * non-existent. */ static void CheckDirectory(string); /** * Return the absolute path. If the given name starts with a * / the name is assumed to be an absolute path already, * otherwise it is assumed to be a path relative to the current * directory, and the absolute path is constructed. */ static void DirectoryAppend(string &); /** * Set the current directory to \a name. \a name can be aither a * relative or absolute path. The new directory replaces the * previous current directory on the directory stack. * * @throws RepositoryNoDirectory if the directory is non-existent. */ static void ChangeDirectory(string name); /** * Set the current directory to \a name. \a name can be aither a * relative or absolute path. The new directory is pushed onto the * directory stack. * * @throws RepositoryNoDirectory if the directory is non-existent. */ static void PushDirectory(string name); /** * Pop the directory stack. Leave the current directory and set the * directory which is on top of the popped directory stack. */ static void PopDirectory(); /** * A list of all globally loaded libraries. */ static vector & globalLibraries(); //@} /** @name Information on where to read input files. */ //@{ protected: /** * The stack of directories used by the "read" command. */ static stack & currentReadDirStack(); /** * List of directories to search for files for the "read" command. */ static vector & readDirs(); public: /** * Access to list of directories to search for files for the "read" command. */ static const vector & listReadDirs(); /** * Add a directory to readDirs(). */ static void prependReadDir(string); /** * Add a string vector with directories to readDirs(). */ static void prependReadDir(const std::vector& dirs); /** * Add a directory to readDirs(). */ static void appendReadDir(string); /** * Add a string vector with directories to readDirs(). */ static void appendReadDir(const std::vector& dirs); //@} /** @name Access objects in the repository. */ //@{ /** * Return a reference counted pointer to the given object. This * currently not needed when ThePEG is used with the * ThePEG::Pointer::RCPtr class of pointers. */ template static typename Ptr::pointer GetPtr(const T &); /** * Return a pointer of the specified type to an object with the * given name. If such an object does not exist, GetPtr will return * a null pointer. */ template static PtrType GetPtr(string); /** * Return a pointer of the specified type to an object with the * given name. If such an object does not exist an exception will be * thrown. * @throws RepositoryNotFound if the object was not found. * @throws RepositoryClassMisMatch if the object exists but is of * the wrong class. */ template static PtrType GetObject(string); /** * Return a pointer to an object with the given name or null if no * such object exists. */ static IBPtr GetPointer(string); /** * Return all objects in the directory \a name. Optionally only return * objects of class \a className or of a sub-class thereof. */ static IVector SearchDirectory(string name, string className = ""); /** * Find an object. If the \a name does not begin with '/', the * current directory is prepended. If the string is on the form * object:interface (or * object:interface[i]) and interface * corresponds to an Reference (or RefVector) interface, the * corresponding referenced object is returned. (also * object:interface:interface is allowed etc.) */ static IBPtr TraceObject(string name); /** * Return a string containing the name of the given class * description and its base classes, one on each line. */ static string GetInterfacedBaseClasses(const ClassDescriptionBase * cdb); /** * Get an object. Decompose a string of the form * object:interface or * object:vector-interface[pos]. Retrun a pointer to * the corresponding object. */ static IBPtr getObjectFromNoun(string noun); //@} /** @name Access references between object in the repository. */ //@{ /** * Get referring objects. Return all object which refers to the * given object through a Reference of RefVector interface. */ static IVector GetObjectsReferringTo(IBPtr); /** * Get direct references. Return all objects the given object refers * to directly through a Reference of RefVector interface. */ static IVector DirectReferences(IBPtr); /** * Get all references. If \a obj contains references to other objects, * either through a Reference or RefVector interface or through the * virtual getReferences member function, add these to refs. Do the * same to the references recursively. */ static void addReferences(tIBPtr obj, ObjectSet & refs); //@} /** @name Access the interfaces of the objects in the repository. */ //@{ /** * Get interfaces. Return the interfaces defined for the * InterfacedBase class with the given type_info, \a ti, mapped to * their name. If several interfaces with the same name exists only * the one which correspond to the most derived class will be given, * except if \a all is true in which case all interfaces are given * (prefixed by '+'s to become unique). */ static InterfaceMap getInterfaces(const type_info & ti, bool all = true); /** * Return an interface with the given \a name to the given \a object. */ static const InterfaceBase * FindInterface(IBPtr object, string name); /** * Get an interface name. Decompose a string of the form * object:interface or * object:vector-interface[pos]. Return the interface * name (without the [pos]). */ static string getInterfaceFromNoun(string noun); /** * Get interface index. Decompose a string of the form * object:interface or * object:vector-interface[pos]. Return the * pos part or empty string if not present. */ static string getPosArgFromNoun(string noun); /** * Return a list of the interfaces which do not have their default * values for the given objects. */ template static vector< pair > getNonDefaultInterfaces(const Cont &); //@} /** @name Manipulate objects in the repository. */ //@{ /** * Call the InterfacedBase::update() function of all objects. */ static void update(); /** * Clear the InterfacedBase::touched() flag in all objects in the * given container. */ template static void clearAll(const Cont & c) { - for_each(c, mem_fun(&InterfacedBase::clear)); + for_each(c, mem_fn(&InterfacedBase::clear)); } /** * Set the status of all objects in the given container to * InterfacedBase::uninitialized. */ template static void resetAll(const Cont & c) { - for_each(c, mem_fun(&InterfacedBase::reset)); + for_each(c, mem_fn(&InterfacedBase::reset)); } /** * Setup an object. Execute the InterfacedBase::readSetup() method * of \a ip with the stream \a is as argument. */ static void readSetup(tIBPtr ip, istream & is); /** * Lock the given object. Locked objects cannot be * changed through an interface. */ static void lock(tIBPtr ip) { ip->lock(); } /** * Unlock the given object. Locked objects cannot be changed through * an interface. */ static void unlock(tIBPtr ip) { ip->unlock(); } //@} /** @name Access the documentation of objects. */ //@{ /** * Return the class documentation of a given object */ static const ClassDocumentationBase * getDocumentation(tcIBPtr ip); /** * Get the description for the model implemented in the class of the * given object. */ static string getModelDescription(tcIBPtr ip); /** * Get the references for the model implemented in the class of the * given object. */ static string getModelReferences(tcIBPtr ip); //@} /** @name Manipulate the output streams of the repository. */ //@{ /** * Set the standard output stream */ static void cout(ostream & os) { coutp() = &os; } /** * Get the standard output stream */ static ostream & cout() { return *coutp(); } /** * Set the standard error stream */ static void cerr(ostream & os) { cerrp() = &os; } /** * Get the standard error stream */ static ostream & cerr() { return *cerrp(); } /** * Set the standard log stream */ static void clog(ostream & os) { clogp() = &os; } /** * Get the standard log stream */ static ostream & clog() { return *clogp(); } //@} protected: /** @name Access standard InterfacedBase functions. */ //@{ /** * Return a clone of the given object. Calls the * InterfacedBase::clone() function of \a t and casts the resulting * pointer to the correct type. */ template static typename Ptr::pointer clone(const T & t); /** * Return a clone of the given object. Calls the * InterfacedBase::fullclone() function of \a t and casts the * resulting pointer to the correct type. */ template static typename Ptr::pointer fullclone(const T & t); /** * Rebind references. For all objects directly referenced by \a obj, * replace them with the translation found in \a trans. If \a obj has a * Reference or a member of a RefVector interface which is null, and * the corresponding interface has the RefInterfaceBase::defaultIfNull() flag set, * translate the null pointer to the first acceptable object in * defaults. */ static void rebind(InterfacedBase & obj, const TranslationMap & trans, const IVector & defaults); //@} /** * Add interfaces to the given map for the class with the given * class description. Recursively do the same with the base classes. */ static void addInterfaces(const ClassDescriptionBase &, InterfaceMap &, bool all = true); /** @name Functions containing the static instances of objects used by the repository. */ //@{ /** * All InterfacedBase objects mapped to their name. */ static ObjectMap & objects(); /** * All InterfacedBase objects. */ static ObjectSet & allObjects(); /** * Sets of InterfaceBase objects mapped to the class description of * the class for which they are defined. */ static TypeInterfaceMap & interfaces(); /** * Sets of ClassDocumentationBase objects mapped to the class * description of the class for which they are defined. */ static TypeDocumentationMap & documentations(); /** * All defined directories. */ static DirectorySet & directories(); /** * The current directory stack. */ static StringVector & directoryStack(); /** * Flag to say if we are in the middle of an update procedure. */ static bool & updating(); /** * The current current standard output stream. */ static ostream *& coutp(); /** * The current current standard error stream. */ static ostream *& cerrp(); /** * The current current standard log stream. */ static ostream *& clogp(); //@} }; } #ifndef ThePEG_TEMPLATES_IN_CC_FILE #include "BaseRepository.tcc" #endif #endif /* ThePEG_BaseRepository_H */ diff --git a/Repository/EventGenerator.cc b/Repository/EventGenerator.cc --- a/Repository/EventGenerator.cc +++ b/Repository/EventGenerator.cc @@ -1,1415 +1,1415 @@ // -*- C++ -*- // // EventGenerator.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 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 EventGenerator class. // #include "EventGenerator.h" #include "EventGenerator.xh" #include "ThePEG/Handlers/EventHandler.h" #include "Repository.h" #include "ThePEG/Utilities/HoldFlag.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/Utilities/DebugItem.h" #include "ThePEG/Interface/Interfaced.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/RefVector.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/PDT/ParticleData.h" #include "ThePEG/PDT/MatcherBase.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/StandardModel/StandardModelBase.h" #include "ThePEG/Repository/Strategy.h" #include "ThePEG/Repository/CurrentGenerator.h" #include "ThePEG/Handlers/AnalysisHandler.h" #include "ThePEG/Analysis/FactoryBase.h" #include "ThePEG/Handlers/EventManipulator.h" #include "ThePEG/Handlers/LuminosityFunction.h" #include "ThePEG/MatrixElement/MEBase.h" #include "ThePEG/EventRecord/Event.h" #include "ThePEG/Handlers/SubProcessHandler.h" #include "ThePEG/Handlers/CascadeHandler.h" #include "ThePEG/Handlers/HadronizationHandler.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Config/algorithm.h" #include "ThePEG/Utilities/DynamicLoader.h" #include #include "ThePEG/Repository/Main.h" #include #ifdef ThePEG_TEMPLATES_IN_CC_FILE #include "EventGenerator.tcc" #endif using namespace ThePEG; namespace { volatile sig_atomic_t THEPEG_SIGNAL_STATE = 0; } // signal handler function // very restricted in what it is allowed do // without causing undefined behaviour extern "C" { void thepegSignalHandler(int id) { THEPEG_SIGNAL_STATE=id; signal(id,SIG_DFL); } } void EventGenerator::checkSignalState() { if (THEPEG_SIGNAL_STATE) { log() << "Caught signal " << THEPEG_SIGNAL_STATE << ". Exiting ..." << std::endl; finalize(); exit(0); } } EventGenerator::EventGenerator() : thePath("."), theNumberOfEvents(1000), theQuickSize(7000), preinitializing(false), ieve(0), weightSum(0.0), theDebugLevel(0), logNonDefault(-1), printEvent(0), dumpPeriod(0), keepAllDumps(false), debugEvent(0), maxWarnings(10), maxErrors(10), theCurrentRandom(0), theCurrentGenerator(0), useStdout(false), theIntermediateOutput(false) {} EventGenerator::EventGenerator(const EventGenerator & eg) : Interfaced(eg), theDefaultObjects(eg.theDefaultObjects), theLocalParticles(eg.theLocalParticles), theStandardModel(eg.theStandardModel), theStrategy(eg.theStrategy), theRandom(eg.theRandom), theEventHandler(eg.theEventHandler), theAnalysisHandlers(eg.theAnalysisHandlers), theHistogramFactory(eg.theHistogramFactory), theEventManipulator(eg.theEventManipulator), thePath(eg.thePath), theRunName(eg.theRunName), theNumberOfEvents(eg.theNumberOfEvents), theObjects(eg.theObjects), theObjectMap(eg.theObjectMap), theParticles(eg.theParticles), theQuickParticles(eg.theQuickParticles), theQuickSize(eg.theQuickSize), preinitializing(false), theMatchers(eg.theMatchers), usedObjects(eg.usedObjects), ieve(eg.ieve), weightSum(eg.weightSum), theDebugLevel(eg.theDebugLevel), logNonDefault(eg.logNonDefault), printEvent(eg.printEvent), dumpPeriod(eg.dumpPeriod), keepAllDumps(eg.keepAllDumps), debugEvent(eg.debugEvent), maxWarnings(eg.maxWarnings), maxErrors(eg.maxErrors), theCurrentRandom(0), theCurrentGenerator(0), theCurrentEventHandler(eg.theCurrentEventHandler), theCurrentStepHandler(eg.theCurrentStepHandler), useStdout(eg.useStdout), theIntermediateOutput(eg.theIntermediateOutput) {} EventGenerator::~EventGenerator() { if ( theCurrentRandom ) delete theCurrentRandom; if ( theCurrentGenerator ) delete theCurrentGenerator; } IBPtr EventGenerator::clone() const { return new_ptr(*this); } IBPtr EventGenerator::fullclone() const { return new_ptr(*this); } tcEventPtr EventGenerator::currentEvent() const { return eventHandler()->currentEvent(); } CrossSection EventGenerator::histogramScale() const { return eventHandler()->histogramScale(); } CrossSection EventGenerator::integratedXSec() const { return eventHandler()->integratedXSec(); } CrossSection EventGenerator::integratedXSecErr() const { return eventHandler()->integratedXSecErr(); } void EventGenerator::setSeed(long seed) { random().setSeed(seed); ostringstream s; s << seed; const InterfaceBase * ifb = BaseRepository::FindInterface(theRandom, "Seed"); ifb->exec(*theRandom, "set", s.str()); } void EventGenerator::setup(string newRunName, ObjectSet & newObjects, ParticleMap & newParticles, MatcherSet & newMatchers) { HoldFlag debug(Debug::level, Debug::isset? Debug::level: theDebugLevel); theRunName = newRunName; theObjects.swap(newObjects); theParticles.swap(newParticles); theMatchers.swap(newMatchers); theObjectMap.clear(); for ( ObjectSet::const_iterator it = objects().begin(); it != objects().end(); ++it ) theObjectMap[(**it).fullName()] = *it; UseRandom currentRandom(theRandom); CurrentGenerator currentGenerator(this); // Force update of all objects and then reset. touch(); - for_each(theObjects, mem_fun(&InterfacedBase::touch)); + for_each(theObjects, mem_fn(&InterfacedBase::touch)); update(); - for_each(theObjects, mem_fun(&InterfacedBase::update)); + for_each(theObjects, mem_fn(&InterfacedBase::update)); clear(); BaseRepository::clearAll(theObjects); init(); } IBPtr EventGenerator::getPointer(string name) const { ObjectMap::const_iterator it = objectMap().find(name); if ( it == objectMap().end() ) return IBPtr(); else return it->second; } void EventGenerator::openOutputFiles() { if ( !useStdout ) { logfile().open((filename() + ".log").c_str()); theOutFileName = filename() + ".out"; outfile().open(theOutFileName.c_str()); outfile().close(); theOutStream.str(""); } out() << Repository::banner() << endl; log() << Repository::banner() << endl; } void EventGenerator::closeOutputFiles() { flushOutputFile(); if ( !useStdout ) logfile().close(); } void EventGenerator::flushOutputFile() { if ( !useStdout ) { outfile().open(theOutFileName.c_str(), ios::out|ios::app); outfile() << theOutStream.str(); outfile().close(); } else BaseRepository::cout() << theOutStream.str(); theOutStream.str(""); } void EventGenerator::doinit() { HoldFlag debug(Debug::level, Debug::isset? Debug::level: theDebugLevel); // First initialize base class and random number generator. Interfaced::doinit(); random().init(); // Make random generator and this available in standard static // classes. UseRandom useRandom(theRandom); CurrentGenerator currentGenerator(this); // First initialize all objects which have requested this by // implementing a InterfacedBase::preInitialize() function which // returns true. while ( true ) { HoldFlag hold(preinitializing, true); ObjectSet preinits; for ( ObjectSet::iterator it = objects().begin(); it != objects().end(); ++it ) if ( (**it).preInitialize() && (**it).state() == InterfacedBase::uninitialized ) preinits.insert(*it); if ( preinits.empty() ) break; - for_each(preinits, mem_fun(&InterfacedBase::init)); + for_each(preinits, std::mem_fn(&InterfacedBase::init)); } // Initialize the quick access to particles. theQuickParticles.clear(); theQuickParticles.resize(2*theQuickSize); for ( ParticleMap::const_iterator pit = theParticles.begin(); pit != theParticles.end(); ++pit ) if ( abs(pit->second->id()) < theQuickSize ) theQuickParticles[pit->second->id()+theQuickSize] = pit->second; // Then call the init method for all objects. Start with the // standard model and the strategy. standardModel()->init(); if ( strategy() ) strategy()->init(); eventHandler()->init(); // initialize particles first for(ParticleMap::const_iterator pit = particles().begin(); pit != particles().end(); ++pit) pit->second->init(); - for_each(objects(), mem_fun(&InterfacedBase::init)); + for_each(objects(), std::mem_fn(&InterfacedBase::init)); // Then initialize the Event Handler calculating initial cross // sections and stuff. eventHandler()->initialize(); } void EventGenerator::doinitrun() { HoldFlag debug(Debug::level, Debug::isset? Debug::level: theDebugLevel); signal(SIGHUP, thepegSignalHandler); signal(SIGINT, thepegSignalHandler); signal(SIGTERM,thepegSignalHandler); currentEventHandler(eventHandler()); Interfaced::doinitrun(); random().initrun(); // Then call the init method for all objects. Start with the // standard model and the strategy. standardModel()->initrun(); if ( strategy() ) { strategy()->initrun(); if ( ! strategy()->versionstring().empty() ) { out() << ">> " << strategy()->versionstring() << '\n' << endl; log() << ">> " << strategy()->versionstring() << '\n' << endl; } } // initialize particles first for(ParticleMap::const_iterator pit = particles().begin(); pit != particles().end(); ++pit) { pit->second->initrun(); } eventHandler()->initrun(); - for_each(objects(), mem_fun(&InterfacedBase::initrun)); + for_each(objects(), std::mem_fn(&InterfacedBase::initrun)); if ( logNonDefault > 0 || ( ThePEG_DEBUG_LEVEL && logNonDefault == 0 ) ) { vector< pair > changed = Repository::getNonDefaultInterfaces(objects()); if ( changed.size() ) { log() << string(78, '=') << endl << "The following interfaces have non-default values (default):" << endl << string(78, '-') << endl; for ( int i = 0, N = changed.size(); i < N; ++i ) { log() << changed[i].first->fullName() << ":" << changed[i].second->name() << " = " << changed[i].second->exec(*changed[i].first, "notdef", "") << endl; } log() << string(78,'=') << endl; } } weightSum = 0.0; } PDPtr EventGenerator::getParticleData(PID id) const { long newId = id; if ( abs(newId) < theQuickSize && theQuickParticles.size() ) return theQuickParticles[newId+theQuickSize]; ParticleMap::const_iterator it = theParticles.find(newId); if ( it == theParticles.end() ) return PDPtr(); return it->second; } PPtr EventGenerator::getParticle(PID newId) const { tcPDPtr pd = getParticleData(newId); if ( !pd ) return PPtr(); return pd->produceParticle(); } void EventGenerator::finalize() { UseRandom currentRandom(theRandom); CurrentGenerator currentGenerator(this); finish(); finally(); } void EventGenerator::dofinish() { HoldFlag debug(Debug::level, Debug::isset? Debug::level: theDebugLevel); // first write out statistics from the event handler. eventHandler()->statistics(out()); // Call the finish method for all other objects. - for_each(objects(), mem_fun(&InterfacedBase::finish)); + for_each(objects(), std::mem_fn(&InterfacedBase::finish)); if ( theExceptions.empty() ) { log() << "No exceptions reported in this run.\n"; } else { log() << "\nThe following exception classes were reported in this run:\n"; for ( ExceptionMap::iterator it = theExceptions.begin(); it != theExceptions.end(); ++it ) { string severity; switch ( it->first.second ) { case Exception::info : severity="info"; break; case Exception::warning : severity="warning"; break; case Exception::setuperror : severity="setuperror"; break; case Exception::eventerror : severity="eventerror"; break; case Exception::runerror : severity="runerror"; break; case Exception::maybeabort : severity="maybeabort"; break; case Exception::abortnow : severity="abortnow"; break; default : severity="unknown"; } log() << it->first.first << ' ' << severity << " (" << it->second << " times)\n"; } } theExceptions.clear(); const string & msg = theMiscStream.str(); if ( ! msg.empty() ) { log() << endl << "Miscellaneous output from modules to the standard output:\n\n" << msg; theMiscStream.str(""); } flushOutputFile(); } void EventGenerator::finally() { generateReferences(); closeOutputFiles(); if ( theCurrentRandom ) delete theCurrentRandom; if ( theCurrentGenerator ) delete theCurrentGenerator; theCurrentRandom = 0; theCurrentGenerator = 0; } void EventGenerator::initialize(bool initOnly) { UseRandom currentRandom(theRandom); CurrentGenerator currentGenerator(this); doInitialize(initOnly); } bool EventGenerator::loadMain(string file) { initialize(); UseRandom currentRandom(theRandom); CurrentGenerator currentGenerator(this); Main::eventGenerator(this); bool ok = DynamicLoader::load(file); finish(); finally(); return ok; } void EventGenerator::go(long next, long maxevent, bool tics) { UseRandom currentRandom(theRandom); CurrentGenerator currentGenerator(this); doGo(next, maxevent, tics); } EventPtr EventGenerator::shoot() { static DebugItem debugfpu("ThePEG::FPU", 1); if ( debugfpu ) Debug::unmaskFpuErrors(); UseRandom currentRandom(theRandom); CurrentGenerator currentGenerator(this); checkSignalState(); EventPtr event = doShoot(); if ( event ) weightSum += event->weight(); DebugItem::tic(); return event; } EventPtr EventGenerator::doShoot() { EventPtr event; if ( N() >= 0 && ++ieve > N() ) return event; HoldFlag debug(Debug::level, Debug::isset? Debug::level: theDebugLevel); do { int state = 0; int loop = 1; eventHandler()->clearEvent(); try { do { // Generate a full event or part of an event if ( eventHandler()->empty() ) event = eventHandler()->generateEvent(); else event = eventHandler()->continueEvent(); if ( eventHandler()->empty() ) loop = -loop; // Analyze the possibly uncomplete event for ( AnalysisVector::iterator it = analysisHandlers().begin(); it != analysisHandlers().end(); ++it ) (**it).analyze(event, ieve, loop, state); // Manipulate the current event, possibly deleting some steps // and telling the event handler to redo them. if ( manipulator() ) state = manipulator()->manipulate(eventHandler(), event); // If the event was not completed, continue generation and continue. loop = abs(loop) + 1; } while ( !eventHandler()->empty() ); } catch (Exception & ex) { if ( logException(ex, eventHandler()->currentEvent()) ) throw; } catch (...) { event = eventHandler()->currentEvent(); if ( event ) log() << *event; else log() << "An exception occurred before any event object was created!"; log() << endl; if ( ThePEG_DEBUG_LEVEL ) dump(); throw; } if ( ThePEG_DEBUG_LEVEL ) { if ( ( ThePEG_DEBUG_LEVEL == Debug::printEveryEvent || ieve < printEvent ) && event ) log() << *event; if ( debugEvent > 0 && ieve + 1 >= debugEvent ) Debug::level = Debug::full; } } while ( !event ); // If scheduled, dump a clean state between events if ( ThePEG_DEBUG_LEVEL && dumpPeriod > 0 && ieve%dumpPeriod == 0 ) { eventHandler()->clearEvent(); eventHandler()->clean(); dump(); } return event; } EventPtr EventGenerator::doGenerateEvent(tEventPtr e) { if ( N() >= 0 && ++ieve > N() ) return EventPtr(); EventPtr event = e; try { event = eventHandler()->generateEvent(e); } catch (Exception & ex) { if ( logException(ex, eventHandler()->currentEvent()) ) throw; } catch (...) { event = eventHandler()->currentEvent(); if ( !event ) event = e; log() << *event << endl; dump(); throw; } return event; } EventPtr EventGenerator::doGenerateEvent(tStepPtr s) { if ( N() >= 0 && ++ieve > N() ) return EventPtr(); EventPtr event; try { event = eventHandler()->generateEvent(s); } catch (Exception & ex) { if ( logException(ex, eventHandler()->currentEvent()) ) throw; } catch (...) { event = eventHandler()->currentEvent(); if ( event ) log() << *event << endl; dump(); throw; } return event; } EventPtr EventGenerator::generateEvent(Event & e) { UseRandom currentRandom(theRandom); CurrentGenerator currentGenerator(this); EventPtr event = doGenerateEvent(tEventPtr(&e)); if ( event ) weightSum += event->weight(); return event; } EventPtr EventGenerator::generateEvent(Step & s) { UseRandom currentRandom(theRandom); CurrentGenerator currentGenerator(this); EventPtr event = doGenerateEvent(tStepPtr(&s)); if ( event ) weightSum += event->weight(); return event; } Energy EventGenerator::maximumCMEnergy() const { tcEHPtr eh = eventHandler(); return eh->lumiFnPtr()? eh->lumiFn().maximumCMEnergy(): ZERO; } void EventGenerator::doInitialize(bool initOnly) { if ( !initOnly ) openOutputFiles(); init(); if ( !initOnly ) initrun(); if ( !ThePEG_DEBUG_LEVEL ) Exception::noabort = true; } void EventGenerator::doGo(long next, long maxevent, bool tics) { if ( maxevent >= 0 ) N(maxevent); if ( next >= 0 ) { if ( tics ) cerr << "event> " << setw(9) << "init\r" << flush; initialize(); ieve = next-1; } else { openOutputFiles(); } if ( tics ) tic(); try { while ( shoot() ) { if ( tics ) tic(); } } catch ( ... ) { finish(); throw; } finish(); finally(); } void EventGenerator::tic(long currev, long totev) const { if ( !currev ) currev = ieve; if ( !totev ) totev = N(); long i = currev; long n = totev; bool skip = currev%(max(totev/100, 1L)); if ( i > n/2 ) i = n-i; while ( skip && i >= 10 && !(i%10) ) i /= 10; if ( i == 1 || i == 2 || i == 5 ) skip = false; if (!theIntermediateOutput) { //default if ( skip ) return; cerr << "event> " << setw(8) << currev << " " << setw(8) << totev << "\r"; } else if (theIntermediateOutput) { if ( skip && currev%10000!=0) return; cerr << "event> " << setw(9) << right << currev << "/" << totev << "; xs = " << integratedXSec()/picobarn << " pb +- " << integratedXSecErr()/picobarn << " pb" << endl; } cerr.flush(); if ( currev == totev ) cerr << endl; } void EventGenerator::dump() const { if ( dumpPeriod > -1 ) { string dumpfile; if ( keepAllDumps ) { ostringstream number; number << ieve; dumpfile = filename() + "-" + number.str() + ".dump"; } else dumpfile = filename() + ".dump"; PersistentOStream file(dumpfile, globalLibraries()); file << tcEGPtr(this); } } void EventGenerator::use(const Interfaced & i) { IBPtr ip = getPtr(i); if ( ip ) usedObjects.insert(ip); } void EventGenerator::generateReferences() { typedef map StringMap; StringMap references; // First get all model descriptions and model references from the // used objects. Put them in a map indexed by the description to // avoid duplicates. for ( ObjectSet::iterator it = usedObjects.begin(); it != usedObjects.end(); ++it ) { if ( *it == strategy() ) continue; string desc = Repository::getModelDescription(*it); if ( desc.empty() ) continue; if ( dynamic_ptr_cast(*it) ) desc = "A " + desc; else if ( dynamic_ptr_cast(*it) ) desc = "B " + desc; else if ( dynamic_ptr_cast(*it) ) desc = "C " + desc; else if ( dynamic_ptr_cast(*it) ) desc = "D " + desc; else if ( dynamic_ptr_cast(*it) ) desc = "E " + desc; else if ( dynamic_ptr_cast(*it) ) desc = "F " + desc; else if ( dynamic_ptr_cast(*it) ) desc = "Y " + desc; else if ( dynamic_ptr_cast(*it) ) desc = "Z " + desc; else if ( dynamic_ptr_cast::const_pointer>(*it) ) desc = "G " + desc; else desc = "H " + desc; references[desc] = Repository::getModelReferences(*it); } // Now get the main strategy description which should put first and // remove it from the map. string stratdesc; string stratref; if ( strategy() ) { stratdesc = Repository::getModelDescription(strategy()); stratref = Repository::getModelReferences(strategy()); references.erase(stratdesc); } // Open the file and write out an appendix header if ( !useStdout ) reffile().open((filename() + ".tex").c_str()); ref() << "\\documentclass{article}\n" << "\\usepackage{graphics}\n" << "\\begin{document}\n" << "\\appendix\n" << "\\section[xxx]{\\textsc{ThePEG} version " << Repository::version() << " \\cite{ThePEG} Run Information}\n" << "Run name: \\textbf{" << runName() << "}:\\\\\n"; if ( !stratdesc.empty() ) ref() << "This run was generated using " << stratdesc << " and the following models:\n"; else ref() << "The following models were used:\n"; ref() << "\\begin{itemize}\n"; // Write out all descriptions. for ( StringMap::iterator it = references.begin(); it != references.end(); ++it ) ref() << "\\item " << it->first.substr(2) << endl; // Write out thebibliography header and all references. ref() << "\\end{itemize}\n\n" << "\\begin{thebibliography}{99}\n" << "\\bibitem{ThePEG} L.~L\\\"onnblad, " << "Comput.~Phys.~Commun.\\ {\\bf 118} (1999) 213.\n"; if ( !stratref.empty() ) ref() << stratref << '\n'; for ( StringMap::iterator it = references.begin(); it != references.end(); ++it ) ref() << it->second << '\n'; ref() << "\\end{thebibliography}\n" << "\\end{document}" << endl; if ( !useStdout ) reffile().close(); } void EventGenerator::strategy(StrategyPtr s) { theStrategy = s; } int EventGenerator::count(const Exception & ex) { return ++theExceptions[make_pair(StringUtils::typeName(typeid(ex)), ex.severity())]; } void EventGenerator::printException(const Exception & ex) { switch ( ex.severity() ) { case Exception::info: log() << "* An information"; break; case Exception::warning: log() << "* A warning"; break; case Exception::setuperror: log() << "** A setup"; break; case Exception::eventerror: log() << "** An event"; break; case Exception::runerror: log() << "*** An run"; break; case Exception::maybeabort: case Exception::abortnow: log() << "**** A serious"; break; default: log() << "**** An unknown"; break; } if ( ieve > 0 ) log() << " exception of type " << StringUtils::typeName(typeid(ex)) << " occurred while generating event number " << ieve << ": \n" << ex.message() << endl; else log() << " exception occurred in the initialization of " << name() << ": \n" << ex.message() << endl; if ( ex.severity() == Exception::eventerror ) log() << "The event will be discarded." << endl; } void EventGenerator::logWarning(const Exception & ex) { if ( ex.severity() != Exception::info && ex.severity() != Exception::warning ) throw ex; ex.handle(); int c = count(ex); if ( c > maxWarnings ) return; printException(ex); if ( c == maxWarnings ) log() << "No more warnings of this kind will be reported." << endl; } bool EventGenerator:: logException(const Exception & ex, tcEventPtr event) { bool noEvent = !event; ex.handle(); int c = count(ex); if ( c <= maxWarnings ) { printException(ex); if ( c == maxWarnings ) log() << "No more warnings of this kind will be reported." << endl; } if ( ex.severity() == Exception::info || ex.severity() == Exception::warning ) { ex.handle(); return false; } if ( ex.severity() == Exception::eventerror ) { if ( c < maxErrors || maxErrors <= 0 ) { ex.handle(); if ( ThePEG_DEBUG_LEVEL > 0 && !noEvent ) log() << *event; return false; } if ( c > maxErrors ) printException(ex); log() << "Too many (" << c << ") exceptions of this kind has occurred. " "Execution will be stopped.\n"; } else { log() << "This exception is too serious. Execution will be stopped.\n"; } if ( !noEvent ) log() << *event; else log() << "An exception occurred before any event object was created!\n"; dump(); return true; } struct MatcherOrdering { bool operator()(tcPMPtr m1, tcPMPtr m2) const { return m1->name() < m2->name() || ( m1->name() == m2->name() && m1->fullName() < m2->fullName() ); } }; struct ObjectOrdering { bool operator()(tcIBPtr i1, tcIBPtr i2) const { return i1->fullName() < i2->fullName(); } }; void EventGenerator::persistentOutput(PersistentOStream & os) const { set match(theMatchers.begin(), theMatchers.end()); set usedset(usedObjects.begin(), usedObjects.end()); os << theDefaultObjects << theLocalParticles << theStandardModel << theStrategy << theRandom << theEventHandler << theAnalysisHandlers << theHistogramFactory << theEventManipulator << thePath << theRunName << theNumberOfEvents << theObjectMap << theParticles << theQuickParticles << theQuickSize << match << usedset << ieve << weightSum << theDebugLevel << logNonDefault << printEvent << dumpPeriod << keepAllDumps << debugEvent << maxWarnings << maxErrors << theCurrentEventHandler << theCurrentStepHandler << useStdout << theIntermediateOutput << theMiscStream.str() << Repository::listReadDirs(); } void EventGenerator::persistentInput(PersistentIStream & is, int) { string dummy; vector readdirs; theGlobalLibraries = is.globalLibraries(); is >> theDefaultObjects >> theLocalParticles >> theStandardModel >> theStrategy >> theRandom >> theEventHandler >> theAnalysisHandlers >> theHistogramFactory >> theEventManipulator >> thePath >> theRunName >> theNumberOfEvents >> theObjectMap >> theParticles >> theQuickParticles >> theQuickSize >> theMatchers >> usedObjects >> ieve >> weightSum >> theDebugLevel >> logNonDefault >> printEvent >> dumpPeriod >> keepAllDumps >> debugEvent >> maxWarnings >> maxErrors >> theCurrentEventHandler >> theCurrentStepHandler >> useStdout >> theIntermediateOutput >> dummy >> readdirs; theMiscStream.str(dummy); theMiscStream.seekp(0, std::ios::end); theObjects.clear(); for ( ObjectMap::iterator it = theObjectMap.begin(); it != theObjectMap.end(); ++it ) theObjects.insert(it->second); Repository::appendReadDir(readdirs); } void EventGenerator::setLocalParticles(PDPtr pd, int) { localParticles()[pd->id()] = pd; } void EventGenerator::insLocalParticles(PDPtr pd, int) { localParticles()[pd->id()] = pd; } void EventGenerator::delLocalParticles(int place) { ParticleMap::iterator it = localParticles().begin(); while ( place-- && it != localParticles().end() ) ++it; if ( it != localParticles().end() ) localParticles().erase(it); } vector EventGenerator::getLocalParticles() const { vector ret; for ( ParticleMap::const_iterator it = localParticles().begin(); it != localParticles().end(); ++it ) ret.push_back(it->second); return ret; } void EventGenerator::setPath(string newPath) { if ( std::system(("mkdir -p " + newPath).c_str()) ) throw EGNoPath(newPath); if ( std::system(("touch " + newPath + "/.ThePEG").c_str()) ) throw EGNoPath(newPath); if ( std::system(("rm -f " + newPath + "/.ThePEG").c_str()) ) throw EGNoPath(newPath); thePath = newPath; } string EventGenerator::defPath() const { char * env = std::getenv("ThePEG_RUN_DIR"); if ( env ) return string(env); return string("."); } ostream & EventGenerator::out() { return theOutStream; } ostream & EventGenerator::log() { return logfile().is_open()? logfile(): BaseRepository::cout(); } ostream & EventGenerator::ref() { return reffile().is_open()? reffile(): BaseRepository::cout(); } string EventGenerator::doSaveRun(string runname) { runname = StringUtils::car(runname); if ( runname.empty() ) runname = theRunName; if ( runname.empty() ) runname = name(); EGPtr eg = Repository::makeRun(this, runname); string file = eg->filename() + ".run"; PersistentOStream os(file); os << eg; if ( !os ) return "Error: Save failed! (I/O error)"; return ""; } string EventGenerator::doMakeRun(string runname) { runname = StringUtils::car(runname); if ( runname.empty() ) runname = theRunName; if ( runname.empty() ) runname = name(); Repository::makeRun(this, runname); return ""; } bool EventGenerator::preinitRegister(IPtr obj, string fullname) { if ( !preinitializing ) throw InitException() << "Tried to register a new object in the initialization of an " << "EventGenerator outside of the pre-initialization face. " << "The preinitRegister() can only be called from a doinit() function " << "in an object for which preInitialize() returns true."; if ( objectMap().find(fullname) != objectMap().end() ) return false; obj->name(fullname); objectMap()[fullname] = obj; objects().insert(obj); obj->theGenerator = this; PDPtr pd = dynamic_ptr_cast(obj); if ( pd ) theParticles[pd->id()] = pd; PMPtr pm = dynamic_ptr_cast(obj); if ( pm ) theMatchers.insert(pm); return true; } bool EventGenerator::preinitRemove(IPtr obj) { bool deleted=true; if(theObjects.find(obj)!=theObjects.end()) { theObjects.erase(obj); } else deleted = false; if(theObjectMap.find(obj->fullName())!=theObjectMap.end()) theObjectMap.erase(obj->fullName()); else deleted = false; return deleted; } IPtr EventGenerator:: preinitCreate(string classname, string fullname, string libraries) { if ( !preinitializing ) throw InitException() << "Tried to create a new object in the initialization of an " << "EventGenerator outside of the pre-initialization face. " << "The preinitCreate() can only be called from a doinit() function " << "in an object for which preInitialize() returns true."; if ( objectMap().find(fullname) != objectMap().end() ) return IPtr(); const ClassDescriptionBase * db = DescriptionList::find(classname); while ( !db && libraries.length() ) { string library = StringUtils::car(libraries); libraries = StringUtils::cdr(libraries); DynamicLoader::load(library); db = DescriptionList::find(classname); } if ( !db ) return IPtr(); IPtr obj = dynamic_ptr_cast(db->create()); if ( !obj ) return IPtr(); if ( !preinitRegister(obj, fullname) ) return IPtr(); return obj; } string EventGenerator:: preinitInterface(IPtr obj, string ifcname, string cmd, string value) { if ( !preinitializing ) throw InitException() << "Tried to manipulate an external object in the initialization of an " << "EventGenerator outside of the pre-initialization face. " << "The preinitSet() can only be called from a doinit() function " << "in an object for which preInitialize() returns true."; if ( !obj ) return "Error: No object found."; const InterfaceBase * ifc = Repository::FindInterface(obj, ifcname); if ( !ifc ) return "Error: No such interface found."; try { return ifc->exec(*obj, cmd, value); } catch ( const InterfaceException & ex) { ex.handle(); return "Error: " + ex.message(); } } string EventGenerator:: preinitInterface(IPtr obj, string ifcname, int index, string cmd, string value) { ostringstream os; os << index; return preinitInterface(obj, ifcname, cmd, os.str() + " " + value); } string EventGenerator:: preinitInterface(string fullname, string ifcname, string cmd, string value) { return preinitInterface(getObject(fullname), ifcname, cmd, value); } string EventGenerator:: preinitInterface(string fullname, string ifcname, int index, string cmd, string value) { return preinitInterface(getObject(fullname), ifcname, index, cmd, value); } tDMPtr EventGenerator::findDecayMode(string tag) const { for ( ObjectSet::const_iterator it = objects().begin(); it != objects().end(); ++it ) { tDMPtr dm = dynamic_ptr_cast(*it); if ( dm && dm->tag() == tag ) return dm; } return tDMPtr(); } tDMPtr EventGenerator::preinitCreateDecayMode(string tag) { return constructDecayMode(tag); } DMPtr EventGenerator::constructDecayMode(string & tag) { 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; } } rdm = findDecayMode(tag.substr(0,end)); if ( rdm ) return rdm; string::size_type next = tag.find("->"); if ( next == string::npos ) return rdm; if ( tag.find(';') == string::npos ) return rdm; tPDPtr pd = getObject(tag.substr(0,next)); if ( !pd ) pd = findParticle(tag.substr(0,next)); if ( !pd ) return rdm; rdm = ptr_new(); rdm->parent(pd); if ( pd->CC() ) { adm = ptr_new(); 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); tDMPtr cdm = constructDecayMode(tag); if ( cdm ) rdm->addCascadeProduct(cdm); else error = true; } break; case '=': dolink = true; [[fallthrough]]; case ',': case ']': tag = tag.substr(1); break; case '?': { next = min(tag.find(','), tag.find(';')); tPMPtr pm = 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 = 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 = 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 = 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); DMPtr ndm = findDecayMode(rdm->tag()); if ( ndm ) return ndm; pd->addDecayMode(rdm); if ( !preinitRegister(rdm, pd->fullName() + "/" + rdm->tag()) ) return DMPtr(); if ( adm ) { preinitRegister(adm, pd->CC()->fullName() + "/" + adm->tag()); rdm->CC(adm); adm->CC(rdm); } return rdm; } tPDPtr EventGenerator::findParticle(string pdgname) const { for ( ParticleMap::const_iterator it = particles().begin(); it != particles().end(); ++it ) if ( it->second->PDGName() == pdgname ) return it->second; return tPDPtr(); } tPMPtr EventGenerator::findMatcher(string name) const { for ( MatcherSet::const_iterator it = matchers().begin(); it != matchers().end(); ++it ) if ( (**it).name() == name ) return *it; return tPMPtr(); } ClassDescription EventGenerator::initEventGenerator; void EventGenerator::Init() { static ClassDocumentation documentation ("This is the main class used to administer an event generation run. " "The actual generation of each event is handled by the assigned " "EventHandler object. When the event generator" "is properly set up it can be initialized with the command " "MakeRun and/or saved to a file with the command " "SaveRun. If saved to a file, the event generator " "can be read into another program to produce events. The file can also " "be read into the runThePEG program where a number of events " "determined by the parameter NumberOfEvents is " "generated with each event analysed by the list of assigned " "AnalysisHandlers."); static Reference interfaceStandardModel ("StandardModelParameters", "The ThePEG::StandardModelBase object to be used to access standard " "model parameters in this run.", &EventGenerator::theStandardModel, false, false, true, false); static Reference interfaceEventHandler ("EventHandler", "The ThePEG::EventHandler object to be used to generate the " "individual events in this run.", &EventGenerator::theEventHandler, false, false, true, false); static RefVector interfaceAnalysisHandlers ("AnalysisHandlers", "ThePEG::AnalysisHandler objects to be used to analyze the produced " "events in this run.", &EventGenerator::theAnalysisHandlers, 0, true, false, true, false); static Reference interfaceHistogramFactory ("HistogramFactory", "An associated factory object for handling histograms to be used by " "AnalysisHandlers.", &EventGenerator::theHistogramFactory, true, false, true, true, true); static Reference interfaceEventManip ("EventManipulator", "An ThePEG::EventManipulator called each time the generation of an " "event is stopped. The ThePEG::EventManipulator object is able to " "manipulate the generated event, as opposed to an " "ThePEG::AnalysisHandler which may only look at the event.", &EventGenerator::theEventManipulator, true, false, true, true); static RefVector interfaceLocalParticles ("LocalParticles", "Special versions of ThePEG::ParticleData objects to be used " "in this run. Note that to delete an object, its number in the list " "should be given, rather than its id number.", 0, 0, false, false, true, false, &EventGenerator::setLocalParticles, &EventGenerator::insLocalParticles, &EventGenerator::delLocalParticles, &EventGenerator::getLocalParticles); static RefVector interfaceDefaultObjects ("DefaultObjects", "A vector of pointers to default objects. In a ThePEG::Reference or " "ThePEG::RefVector interface with the defaultIfNull() flag set, if a " "null pointer is encountered this vector is gone through until an " "acceptable object is found in which case the null pointer is replaced " "by a pointer to this object.", &EventGenerator::theDefaultObjects, 0, true, false, true, false, false); static Reference interfaceStrategy ("Strategy", "An ThePEG::Strategy with additional ThePEG::ParticleData objects to " "be used in this run.", &EventGenerator::theStrategy, false, false, true, true); static Reference interfaceRandomGenerator ("RandomNumberGenerator", "An ThePEG::RandomGenerator object which should typically interaface to " "a CLHEP Random object. This will be the default random number generator " "for the run, but individual objects may use their own random generator " "if they wish.", &EventGenerator::theRandom, true, false, true, false); static Parameter interfacePath ("Path", "The directory where the output files are put.", &EventGenerator::thePath, ".", true, false, &EventGenerator::setPath, 0, &EventGenerator::defPath); interfacePath.directoryType(); static Parameter interfaceRunName ("RunName", "The name of this run. This name will be used in the output filenames. " "The files wil be placed in the directory specified by the " "Path parameter" "If empty the name of the event generator will be used instead.", &EventGenerator::theRunName, "", true, false, 0, 0, &EventGenerator::name); static Parameter interfaceNumberOfEvents ("NumberOfEvents", "The number of events to be generated in this run. If less than zero, " "the number of events is unlimited", &EventGenerator::theNumberOfEvents, 1000, -1, Constants::MaxInt, true, false, Interface::lowerlim); static Parameter interfaceDebugLevel ("DebugLevel", "The level of debug information sent out to the log file in the run. " "Level 0 only gives a limited ammount of warnings and error messages. " "Level 1 will print the first few events. " "Level 5 will print every event. " "Level 9 will print every step in every event.", &EventGenerator::theDebugLevel, 0, 0, 9, true, false, true); static Parameter interfacePrintEvent ("PrintEvent", "If the debug level is above zero, print the first 'PrintEvent' events.", &EventGenerator::printEvent, 0, 0, 1000, true, false, Interface::lowerlim); static Parameter interfaceDumpPeriod ("DumpPeriod", "If the debug level is above zero, dump the full state of the run every " "'DumpPeriod' events. Set it to -1 to disable dumping even in the case of errors.", &EventGenerator::dumpPeriod, 0, -1, Constants::MaxInt, true, false, Interface::lowerlim); static Switch interfaceKeepAllDumps ("KeepAllDumps", "Whether all dump files should be kept, labelled by event number.", &EventGenerator::keepAllDumps, false, true, false); static SwitchOption interfaceKeepAllDumpsYes (interfaceKeepAllDumps, "Yes", "Keep all dump files, labelled by event number.", true); static SwitchOption interfaceKeepAllDumpsNo (interfaceKeepAllDumps, "No", "Keep only the latest dump file.", false); static Parameter interfaceDebugEvent ("DebugEvent", "If the debug level is above zero, step up to the highest debug level " "befor event number 'DebugEvent'.", &EventGenerator::debugEvent, 0, 0, Constants::MaxInt, true, false, Interface::lowerlim); static Parameter interfaceMaxWarnings ("MaxWarnings", "The maximum number of warnings of each type which will be printed.", &EventGenerator::maxWarnings, 10, 1, 100, true, false, Interface::lowerlim); static Parameter interfaceMaxErrors ("MaxErrors", "The maximum number of errors of each type which will be tolerated. " "If more errors are reported, the run will be aborted.", &EventGenerator::maxErrors, 10, -1, 100000, true, false, Interface::lowerlim); static Parameter interfaceQuickSize ("QuickSize", "The max absolute id number of particle data objects which are accessed " "quickly through a vector indexed by the id number.", &EventGenerator::theQuickSize, 7000, 0, 50000, true, false, Interface::lowerlim); static Command interfaceSaveRun ("SaveRun", "Isolate, initialize and save this event generator to a file, from which " "it can be read in and run in another program. If an agument is given " "this is used as the run name, otherwise the run name is taken from the " "RunName parameter.", &EventGenerator::doSaveRun, true); static Command interfaceMakeRun ("MakeRun", "Isolate and initialize this event generator and give it a run name. " "If no argument is given, the run name is taken from the " "RunName parameter.", &EventGenerator::doMakeRun, true); interfaceEventHandler.rank(11.0); interfaceSaveRun.rank(10.0); interfaceMakeRun.rank(9.0); interfaceRunName.rank(8.0); interfaceNumberOfEvents.rank(7.0); interfaceAnalysisHandlers.rank(6.0); static Switch interfaceUseStdout ("UseStdout", "Redirect the logging and output to stdout instead of files.", &EventGenerator::useStdout, false, true, false); static SwitchOption interfaceUseStdoutYes (interfaceUseStdout, "Yes", "Use stdout instead of log files.", true); static SwitchOption interfaceUseStdoutNo (interfaceUseStdout, "No", "Use log files.", false); static Switch interfaceLogNonDefault ("LogNonDefault", "Controls the printout of important interfaces which has been changed from their default values.", &EventGenerator::logNonDefault, -1, true, false); static SwitchOption interfaceLogNonDefaultYes (interfaceLogNonDefault, "Yes", "Always print changed interfaces.", 1); static SwitchOption interfaceLogNonDefaultOnDebug (interfaceLogNonDefault, "OnDebug", "Only print changed interfaces if debugging is turned on.", 0); static SwitchOption interfaceLogNonDefaultNo (interfaceLogNonDefault, "No", "Don't print changed interfaces.", -1); interfaceLogNonDefault.setHasDefault(false); static Switch interfaceIntermediateOutput ("IntermediateOutput", "Modified event number count with the number of events processed so far, " "which updates at least every 10000 events, together with the corresponding " "intermediate estimate for the cross section plus the integration error.", &EventGenerator::theIntermediateOutput, false, true, false); static SwitchOption interfaceIntermediateOutputYes (interfaceIntermediateOutput, "Yes", "Show the modified event number count with the number of events processed so far, " "plus further information on the intermediate cross section estimate.", true); static SwitchOption interfaceIntermediateOutputNo (interfaceIntermediateOutput, "No", "Show the usual event number count with the number of events processed so far, " "but no further information on the intermediate cross section estimate.", false); } EGNoPath::EGNoPath(string path) { theMessage << "Cannot set the directory path for output files to '" << path << "' because the directory did not exist and could not be " << "created."; severity(warning); } diff --git a/Repository/MultiEventGenerator.cc b/Repository/MultiEventGenerator.cc --- a/Repository/MultiEventGenerator.cc +++ b/Repository/MultiEventGenerator.cc @@ -1,347 +1,347 @@ // -*- C++ -*- // // MultiEventGenerator.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 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 MultiEventGenerator class. // #include "MultiEventGenerator.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Command.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Repository/BaseRepository.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/Rebinder.h" #include "ThePEG/EventRecord/Event.h" #include "ThePEG/Config/algorithm.h" #include "ThePEG/Utilities/StringUtils.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/Repository/CurrentGenerator.h" #include using namespace ThePEG; MultiEventGenerator::~MultiEventGenerator() {} IBPtr MultiEventGenerator::clone() const { return new_ptr(*this); } IBPtr MultiEventGenerator::fullclone() const { return new_ptr(*this); } string MultiEventGenerator::removeInterface(string cmd) { string noun = StringUtils::car(cmd); IBPtr ip = BaseRepository::getObjectFromNoun(noun); const InterfaceBase * ifb = BaseRepository:: FindInterface(ip, BaseRepository::getInterfaceFromNoun(noun)); string posarg = BaseRepository::getPosArgFromNoun(noun); for ( string::size_type i = 0; i < theObjects.size(); ++i ) { if ( theObjects[i] == ip && theInterfaces[i] == ifb->name() && thePosArgs[i] == posarg ) { theObjects.erase(theObjects.begin() + i); theInterfaces.erase(theInterfaces.begin() + i); thePosArgs.erase(thePosArgs.begin() + i); theValues.erase(theValues.begin() + i); return ""; } } return "No such object/interface defined for this MultiEventGenerator."; } string MultiEventGenerator::addInterface(string cmd) { return addInterface(cmd, false); } string MultiEventGenerator::addRndInterface(string cmd) { return addInterface(cmd, true); } string MultiEventGenerator::addInterface(string cmd, bool rnd) { breakThePEG(); string noun = StringUtils::car(cmd); IBPtr ip = BaseRepository::getObjectFromNoun(noun); const InterfaceBase * ifb = BaseRepository:: FindInterface(ip, BaseRepository::getInterfaceFromNoun(noun)); string posarg = BaseRepository::getPosArgFromNoun(noun); cmd = StringUtils::cdr(cmd); if ( cmd.empty() ) return "Error: empty argument list."; string ret; string oldvalue = ifb->exec(*ip, "get", posarg); StringVector args; try { if ( rnd ) { do { args.push_back(StringUtils::car(cmd)); cmd = StringUtils::cdr(cmd); } while ( !cmd.empty() ); if ( args.size() < 3 ) return "Error: Argument list should be 'N min max mean width'."; int N = atoi(args[0].c_str()); string vmin = args[1]; string vmax = args[2]; string vmean = "0"; if ( args.size() > 3 ) vmean = args[3]; string vwidth = "0"; if ( args.size() > 4 ) vwidth = args[4]; string arg = "RND " + vmin + " " + vmax + " " + vmean + " " + vwidth; args = vector(N, arg); ifb->exec(*ip, "set", vmin); ifb->exec(*ip, "set", vmax); ifb->exec(*ip, "set", posarg + " " + oldvalue); } else { do { args.push_back(StringUtils::car(cmd, ",")); cmd = StringUtils::cdr(cmd, ","); } while ( !cmd.empty() ); for ( string::size_type i = 0; i < args.size(); ++i ) ifb->exec(*ip, "set", args[i]); } } catch (const Exception & e) { e.handle(); ret = "Error: " + e.message(); } ifb->exec(*ip, "set", posarg + " " + oldvalue); if ( !ret.empty() ) return ret; for ( string::size_type i = 0; i < theObjects.size(); ++i ) { if ( theObjects[i] == ip && theInterfaces[i] == ifb->name() && thePosArgs[i] == posarg ) { if ( rnd || theValues[i][0].substr(0,3) == "RND" ) theValues[i] = args; else theValues[i].insert(theValues[i].end(), args.begin(), args.end()); return ""; } } theObjects.push_back(ip); theInterfaces.push_back(ifb->name()); thePosArgs.push_back(posarg); theValues.push_back(args); return ""; } void MultiEventGenerator::addTag(string tag) { if ( tag[0] == '#' ) { string::size_type dash = tag.find('-'); if ( dash == string::npos ) firstSubrun = lastSubrun = atoi(tag.substr(1).c_str()); else { firstSubrun = atoi(tag.substr(1, dash - 1).c_str()); lastSubrun = atoi(tag.substr(dash + 1).c_str()); } } EventGenerator::addTag(tag); } void MultiEventGenerator::doGo(long next, long maxevent, bool tics) { if ( theObjects.empty() || next < 0 ) { EventGenerator::doGo(next, maxevent, tics); return; } if ( maxevent >= 0 ) N(maxevent); vector interfaces; long nargs = 1; for ( string::size_type i = 0; i < theObjects.size(); ++i ) { nargs *= theValues[i].size(); interfaces.push_back(BaseRepository::FindInterface(theObjects[i], theInterfaces[i])); } if ( theSeparateRandom ) { theSeparateRandom->init(); theSeparateRandom->initrun(); const InterfaceBase * ifb = BaseRepository::FindInterface(theSeparateRandom, "Seed"); ifb->exec(*theSeparateRandom, "set", "0"); } openOutputFiles(); string baseName = runName(); if ( tics ) tic(next - 1, nargs*N()); for ( long iargs = 0; iargs < nargs; ++iargs ) { ostringstream subname; subname << baseName << ":" << iargs + 1; runName(subname.str()); string head = heading(iargs, interfaces, baseName); if ( ( firstSubrun > 0 && iargs + 1 < firstSubrun ) || ( lastSubrun > 0 && iargs + 1 > lastSubrun ) ) { if ( theSeparateRandom ) { // This is needed to ensure the same random settings for a // given sub-run irrespectively if previous sub-runs have been // included or not. theSeparateRandom->reset(); theSeparateRandom->init(); theSeparateRandom->initrun(); } continue; } log() << head; out() << head; reset(); - for_each(objects(), mem_fun(&InterfacedBase::reset)); + for_each(objects(), mem_fn(&InterfacedBase::reset)); init(); initrun(); ieve = next-1; try { while ( shoot() ) { if ( tics ) tic(ieve + iargs*N(), nargs*N()); } } catch ( ... ) { finish(); throw; } finish(); } runName(baseName); finally(); } string MultiEventGenerator:: heading(long iargs, const vector & interfaces, string baseName) const { ostringstream os; long div = 1; if ( iargs > 0 ) os << endl; os << ">> " << baseName << " sub-run number " << iargs + 1 << " using the following interface values:" << endl; for ( string::size_type i = 0; i < theObjects.size(); ++i ) { long iarg = (iargs/div)%theValues[i].size(); string sval = theValues[i][iarg]; if ( theValues[i][iarg].substr(0,3) == "RND" ) { double vmin, vmax, vmean, vwidth; istringstream is(theValues[i][iarg].substr(3)); is >> vmin >> vmax >> vmean >> vwidth; double val = randomArg().rnd(vmin, vmax); if ( vwidth > 0.0 ) do { val = randomArg().rndGauss(vwidth, vmean); } while ( val < vmin || val > vmax ); ostringstream ssv; ssv << val; sval = ssv.str(); } interfaces[i]->exec(*theObjects[i], "set", thePosArgs[i] + " " + sval); os << " set " << theObjects[i]->name() << ":" << theInterfaces[i]; if ( !thePosArgs[i].empty() ) os << "[" << thePosArgs[i] << "]"; os << " " << sval << endl; div *= theValues[i].size(); } os << endl; return os.str(); } void MultiEventGenerator::persistentOutput(PersistentOStream & os) const { os << theObjects << theInterfaces << thePosArgs << theValues << firstSubrun << lastSubrun << theSeparateRandom; } void MultiEventGenerator::persistentInput(PersistentIStream & is, int) { is >> theObjects >> theInterfaces >> thePosArgs >> theValues >> firstSubrun >> lastSubrun >> theSeparateRandom; } IVector MultiEventGenerator::getReferences() { IVector ret = EventGenerator::getReferences(); ret.insert(ret.end(), theObjects.begin(), theObjects.end()); return ret; } void MultiEventGenerator::rebind(const TranslationMap & trans) { for ( string::size_type i = 0; i < theObjects.size(); ++i ) theObjects[i] = trans.translate(theObjects[i]); EventGenerator::rebind(trans); } ClassDescription MultiEventGenerator::initMultiEventGenerator; // Definition of the static class description member. void MultiEventGenerator::Init() { static ClassDocumentation documentation ("The ThePEG::MultiEventGenerator class is derived from the " "ThePEG::EventGenerator and is capable of making " "several runs with a pre-defined set of parameter and switch values."); static Command interfaceAddInterface ("AddInterface", "If arguments are given on the form 'object-name:interface-name arg1, " "arg2, arg3' or 'object-name:vectorinterface-name[pos] arg1, arg2, arg3' " "the generator will be run three times with the corresonding interface of " "the given object set to arg1, arg2, arg3 in each run respectively. If " "another interface with e.g. 4 different arguments, the generator will " "be run 12 times once for each combination of arguments. If called with " "an object and interface wich has already been given in a previous call, " "the new arguments will be added to the previously specified list without " "checking if any argument is doubled.", &MultiEventGenerator::addInterface); static Command interfaceAddRndInterface ("AddRndInterface", "If arguments are given on the form 'object-name:interface-name N min max " "mean width'"" or 'object-name:vectorinterface-name[pos] N min max mean width' " "the generator will be run N times with the corresonding interface of " "the given object set to a random value between min and max according to a " "Gaussian distribution with the given mean and width (if the width is absent " "or zero a flat distribution between min and max will be used instead) . If " "another interface with e.g. 4 different arguments, the generator will " "be run N*4 times once for each combination of arguments and the specified " "interface will get a new random value each time.. If called with " "an object and interface wich has already been given in a previous call to " "AddInterface or AddRndInterface " "the previous call will be ignored.", &MultiEventGenerator::addRndInterface); static Command interfaceRemoveInterface ("RemoveInterface", "If arguments are given on the form 'object-name:interface-name' and " "the same interface and object was previously with an " "AddInterface}, the corresponding arguments are " "removed and the interfaced will be left unchanged during the generation.", &MultiEventGenerator::removeInterface); static Reference interfaceSeparateRandom ("SeparateRandom", "A separate random number generator used for AddRndInterface" " to ensure reproducible sequences of interface values. If null, the standard " "random generator will be used instead.", &MultiEventGenerator::theSeparateRandom, true, false, true, true, false); interfaceAddInterface.rank(10.7); interfaceRemoveInterface.rank(10.5); } diff --git a/Repository/Repository.cc b/Repository/Repository.cc --- a/Repository/Repository.cc +++ b/Repository/Repository.cc @@ -1,1131 +1,1131 @@ // -*- C++ -*- // // Repository.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 1999-2017 Leif Lonnblad // // ThePEG is licenced under version 3 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 Repository class. // // macro is passed in from -D compile flag #ifndef THEPEG_PKGLIBDIR #error Makefile.am needs to define THEPEG_PKGLIBDIR #endif #include "Repository.h" #include "ThePEG/Utilities/Rebinder.h" #include "ThePEG/Handlers/EventHandler.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Repository/Strategy.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Utilities/Debug.h" #include "ThePEG/Config/algorithm.h" #include "ThePEG/Utilities/DynamicLoader.h" #include "ThePEG/Utilities/StringUtils.h" #include #include #include // readline options taken from // http://autoconf-archive.cryp.to/vl_lib_readline.html // Copyright © 2008 Ville Laurikari // Copying and distribution of this file, with or without // modification, are permitted in any medium without royalty provided // the copyright notice and this notice are preserved. #ifdef HAVE_LIBREADLINE # if defined(HAVE_READLINE_READLINE_H) # include # elif defined(HAVE_READLINE_H) # include # else extern "C" char *readline (const char *); # endif #endif #ifdef HAVE_READLINE_HISTORY # if defined(HAVE_READLINE_HISTORY_H) # include # elif defined(HAVE_HISTORY_H) # include # else extern "C" void add_history (const char *); # endif #endif using namespace ThePEG; ParticleMap & Repository::defaultParticles() { static ParticleMap theMap; return theMap; } ParticleDataSet & Repository::particles() { static ParticleDataSet theSet; return theSet; } MatcherSet & Repository::matchers() { static MatcherSet theSet; return theSet; } Repository::GeneratorMap & Repository::generators() { static GeneratorMap theMap;; return theMap; } string & Repository::currentFileName() { static string theCurrentFileName; return theCurrentFileName; } int & Repository::exitOnError() { static int exitonerror = 0; return exitonerror; } void Repository::cleanup() { generators().clear(); } void Repository::Register(IBPtr ip) { BaseRepository::Register(ip); registerParticle(dynamic_ptr_cast(ip)); registerMatcher(dynamic_ptr_cast(ip)); } void Repository::Register(IBPtr ip, string newName) { DirectoryAppend(newName); BaseRepository::Register(ip, newName); registerParticle(dynamic_ptr_cast(ip)); registerMatcher(dynamic_ptr_cast(ip)); } void Repository::registerParticle(tPDPtr pd) { if ( !pd ) return; if ( !member(particles(), pd) ) { particles().insert(pd); CreateDirectory(pd->fullName()); } if ( pd->id() == 0 ) return; if ( !member(defaultParticles(), pd->id()) ) defaultParticles()[pd->id()] = pd; for ( MatcherSet::iterator it = matchers().begin(); it != matchers().end(); ++it) (*it)->addPIfMatch(pd); } void Repository::registerMatcher(tPMPtr pm) { if ( !pm || member(matchers(), pm) ) return; pm->addPIfMatchFrom(particles()); for ( MatcherSet::iterator it = matchers().begin(); it != matchers().end(); ++it) { (*it)->addMIfMatch(pm); pm->addMIfMatch(*it); } matchers().insert(pm); } tPDPtr Repository::findParticle(string name) { tPDPtr pd; string path = name; DirectoryAppend(path); pd = dynamic_ptr_cast(GetPointer(path)); if ( pd ) return pd; for ( ParticleMap::iterator pit = defaultParticles().begin(); pit != defaultParticles().end(); ++pit ) if ( pit->second->PDGName() == name ) return pit->second; for ( ParticleDataSet::iterator pit = particles().begin(); pit != particles().end(); ++pit ) if ( (**pit).PDGName() == name ) return *pit; return pd; } tPMPtr Repository::findMatcher(string name) { for ( MatcherSet::iterator mit = matchers().begin(); mit != matchers().end(); ++mit ) if ( name == (**mit).name() ) return *mit; return tPMPtr(); } void Repository::saveRun(string EGname, string name, string filename) { EGPtr eg = BaseRepository::GetObject(EGname); EGPtr run = makeRun(eg, name); PersistentOStream os(filename, globalLibraries()); if ( ThePEG_DEBUG_ITEM(3) ) clog() << "Saving event generator '" << name << "'... " << flush; os << run; if ( ThePEG_DEBUG_ITEM(3) ) clog() << "done" << endl; } EGPtr Repository::makeRun(tEGPtr eg, string name) { // Clone all objects relevant for the EventGenerator. This is // the EventGenerator itself, all particles and all particle // matchers. 'localObject' is the set of all object refered to by // the generator particles and matcher and in the end these are // cloned as well. // Clone all Particle matchers if ( ThePEG_DEBUG_ITEM(3) ) clog() << "Making event generator '" << name << "':" << endl << "Updating all objects... " << flush; if ( ThePEG_DEBUG_ITEM(3) ) clog() << "done\nCloning matchers and particles... " << flush; MatcherSet localMatchers; ObjectSet localObjects; ObjectSet clonedObjects; TranslationMap trans; for ( MatcherSet::iterator mit = matchers().begin(); mit != matchers().end(); ++mit ) { PMPtr pm = clone(**mit); pm->clear(); trans[*mit] = pm; localMatchers.insert(pm); clonedObjects.insert(pm); localObjects.insert(*mit); addReferences(*mit, localObjects); } // Clone the particles. But only the ones which should be // used. First select the localParticles of the EventGenerator, then // add particles from the strategy of the EventGenerator which have // not already been selected. Finally add particles from the global // default if no default directories has been specified in the // strategy which have not already been selected. PDVector allParticles; for ( ParticleMap::const_iterator pit = eg->localParticles().begin(); pit != eg->localParticles().end(); ++pit ) allParticles.push_back(pit->second); if ( eg->strategy() ) { tcStrategyPtr strat = eg->strategy(); for ( ParticleMap::const_iterator pit = strat->particles().begin(); pit != strat->particles().end(); ++pit ) allParticles.push_back(pit->second); vector pdirs; if ( eg->strategy()->localParticlesDir().length() ) pdirs.push_back(eg->strategy()->localParticlesDir()); pdirs.insert(pdirs.end(), eg->strategy()->defaultParticlesDirs().begin(), eg->strategy()->defaultParticlesDirs().end()); for ( int i = 0, N = pdirs.size(); i < N; ++i ) { string dir = pdirs[i]; for ( ParticleDataSet::iterator pit = particles().begin(); pit != particles().end(); ++pit ) if ( (**pit).fullName().substr(0, dir.length()) == dir ) allParticles.push_back(*pit); } } if ( !eg->strategy() || eg->strategy()->defaultParticlesDirs().empty() ) for ( ParticleMap::iterator pit = defaultParticles().begin(); pit != defaultParticles().end(); ++pit ) allParticles.push_back(pit->second); for ( ParticleDataSet::iterator pit = particles().begin(); pit != particles().end(); ++pit ) allParticles.push_back(*pit); ParticleMap localParticles; set pdgnames; for ( PDVector::iterator pit = allParticles.begin(); pit != allParticles.end(); ++pit ) { ParticleMap::iterator it = localParticles.find((**pit).id()); if ( it == localParticles.end() ) { PDPtr pd = clone(**pit); trans[*pit] = pd; localParticles[pd->id()] = pd; clonedObjects.insert(pd); localObjects.insert(*pit); addReferences(*pit, localObjects); if ( pdgnames.find(pd->PDGName()) != pdgnames.end() ) std::cerr << "Using duplicate PDGName " << pd->PDGName() << " for a new particle.\n This can cause problems and is not " << "recommended.\n If this second particle is a new particle " << "in a BSM Model we recommend you change the name of the particle.\n"; else pdgnames.insert(pd->PDGName()); } else { trans[*pit] = it->second; } } if ( ThePEG_DEBUG_ITEM(3) ) clog() << "done\nCloning other objects... " << flush; // Clone the OldEventGenerator object to be used: localObjects.insert(eg); addReferences(eg, localObjects); EGPtr egrun = clone(*eg); clonedObjects.insert(egrun); trans[eg] = egrun; for ( ObjectSet::iterator it = localObjects.begin(); it != localObjects.end(); ++it ) { if ( member(trans.map(), *it) ) continue; IBPtr ip = clone(**it); trans[*it] = ip; clonedObjects.insert(ip); } if ( ThePEG_DEBUG_ITEM(3) ) clog() << "done\nRebind references... " << flush; IVector defaults; trans.translate(inserter(defaults), eg->defaultObjects().begin(), eg->defaultObjects().end()); if ( eg->strategy() ) trans.translate(inserter(defaults), eg->strategy()->defaultObjects().begin(), eg->strategy()->defaultObjects().end()); for ( ObjectSet::iterator it = clonedObjects.begin(); it != clonedObjects.end(); ++it ) { dynamic_cast(**it).theGenerator = egrun; rebind(**it, trans, defaults); } // Now, dependencies may have changed, so we do a final round of // updates. if ( ThePEG_DEBUG_ITEM(3) ) clog() << "done\nUpdating cloned objects... " << flush; if ( ThePEG_DEBUG_ITEM(3) ) clog() << "done\nInitializing... " << flush; clonedObjects.erase(egrun); egrun->setup(name, clonedObjects, localParticles, localMatchers); if ( ThePEG_DEBUG_ITEM(3) ) clog() << "done" << endl; generators()[name] = egrun; return egrun; } PDPtr Repository::defaultParticle(PID id) { ParticleMap::iterator pit = defaultParticles().find(id); return pit == defaultParticles().end()? PDPtr(): pit->second; } void Repository::defaultParticle(tPDPtr pdp) { if ( pdp ) defaultParticles()[pdp->id()] = pdp; } struct ParticleOrdering { bool operator()(tcPDPtr p1, tcPDPtr p2) const { 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 MatcherOrdering { bool operator()(tcPMPtr m1, tcPMPtr m2) const { return m1->name() < m2->name() || ( m1->name() == m2->name() && m1->fullName() < m2->fullName() ); } }; struct InterfaceOrdering { bool operator()(tcIBPtr i1, tcIBPtr i2) const { return i1->fullName() < i2->fullName(); } }; void Repository::save(string filename) { if ( ThePEG_DEBUG_ITEM(3) ) clog() << "saving '" << filename << "'... " << flush; PersistentOStream os(filename, globalLibraries()); set part(particles().begin(), particles().end()); set match(matchers().begin(), matchers().end()); os << objects().size(); for ( ObjectMap::iterator it = objects().begin(); it != objects().end(); ++it ) os << it->second; os << defaultParticles() << part << match << generators() << directories() << directoryStack() << globalLibraries() << readDirs(); if ( ThePEG_DEBUG_ITEM(3) ) clog() << "(" << objects().size() << " objects in " << directories().size() << " directories) done" << endl; } string Repository::load(string filename) { if ( ThePEG_DEBUG_ITEM(3) ) clog() << "loading '" << filename << "'... " << flush; currentFileName() = filename; PersistentIStream * is = new PersistentIStream(filename); if ( !*is ) { delete is; // macro is passed in from -D compile flag string fullpath = string(THEPEG_PKGLIBDIR) + '/' + filename; is = new PersistentIStream(fullpath); if ( !*is ) { delete is; return "Error: Could not find repository '" + filename + "'."; } } *is >> allObjects() >> defaultParticles() >> particles() >> matchers() >> generators() >> directories() >> directoryStack() >> globalLibraries() >> readDirs(); delete is; objects().clear(); for ( ObjectSet::iterator it = allObjects().begin(); it != allObjects().end(); ++it ) objects()[(**it).fullName()] = *it; if ( ThePEG_DEBUG_ITEM(3) ) clog() << "(" << objects().size() << " objects in " << directories().size() << " directories) done\nUpdating... " << flush; BaseRepository::resetAll(allObjects()); BaseRepository::update(); if ( ThePEG_DEBUG_ITEM(3) ) clog() << "done" << endl; return ""; } void Repository::stats(ostream & os) { os << "number of objects: " << setw(6) << objects().size() << endl; os << "number of objects (all): " << setw(6) << allObjects().size() << endl; os << "number of particles: " << setw(6) << particles().size() << endl; os << "number of matchers: " << setw(6) << matchers().size() << endl; } string Repository::read(string filename, ostream & os) { ifstream is; string file = filename; if ( file[0] == '/' ) { if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= trying to open " << file << " =)" << endl; is.open(file.c_str()); } else { vector dirs(readDirs().rbegin(), readDirs().rend()); dirs.push_back(currentReadDirStack().top()); if ( ThePEG_DEBUG_LEVEL > 1 ) { os << "(= search path order =)\n(== "; std::copy(dirs.rbegin(), dirs.rend(), std::ostream_iterator(os, " ==)\n(== ")); os << ")" << endl; } while ( dirs.size() ) { string dir = dirs.back(); if ( dir != "" && dir[dir.length() -1] != '/' ) dir += '/'; file = dir + filename; is.clear(); if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= trying to open " << file << " =)" << endl; is.open(file.c_str()); if ( is ) break; if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= no, try next search path =)" << endl; dirs.pop_back(); } } if ( !is ) { return "Error: Could not find input file '" + filename + "'"; } if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= yes =)" << endl; const string dir = StringUtils::dirname(file); if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= pushing <" << dir << "> to stack =)" << endl; currentReadDirStack().push(dir); try { Repository::read(is, os); if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= popping <" << currentReadDirStack().top() << "> from stack =)" << endl; currentReadDirStack().pop(); } catch ( ... ) { if ( ThePEG_DEBUG_LEVEL > 1 ) os << "(= popping <" << currentReadDirStack().top() << "> from stack =)" << endl; currentReadDirStack().pop(); throw; } return ""; } string Repository:: modifyEventGenerator(EventGenerator & eg, string filename, ostream & os, bool initOnly) { ObjectSet objs = eg.objects(); objs.insert(&eg); for ( ObjectSet::iterator it = objs.begin(); it != objs.end(); ++it ) { string name = (**it).fullName(); if ( name.rfind('/') != string::npos ) CreateDirectory(name.substr(0, name.rfind('/') + 1)); objects()[name] = *it; allObjects().insert(*it); } string msg = read(filename, os); if ( !msg.empty() ) return msg; - for_each(objs, mem_fun(&InterfacedBase::reset)); + for_each(objs, mem_fn(&InterfacedBase::reset)); eg.initialize(initOnly); if ( !generators().empty() ) msg += "Warning: new generators were initialized while modifying " + eg.fullName() + ".\n"; return msg; } void Repository::resetEventGenerator(EventGenerator & eg) { ObjectSet objs = eg.objects(); objs.insert(&eg); for ( ObjectSet::iterator it = objs.begin(); it != objs.end(); ++it ) { string name = (**it).fullName(); if ( name.rfind('/') != string::npos ) CreateDirectory(name.substr(0, name.rfind('/') + 1)); objects()[name] = *it; allObjects().insert(*it); } - for_each(objs, mem_fun(&InterfacedBase::reset)); + for_each(objs, mem_fn(&InterfacedBase::reset)); eg.initialize(true); } void Repository::execAndCheckReply(string line, ostream & os) { string reply = exec(line, os); if ( reply.size() ) os << reply; if ( reply.size() && reply[reply.size()-1] != '\n' ) os << endl; if ( exitOnError() && reply.size() >= 7 && reply.substr(0, 7) == "Error: " ) exit(exitOnError()); } void Repository::read(istream & is, ostream & os, string prompt) { #ifdef HAVE_LIBREADLINE if ( &is == &std::cin ) { char * line_read = 0; do { if ( line_read ) { free(line_read); line_read = 0; } line_read = readline(prompt.c_str()); if ( line_read && *line_read ) { string line = line_read; while ( !line.empty() && line[line.size() - 1] == '\\' ) { line[line.size() - 1] = ' '; char * cont_read = readline("... "); if ( cont_read ) { line += cont_read; free(cont_read); } } if ( prompt.empty() && ThePEG_DEBUG_LEVEL > 0 ) os << "(" << line << ")" << endl; #ifdef HAVE_READLINE_HISTORY add_history(line.c_str()); #endif // HAVE_READLINE_HISTORY execAndCheckReply(line, os); } } while ( line_read ); } else { #endif // HAVE_LIBREADLINE string line; if ( prompt.size() ) os << prompt; while ( getline(is, line) ) { while ( !line.empty() && line[line.size() - 1] == '\\' ) { line[line.size() - 1] = ' '; string cont; if ( prompt.size() ) os << "... "; getline(is, cont); line += cont; } if ( prompt.empty() && ThePEG_DEBUG_LEVEL > 0 ) os << "(" << line << ")" << endl; execAndCheckReply(line, os); if ( prompt.size() ) os << prompt; } #ifdef HAVE_LIBREADLINE } #endif if ( prompt.size() ) os << endl; } string Repository::copyParticle(tPDPtr p, string newname) { DirectoryAppend(newname); string newdir = newname.substr(0, newname.rfind('/')+1); newname =newname.substr(newname.rfind('/')+1); if ( newname.empty() ) newname = p->name(); if ( GetPointer(newdir + newname) ) return "Error: Cannot create particle " + newdir + newname + ". Object already exists."; if ( p->CC() && GetPointer(newdir + p->CC()->name()) ) return "Error: Cannot create anti-particle " + newdir + newname + ". Object already exists."; PDPtr pd = p->pdclone(); Register(pd, newdir + newname); pd->theDecaySelector.clear(); pd->theDecayModes.clear(); pd->isStable = true; if ( p->CC() ) { PDPtr apd = p->CC()->pdclone(); Register(apd, newdir + apd->name()); apd->theDecaySelector.clear(); apd->theDecayModes.clear(); apd->isStable = true; pd->theAntiPartner = apd; apd->theAntiPartner = pd; pd->syncAnti = p->syncAnti; apd->syncAnti = p->CC()->syncAnti; } HoldFlag<> dosync(pd->syncAnti, true); for ( DecaySet::const_iterator it = p->theDecayModes.begin(); it != p->theDecayModes.end(); ++it ) pd->addDecayMode(*it); return ""; } void Repository::remove(tIBPtr ip) { ObjectMap::iterator it = objects().find(ip->fullName()); if ( it == objects().end() || ip != it->second ) return; objects().erase(it); allObjects().erase(ip); if ( dynamic_ptr_cast(ip) ) { particles().erase(dynamic_ptr_cast(ip)); defaultParticles().erase(dynamic_ptr_cast(ip)->id()); } if ( dynamic_ptr_cast(ip) ) matchers().erase(dynamic_ptr_cast(ip)); } string Repository::remove(const ObjectSet & rmset) { ObjectSet refset; for ( ObjectMap::const_iterator i = objects().begin(); i != objects().end(); ++i ) { if ( member(rmset, i->second) ) continue; IVector ov = DirectReferences(i->second); for ( int j = 0, M = ov.size(); j < M; ++j ) if ( member(rmset, ov[j]) ) { refset.insert(i->second); break; } } if ( refset.empty() ) { for ( ObjectSet::iterator oi = rmset.begin(); oi != rmset.end(); ++oi ) remove(*oi); return ""; } string ret = "Error: cannot remove the objects because the following " "objects refers to some of them:\n"; for ( ObjectSet::iterator oi = refset.begin(); oi != refset.end(); ++oi ) ret += (**oi).fullName() + "\n"; return ret; } string Repository::exec(string command, ostream & os) { string cpcmd = command; try { string verb = StringUtils::car(command); command = StringUtils::cdr(command); if ( verb == "help" ) { help(command, os); return ""; } if ( verb == "rm" ) { ObjectSet rmset; while ( !command.empty() ) { string name = StringUtils::car(command); DirectoryAppend(name); IBPtr obj = GetPointer(name); if ( !obj ) return "Error: Could not find object named " + name; rmset.insert(obj); command = StringUtils::cdr(command); } return remove(rmset); } if ( verb == "rmdir" || verb == "rrmdir" ) { string dir = StringUtils::car(command); DirectoryAppend(dir); if ( dir[dir.size() - 1] != '/' ) dir += '/'; if ( !member(directories(), dir) ) return verb == "rmdir"? "Error: No such directory.": ""; IVector ov = SearchDirectory(dir); if ( ov.size() && verb == "rmdir" ) return "Error: Cannot remove a non-empty directory. " "(Use rrmdir do remove all object and subdirectories.)"; ObjectSet rmset(ov.begin(), ov.end()); string ret = remove(rmset); if ( !ret.empty() ) return ret; StringVector dirs(directories().begin(), directories().end()); for ( int i = 0, N = dirs.size(); i < N; ++ i ) if ( dirs[i].substr(0, dir.size()) == dir ) directories().erase(dirs[i]); for ( int i = 0, N = directoryStack().size(); i < N; ++i ) if ( directoryStack()[i].substr(0, dir.size()) == dir ) directoryStack()[i] = '/'; return ""; } if ( verb == "cp" ) { string name = StringUtils::car(command); DirectoryAppend(name); tPDPtr p = dynamic_ptr_cast(GetPointer(name)); if ( p ) return copyParticle(p, StringUtils::cdr(command)); return BaseRepository::exec(cpcmd, os); } if ( verb == "setup" ) { string name = StringUtils::car(command); DirectoryAppend(name); IBPtr obj = GetPointer(name); if ( !obj ) return "Error: Could not find object named " + name; istringstream is(StringUtils::cdr(command)); readSetup(obj, is); // A particle may have been registered before but under the wrong id(). PDPtr pd = dynamic_ptr_cast(obj); if(pd) registerParticle(pd); return ""; } if ( verb == "decaymode" ) { string tag = StringUtils::car(command); DMPtr dm = DecayMode::constructDecayMode(tag); if ( !dm ) return "Error: Could not create decay mode from the tag " + StringUtils::car(command); istringstream is(StringUtils::cdr(command)); readSetup(dm, is); if ( !dm->CC() ) return ""; if ( dm->CC()->parent()->synchronized() ) { dm->CC()->synchronize(); return ""; } if ( !dm->CC()->decayer() ) return FindInterface(dm, "Decayer")-> exec(*dm->CC(), "set", dm->decayer()->fullName()); return ""; } if ( verb == "makeanti" ) { string name = StringUtils::car(command); DirectoryAppend(name); tPDPtr p = dynamic_ptr_cast(GetPointer(name)); if ( !p ) return "Error: No particle named " + name; name = StringUtils::car(StringUtils::cdr(command)); DirectoryAppend(name); tPDPtr ap = dynamic_ptr_cast(GetPointer(name)); if ( !ap ) return "Error: No particle named " + name; ParticleData::antiSetup(PDPair(p, ap)); return ""; } if ( verb == "read" ) { // remember directory we're in string cwd = directoryStack().back(); string filename = StringUtils::car(command); string msg = read(filename, os); // Return to the original directory, so that // calling 'read' in an input file will not change the // repository directory you're in ChangeDirectory(cwd); return msg; } if ( verb == "load" ) { return load(StringUtils::car(command)); } if ( verb == "save" ) { save(StringUtils::car(command)); return ""; } if ( verb == "lsruns" ) { string ret; for ( GeneratorMap::iterator ieg = generators().begin(); ieg != generators().end(); ++ieg ) ret += ieg->first + "\n"; return ret; } if ( verb == "makerun" ) { string runname = StringUtils::car(command); string generator = StringUtils::car(StringUtils::cdr(command)); DirectoryAppend(generator); EGPtr eg = BaseRepository::GetObject(generator); makeRun(eg, runname); return ""; } if ( verb == "rmrun" ) { string runname = StringUtils::car(command); generators().erase(runname); return ""; } if ( verb == "saverun" || verb == "saverunfile" || verb == "run" ) { string runname = StringUtils::car(command); string generator = StringUtils::car(StringUtils::cdr(command)); DirectoryAppend(generator); GeneratorMap::iterator ieg = generators().find(runname); EGPtr eg; if ( ieg == generators().end() ) { eg = BaseRepository::GetObject(generator); eg = makeRun(eg, runname); } else eg = ieg->second; if ( !eg ) return "Error: Could not create/find run named'" + runname + "'."; if ( verb == "run" ) eg->go(); else if ( verb == "saverunfile" ) { string file = generator; PersistentOStream os(file, globalLibraries()); os << eg; if ( !os ) return "Save failed! (I/O error)"; } else { string file = eg->filename() + ".run"; PersistentOStream os(file, globalLibraries()); os << eg; if ( !os ) return "Save failed! (I/O error)"; } return ""; } if ( verb == "removerun" ) { string runname = StringUtils::car(command); GeneratorMap::iterator ieg = generators().find(runname); if ( ieg != generators().end() ) { generators().erase(ieg); return ""; } else return "Error: No run named '" + runname + "' available."; } if ( verb == "create" ) { string className = StringUtils::car(command); command = StringUtils::cdr(command); string name = StringUtils::car(command); const ClassDescriptionBase * db = DescriptionList::find(className); command = StringUtils::cdr(command); while ( !db && command.length() ) { string library = StringUtils::car(command); command = StringUtils::cdr(command); DynamicLoader::load(library); db = DescriptionList::find(className); } if ( !db ) { string msg = "Error: " + className + ": No such class found."; if ( !DynamicLoader::lastErrorMessage.empty() ) msg += "\nerror message from dynamic loader:\n" + DynamicLoader::lastErrorMessage; return msg; } IBPtr obj = dynamic_ptr_cast(db->create()); if ( !obj ) return "Error: Could not create object of this class class."; if ( name.empty() ) return "Error: No name specified."; Register(obj, name); return ""; } if ( verb == "defaultparticle" ) { while ( !command.empty() ) { string name = StringUtils::car(command); DirectoryAppend(name); tPDPtr p = dynamic_ptr_cast(GetPointer(name)); if ( !p ) return "Error: No particle named " + name; defaultParticle(p); command = StringUtils::cdr(command); } return ""; } if ( verb == "EXITONERROR" ) { exitOnError() = 1; return ""; } } catch (const Exception & e) { e.handle(); return "Error: " + e.message(); } return BaseRepository::exec(cpcmd, os); } void Repository::help(string cmd, ostream & os) { cmd = StringUtils::car(cmd); if ( cmd == "cd" ) os << "Usage: cd " << endl << "Set the current directory to ." << endl; else if ( cmd == "mkdir" ) os << "Usage: mkdir " << endl << "Create a new directory called with the given path name." << endl; else if ( cmd == "rmdir" ) os << "Usage: rmdir " << endl << "Remove an empty directory." << endl; else if ( cmd == "rrmdir" ) os << "Usage: rrmdir " << endl << "Remove a directory and everything that is in it recursively." << endl << "Will only succeed if no other objects refers to the ones to " << "be deleted." << endl; else if ( cmd == "cp" ) os << "Usage: cp " << endl << "Copy the given object to a new object with the given name." << endl; else if ( cmd == "setup" ) os << "Usage: setup ..." << endl << "Tell a given object to read information given by the arguments." << endl; else if ( cmd == "decaymode" ) os << "Usage: decaymode " << endl << "Construct a decay mode from the given decay tag. The resulting " << "object will be inserted in the directory with the same path as " << "the decaying particle object. The given brancing fraction will " << "be set as well as the given decayer object. If the mode should " << "be switched on by default 1(on) should be specified (otherwise " << "0(off))." << endl; else if ( cmd == "makeanti" ) os << "Usage: makeanti " << endl << "Indicate that the two given particle objects are eachothers " << "anti-partnets." << endl; else if ( cmd == "read" ) os << "Usage: read " << endl << "Read more commands from the given file. The file name can be " << "given relative to the current directory in the shell, or " << "relative to standard directories, or as an absolute path." << endl; else if ( cmd == "load" ) os << "Usage: load " << endl << "Discard everything in the reopsitory and read in a completely " << "new repository from the given file." << endl; else if ( cmd == "save" ) os << "Usage: save " << endl << "Save the complete repository to the given file." << endl; else if ( cmd == "lsruns" ) os << "Usage: lsruns" << endl << "List the run names of all initialized event generators." << endl; else if ( cmd == "makerun" ) os << "Usage: makerun " << endl << "Initialize the given event generator and assign a run name." << endl; else if ( cmd == "rmrun" ) os << "Usage: rmrun " << endl << "Remove the initialized event generator given by the run name." << endl; else if ( cmd == "saverun" ) os << "Usage: saverun " << endl << "Initialize the given event generator and assign a run name " << "and save it to a file named .run" << endl; else if ( cmd == "run" ) os << "Usage: run " << endl << "Run the initialized event generator given b the run name." << endl; else if ( cmd == "create" ) os << "Usage: create {}" << endl << "Create an object of the given class and assign the given name. " << "Optionally supply a dynamically loaded library where the class " << "is included." << endl; else if ( cmd == "pushd" ) os << "Usage: pushd " << endl << "Set the current directory to , but keep the previous " << "working directory on the directory stack." << endl; else if ( cmd == "popd" ) os << "Usage: popd" << endl << "Leave the current working directory and set the current " << "directory to the previous one on the directory stack." << endl; else if ( cmd == "pwd" ) os << "Usage: pwd" << endl << "Print the current working directory." << endl; else if ( cmd == "dirs" ) os << "Usage: dirs" << endl << " Print the contents of the directory stack." << endl; else if ( cmd == "mv" ) os << "Usage: mv " << endl << "Rename the given object to a new path name." << endl; else if ( cmd == "ls" ) os << "Usage: ls {}" << endl << "List the objects and subdirectories in the current or given " << "directory." << endl; else if ( cmd == "library" ) os << "Usage: library " << endl << "Make new classes available to the repository by dynamically " << "linking the given library." << endl; else if ( cmd == "globallibrary" ) os << "Usage: globallibrary " << endl << "Make new classes available to the repository by dynamically " << "linking the given library. If this repository is saved and read " << "in again, this library will be linked in from the beginning." << endl; else if ( cmd == "rmgloballibrary" ) os << "Usage: rmgloballibrary " << endl << "Remove a dynamic library previously added with globallibrary." << endl; else if ( cmd == "appendpath" ) os << "Usage: appendpath " << endl << "Add a search path for dynamic libraries to the end of the " << "search list." << endl; else if ( cmd == "lspaths" ) os << "Usage: lspaths" << endl << "List search paths for dynamic libraries." << endl; else if ( cmd == "prependpath" ) os << "Usage: prependpath " << endl << "Add a search path for dynamic libraries to the beginning of the " << "search list." << endl; else if ( cmd == "doxygendump" ) os << "Usage: doxygendump " << endl << "Extract doxygen documentation of all loaded classes in the " << "given name space and weite it to a file.." << endl; else if ( cmd == "mset" || cmd == "minsert" || cmd == "mdo" ) os << "Usage: " << cmd << " " << endl << "Recursively find in the given directory all objects of the " << "given class and call '" << cmd.substr(1) << "' with the given value for the given interface." << endl; else if ( cmd == "msetdef" || cmd == "mget" || cmd == "mdef" || cmd == "mmin" || cmd == "mmax" || cmd == "merase" ) os << "Usage: " << cmd << " " << endl << "Recursively find in the given directory all objects of the given " << "class and call '" << cmd.substr(1) << "' for the given interface." << endl; else if ( cmd == "set" ) os << "Usage: set : " << endl << "Set the interface for the given object to the given value." << endl; else if ( cmd == "setdef" ) os << "Usage: setdef :" << endl << "Set the interface for the given object to its default value." << endl; else if ( cmd == "insert" ) os << "Usage: insert : " << endl << "Insert a value in the vector interface of the given object." << endl; else if ( cmd == "erase" ) os << "Usage: erase :" << endl << "Erase a value from the vector interface of the given object." << endl; else if ( cmd == "do" ) os << "Usage: do : " << endl << "Call the command interface of the given object with the " << "given arguments." << endl; else if ( cmd == "get" ) os << "Usage: get :" << endl << "Print the value of the interface of the given object." << endl; else if ( cmd == "def" ) os << "Usage: def :" << endl << "Print the default value of the interface of the given object." << endl; else if ( cmd == "min" ) os << "Usage: min :" << endl << "Print the minimum value of the interface of the given object." << endl; else if ( cmd == "max" ) os << "Usage: max :" << endl << "Print the maximum value of the interface of the given object." << endl; else if ( cmd == "describe" ) os << "Usage: describe {:}" << endl << "Describe the given object or an interface of the object." << endl; else if ( cmd == "lsclass" ) os << "Usage: lsclass" << endl << "List all classes available in the repository." << endl; else if ( cmd == "all" ) { os << "Available commands:" << endl << "* cd, mkdir, rmdir, rrmdir, pwd, cp, mv, rm, pushd, popd, dirs, ls:\n" << " Manipulate the repository structure. Analogous to unix " << "shell commands." << endl << "* create, setup, decaymode makeanti:\n" << " Create or setup an object." << endl << "* set, get, insert, erase, do, detdef, def, min, max, describe\n" << " mset, minsert, mdo, msetdef, mdef, mmin, mmax, merase:\n" << " Manipulate interfaces to objects." << endl << "* makerun, saverun, run, lsruns, rmrun:\n" << " Create and handle initialized event genrators which can be run." << endl << "* read, load, library globallibrary, rmgloballibrary,\n" << " appendpath, prependpath, lspaths, doxygendump:\n" << " Handle files external files and libraries." << endl; os << "Do 'help syntax' for help on syntax." << endl << "Do 'help ' for help on a particular command." << endl; } else if ( cmd == "syntax" ) os << "* = '/' | | /" << endl << " = | / | :\n" << " Analogous to a unix file structure, an object can be " << "specified with an\n absolute path or a path relative to " << "the current directory." << endl << "* = |[]" << endl << " An interface can be a parameter (floating point, integer or " << "string),\n a switch (integer, possibly named), a reference to " << "another object in the\n repository or a command which takes " << "an arbitrary string as argument.\n There are also vector interfaces " << "of parameters and references for which\n an index must be supplied." << endl; else { if ( !cmd.empty() ) os << "No command '" << cmd << "' found." << endl; os << "Common commands:" << endl << "* cd, mkdir, rmdir, pwd, cp, mv, rm:\n" << " Manipulate the repository structure. Analogous to unix " << "shell commands." << endl << "* create, setup:\n" << " Create an object." << endl << "set, get, insert, erase, do:\n" << " Manipulate interfaces to objects." << endl << "* makerun, saverun, run, lsruns:\n" << " Create and handle initialized event genrators which can be run." << endl; os << "Do 'help all' for a complete list of commands." << endl << "Do 'help syntax' for help on syntax." << endl << "Do 'help ' for help on a particular command." << endl; } } Repository::Repository() { ++ninstances; } Repository::~Repository() { --ninstances; if ( ninstances <= 0 ) { generators().clear(); } } int Repository::ninstances = 0; namespace { static string version_ = #include "versionstamp.inc" ""; } string Repository::version() { return ::version_; } string Repository::banner() { const auto now = std::chrono::system_clock::now(); const auto now_c = std::chrono::system_clock::to_time_t(now); string time = ">>>> " ; time += StringUtils::stripws(string(std::ctime(&now_c))) + ' '; time += string(max(0,74 - int(time.size())), ' '); time += "<<<<"; string line = ">>>> Toolkit for HEP Event Generation - " + Repository::version() + ' '; line += string(max(0,78 - int(line.size())), '<'); string block = string(78, '>') + '\n' + line + '\n' + time + '\n' + string(78, '<') + '\n'; return block; }