optimizer/Drivers/pisa-opal/pisa-opal.cpp

@@ -33,6 +33,7 @@
 #include "Expression/RadialPeak.h"
 #include "Expression/SumErrSqRadialPeak.h"
 #include "Expression/MaxNormRadialPeak.h"
+#include "Expression/NumberOfPeaks.h"
 #include "Expression/ProbeVariable.h"
 //////////////////////////////////////////////////////////////////////////////
 
@@ -82,6 +83,9 @@ int main(int argc, char** argv) {
     ff = MaxNormRadialPeak();
     funcs.insert(std::pair<std::string, client::function::type>
             ("maxNormRadialPeak", ff));
+    ff = NumberOfPeaks();
+    funcs.insert(std::pair<std::string, client::function::type>
+            ("numberOfPeaks", ff));
     ff = SumErrSqRadialPeak();
     funcs.insert(std::pair<std::string, client::function::type>
             ("sumErrSqRadialPeak", ff));

optimizer/Expression/CMakeLists.txt

@@ -3,6 +3,7 @@ SET (_EXPR_SRCS
     SumErrSq.cpp
     FromFile.cpp
     MaxNormRadialPeak.cpp
+    NumberOfPeaks.h
     RadialPeak.cpp
     SumErrSqRadialPeak.cpp
     Parser/expression.cpp

optimizer/Expression/NumberOfPeaks.cpp

+#include "Expression/NumberOfPeaks.h"
+
+const std::string name = "NumberOfPeaks";

optimizer/Expression/NumberOfPeaks.h

+#ifndef __NUMBER_OF_PEAKS_H__
+#define __NUMBER_OF_PEAKS_H__
+
+#include <string>
+
+#include <cmath>
+
+#include "boost/type_traits/remove_cv.hpp"
+#include "boost/variant/get.hpp"
+#include "boost/variant/variant.hpp"
+#include "boost/smart_ptr.hpp"
+
+#include "Util/Types.h"
+#include "Util/PeakReader.h"
+#include "Expression/Parser/function.hpp"
+
+/**
+ *  A simple expression to check the number of turns in a circular machine. It checks
+ *  probe files (*.peaks) and counts the number of turns.
+ */
+struct NumberOfPeaks {
+
+    static const std::string name;
+
+    Expressions::Result_t operator()(client::function::arguments_t args) {
+
+        sim_filename_  = boost::get<std::string>(args[0]);
+
+        bool is_valid = true;
+
+        boost::scoped_ptr<PeakReader> sim_peaks(new PeakReader(sim_filename_));
+        std::size_t nPeaks = 0;
+        
+        try {
+            sim_peaks->parseFile();
+            nPeaks = sim_peaks->getNumberOfPeaks();
+        } catch (OptPilotException &ex) {
+            std::cout << "Caught exception: " << ex.what() << std::endl;
+            is_valid = false;
+        }
+
+        
+        return boost::make_tuple(nPeaks, is_valid);
+    }
+
+private:
+    std::string sim_filename_;
+    
+    // define a mapping to arguments in argument vector
+    boost::tuple<std::string> argument_types;
+};
+
+#endif
+// vi: set et ts=4 sw=4 sts=4:
+// Local Variables:
+// mode:c
+// c-basic-offset: 4
+// indent-tabs-mode:nil
+// End:

optimizer/Util/PeakReader.cpp

@@ -54,3 +54,8 @@ void PeakReader::getPeak(int nPeak, double& radius) {
                                 "peak not found!");
     }
 }
+
+
+std::size_t PeakReader::getNumberOfPeaks() {
+    return peaks_m.size();
+}
\ No newline at end of file

optimizer/Util/PeakReader.h

@@ -25,6 +25,11 @@ public:
      */
     void getPeak(int nPeak, double& radius);
     
+    /**
+     * @returns the number of peaks in the file
+     */
+    std::size_t getNumberOfPeaks();
+    
 private:
     /// Peak filename
     std::string filename_m;

src/Algorithms/bet/EnvelopeBunch.cpp

@@ -74,12 +74,12 @@ EnvelopeBunch::EnvelopeBunch(const PartData *ref):
 }
 
 
-EnvelopeBunch::EnvelopeBunch(const EnvelopeBunch &rhs):
-    PartBunch(rhs),
-    reference(rhs.reference),
-    numSlices_m(0),
-    numMySlices_m(0)
-{}
+// EnvelopeBunch::EnvelopeBunch(const EnvelopeBunch &rhs):
+//     PartBunch(rhs),
+//     reference(rhs.reference),
+//     numSlices_m(0),
+//     numMySlices_m(0)
+// {}
 
 
 EnvelopeBunch::EnvelopeBunch(const std::vector<OpalParticle> &rhs,

src/Algorithms/bet/EnvelopeBunch.h

@@ -58,7 +58,7 @@ enum EnvelopeBunchShape {
  * implementation
  */
 class EnvelopeBunch : public PartBunch {
-    
+
 public:
     /// Default constructor
     EnvelopeBunch(const PartData *ref);
@@ -67,8 +67,8 @@ public:
     EnvelopeBunch(const std::vector<OpalParticle> &,
                   const PartData *ref);
 
-    /// Copy constructor
-    EnvelopeBunch(const EnvelopeBunch &);
+    /// Copy constructor (implement if needed)
+    EnvelopeBunch(const EnvelopeBunch &) = delete;
 
     virtual ~EnvelopeBunch();
 

src/Classic/AbsBeamline/RFCavity.cpp

@@ -608,6 +608,7 @@ ElementBase::ElementType RFCavity::getType() const {
 }
 
 double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, double mass) {
+    using Physics::pi;
     const double dt = 1e-13;
     const double p0 = Util::getP(E0, mass);
     const double origPhase =getPhasem();

src/Classic/Algorithms/PartBunch.cpp

@@ -61,21 +61,13 @@ PartBunch::PartBunch(const PartData *ref): // Layout is set using setSolver()
 
 }
 
-PartBunch::PartBunch(const std::vector<OpalParticle> &rhs,
-                     const PartData *ref):
-    PartBunchBase<double, 3>(new PartBunch::pbase_t(new Layout_t()), rhs, ref),
-    interpolationCacheSet_m(false)
-{
-    ERRORMSG("should not be here: PartBunch::PartBunch(const std::vector<OpalParticle> &rhs, const PartData *ref):" << endl);
-}
-
-PartBunch::PartBunch(const PartBunch &rhs):
-    PartBunchBase<double, 3>(rhs),
-    interpolationCacheSet_m(rhs.interpolationCacheSet_m)
-{
-    ERRORMSG("should not be here: PartBunch::PartBunch(const PartBunch &rhs):" << endl);
-    std::exit(0);
-}
+// PartBunch::PartBunch(const PartBunch &rhs):
+//     PartBunchBase<double, 3>(rhs),
+//     interpolationCacheSet_m(rhs.interpolationCacheSet_m)
+// {
+//     ERRORMSG("should not be here: PartBunch::PartBunch(const PartBunch &rhs):" << endl);
+//     std::exit(0);
+// }
 
 
 PartBunch::~PartBunch() {

src/Classic/Algorithms/PartBunch.h

@@ -37,13 +37,12 @@ public:
 
     /// Default constructor.
     //  Construct empty bunch.
-    PartBunch(const PartData *ref);
+    explicit PartBunch(const PartData *ref);
 
-    /// Conversion.
-    PartBunch(const std::vector<OpalParticle> &,
-              const PartData *ref);
+    PartBunch() = delete;
+    PartBunch(const PartBunch &) = delete;
+    PartBunch &operator=(const PartBunch &) = delete;
 
-    PartBunch(const PartBunch &);
     ~PartBunch();
 
 //     pbase_t* clone();
@@ -125,10 +124,6 @@ private:
 
     ParticleAttrib<CacheDataCIC<double, 3U> > interpolationCache_m;
 
-
-
-    PartBunch &operator=(const PartBunch &) = delete;
-
     //FIXME
     ParticleLayout<double, 3> & getLayout() {
         return pbase->getLayout();

src/Classic/Algorithms/PartBunchBase.h

@@ -31,10 +31,10 @@ template <class T, unsigned Dim>
 class PartBunchBase
 {
 public:
-    typedef typename AbstractParticle<T, Dim>::ParticlePos_t ParticlePos_t;
+    typedef typename AbstractParticle<T, Dim>::ParticlePos_t   ParticlePos_t;
     typedef typename AbstractParticle<T, Dim>::ParticleIndex_t ParticleIndex_t;
-    typedef typename AbstractParticle<T, Dim>::UpdateFlags UpdateFlags;
-    typedef typename AbstractParticle<T, Dim>::Position_t Position_t;
+    typedef typename AbstractParticle<T, Dim>::UpdateFlags     UpdateFlags_t;
+    typedef typename AbstractParticle<T, Dim>::Position_t      Position_t;
 
     typedef std::pair<Vector_t, Vector_t> VectorPair_t;
 
@@ -44,22 +44,17 @@ public:
 
 public:
 
-    PartBunchBase(AbstractParticle<T, Dim>* pb);
+    explicit PartBunchBase(AbstractParticle<T, Dim>* pb);
 
     virtual ~PartBunchBase() { }
 
-    /*
-     * Bunch common member functions
-     */
-
     PartBunchBase(AbstractParticle<T, Dim>* pb, const PartData *ref);
 
-    /// Conversion.
-    PartBunchBase(AbstractParticle<T, Dim>* pb,
-                  const std::vector<OpalParticle> &,
-                  const PartData *ref); //TODO
+    PartBunchBase(const PartBunchBase &rhs) = delete; // implement if needed
 
-    PartBunchBase(const PartBunchBase &rhs); //TODO
+    /*
+     * Bunch common member functions
+     */
 
     // This is required since we initialize the Layout and the RegionLayout with default constructor
     virtual void initialize(FieldLayout_t *fLayout) = 0;
@@ -413,8 +408,8 @@ public:
     ParticleLayout<T, Dim> & getLayout();
     const ParticleLayout<T, Dim>& getLayout() const;
 
-    bool getUpdateFlag(UpdateFlags f) const;
-    void setUpdateFlag(UpdateFlags f, bool val);
+    bool getUpdateFlag(UpdateFlags_t f) const;
+    void setUpdateFlag(UpdateFlags_t f, bool val);
 
 
     ParticleBConds<Position_t, Dimension>& getBConds() {

src/Classic/Algorithms/PartBunchBase.hpp

@@ -10,69 +10,8 @@
 #include "Utilities/Options.h"
 #include "Utilities/Util.h"
 
-using Physics::pi;
-
 extern Inform *gmsg;
 
-// template <class T, unsigned Dim>
-// PartBunchBase<T, Dim>::PartBunchBase()
-//     : pbase(nullptr),
-//       myNode_m(Ippl::myNode()),
-//       nodes_m(Ippl::getNodes()),
-//       fixed_grid(false),
-//       pbin_m(nullptr),
-//       lossDs_m(nullptr),
-//       pmsg_m(nullptr),
-//       f_stream(nullptr),
-//       unit_state_(units),
-//       stateOfLastBoundP_(unitless),
-//       moments_m(),
-//       dt_m(0.0),
-//       t_m(0.0),
-//       eKin_m(0.0),
-//       dE_m(0.0),
-//       spos_m(0.0),
-//       rmax_m(0.0),
-//       rmin_m(0.0),
-//       rrms_m(0.0),
-//       prms_m(0.0),
-//       rmean_m(0.0),
-//       pmean_m(0.0),
-//       eps_m(0.0),
-//       eps_norm_m(0.0),
-//       rprms_m(0.0),
-//       Dx_m(0.0),
-//       Dy_m(0.0),
-//       DDx_m(0.0),
-//       DDy_m(0.0),
-//       hr_m(-1.0),
-//       nr_m(0),
-//       fs_m(nullptr),
-//       couplingConstant_m(0.0),
-//       qi_m(0.0),
-//       distDump_m(0),
-//       fieldDBGStep_m(0),
-//       dh_m(1e-12),
-//       tEmission_m(0.0),
-//       bingamma_m(nullptr),
-//       binemitted_m(nullptr),
-//       lPath_m(0.0),
-//       stepsPerTurn_m(0),
-//       localTrackStep_m(0),
-//       globalTrackStep_m(0),
-//       numBunch_m(1),
-//       SteptoLastInj_m(0),
-//       globalPartPerNode_m(nullptr),
-//       dist_m(nullptr),
-//       globalMeanR_m(Vector_t(0.0, 0.0, 0.0)),
-//       globalToLocalQuaternion_m(Quaternion_t(1.0, 0.0, 0.0, 0.0)),
-//       lowParticleCount_m(false),
-//       dcBeam_m(false)
-// {
-//     R(*(pbase->R_p));   // undefined behaviour due to reference to null pointer
-//     ID(*(pbase->ID_p));   // undefined behaviour due to reference to null pointer
-// }
-
 template <class T, unsigned Dim>
 PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb)
     : R(*(pb->R_p)),
@@ -163,230 +102,11 @@ PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb)
 
 template <class T, unsigned Dim>
 PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb, const PartData *ref)
-    : R(*(pb->R_p)),
-      ID(*(pb->ID_p)),
-      myNode_m(Ippl::myNode()),
-      nodes_m(Ippl::getNodes()),
-      fixed_grid(false),
-      pbin_m(nullptr),
-      lossDs_m(nullptr),
-      pmsg_m(nullptr),
-      f_stream(nullptr),
-      lowParticleCount_m(false),
-      reference(ref),
-      unit_state_(units),
-      stateOfLastBoundP_(unitless),
-      moments_m(),
-      dt_m(0.0),
-      t_m(0.0),
-      eKin_m(0.0),
-      dE_m(0.0),
-      spos_m(0.0),
-      globalMeanR_m(Vector_t(0.0, 0.0, 0.0)),
-      globalToLocalQuaternion_m(Quaternion_t(1.0, 0.0, 0.0, 0.0)),
-      rmax_m(0.0),
-      rmin_m(0.0),
-      rrms_m(0.0),
-      prms_m(0.0),
-      rmean_m(0.0),
-      pmean_m(0.0),
-      eps_m(0.0),
-      eps_norm_m(0.0),
-      rprms_m(0.0),
-      Dx_m(0.0),
-      Dy_m(0.0),
-      DDx_m(0.0),
-      DDy_m(0.0),
-      hr_m(-1.0),
-      nr_m(0),
-      fs_m(nullptr),
-      couplingConstant_m(0.0),
-      qi_m(0.0),
-      distDump_m(0),
-      fieldDBGStep_m(0),
-      dh_m(1e-12),
-      tEmission_m(0.0),
-      bingamma_m(nullptr),
-      binemitted_m(nullptr),
-      lPath_m(0.0),
-      stepsPerTurn_m(0),
-      localTrackStep_m(0),
-      globalTrackStep_m(0),
-      numBunch_m(1),
-      SteptoLastInj_m(0),
-      globalPartPerNode_m(nullptr),
-      dist_m(nullptr),
-      dcBeam_m(false),
-      periodLength_m(Physics::c / 1e9),
-      pbase(pb)
-{
-    setup(pb);
-
-    boundpTimer_m = IpplTimings::getTimer("Boundingbox");
-    boundpBoundsTimer_m = IpplTimings::getTimer("Boundingbox-bounds");
-    boundpUpdateTimer_m = IpplTimings::getTimer("Boundingbox-update");
-    statParamTimer_m = IpplTimings::getTimer("Compute Statistics");
-    selfFieldTimer_m = IpplTimings::getTimer("SelfField total");
-
-    histoTimer_m = IpplTimings::getTimer("Histogram");
-
-    distrCreate_m = IpplTimings::getTimer("Create Distr");
-    distrReload_m = IpplTimings::getTimer("Load Distr");
-
-    globalPartPerNode_m = std::unique_ptr<size_t[]>(new size_t[Ippl::getNodes()]);
-
-    lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(std::string("GlobalLosses"), !Options::asciidump));
-
-    pmsg_m.release();
-    //    f_stream.release();
-    /*
-      if(Ippl::getNodes() == 1) {
-          f_stream = std::unique_ptr<ofstream>(new ofstream);
-          f_stream->open("data/dist.dat", ios::out);
-          pmsg_m = std::unique_ptr<Inform>(new Inform(0, *f_stream, 0));
-      }
-    */
-}
-
-template <class T, unsigned Dim>
-PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb,
-                                     const std::vector<OpalParticle>& rhs,
-                                     const PartData *ref):
-    R(*(pb->R_p)),
-    ID(*(pb->ID_p)),
-    myNode_m(Ippl::myNode()),
-    nodes_m(Ippl::getNodes()),
-    fixed_grid(false),
-    pbin_m(nullptr),
-    lossDs_m(nullptr),
-    pmsg_m(nullptr),
-    f_stream(nullptr),
-    lowParticleCount_m(false),
-    reference(ref),
-    unit_state_(units),
-    stateOfLastBoundP_(unitless),
-    moments_m(),
-    dt_m(0.0),
-    t_m(0.0),
-    eKin_m(0.0),
-    dE_m(0.0),
-    spos_m(0.0),
-    globalMeanR_m(Vector_t(0.0, 0.0, 0.0)),
-    globalToLocalQuaternion_m(Quaternion_t(1.0, 0.0, 0.0, 0.0)),
-    rmax_m(0.0),
-    rmin_m(0.0),
-    rrms_m(0.0),
-    prms_m(0.0),
-    rmean_m(0.0),
-    pmean_m(0.0),
-    eps_m(0.0),
-    eps_norm_m(0.0),
-    rprms_m(0.0),
-    Dx_m(0.0),
-    Dy_m(0.0),
-    DDx_m(0.0),
-    DDy_m(0.0),
-    hr_m(-1.0),
-    nr_m(0),
-    fs_m(nullptr),
-    couplingConstant_m(0.0),
-    qi_m(0.0),
-    distDump_m(0),
-    fieldDBGStep_m(0),
-    dh_m(1e-12),
-    tEmission_m(0.0),
-    bingamma_m(nullptr),
-    binemitted_m(nullptr),
-    lPath_m(0.0),
-    stepsPerTurn_m(0),
-    localTrackStep_m(0),
-    globalTrackStep_m(0),
-    numBunch_m(1),
-    SteptoLastInj_m(0),
-    globalPartPerNode_m(nullptr),
-    dist_m(nullptr),
-    dcBeam_m(false),
-    periodLength_m(Physics::c / 1e9),
-    pbase(pb)
-{
-
-}
-
-
-template <class T, unsigned Dim>
-PartBunchBase<T, Dim>::PartBunchBase(const PartBunchBase<T, Dim>& rhs):
-    R(rhs.R),
-    ID(rhs.ID),
-    myNode_m(Ippl::myNode()),
-    nodes_m(Ippl::getNodes()),
-    fixed_grid(rhs.fixed_grid),
-    pbin_m(nullptr),
-    lossDs_m(nullptr),
-    pmsg_m(nullptr),
-    f_stream(nullptr),
-    lowParticleCount_m(rhs.lowParticleCount_m),
-    reference(rhs.reference),
-    unit_state_(rhs.unit_state_),
-    stateOfLastBoundP_(rhs.stateOfLastBoundP_),
-    moments_m(rhs.moments_m),
-    dt_m(rhs.dt_m),
-    t_m(rhs.t_m),
-    eKin_m(rhs.eKin_m),
-    dE_m(rhs.dE_m),
-    spos_m(0.0),
-    globalMeanR_m(Vector_t(0.0, 0.0, 0.0)),
-    globalToLocalQuaternion_m(Quaternion_t(1.0, 0.0, 0.0, 0.0)),
-    rmax_m(rhs.rmax_m),
-    rmin_m(rhs.rmin_m),
-    rrms_m(rhs.rrms_m),
-    prms_m(rhs.prms_m),
-    rmean_m(rhs.rmean_m),
-    pmean_m(rhs.pmean_m),
-    eps_m(rhs.eps_m),
-    eps_norm_m(rhs.eps_norm_m),
-    rprms_m(rhs.rprms_m),
-    Dx_m(rhs.Dx_m),
-    Dy_m(rhs.Dy_m),
-    DDx_m(rhs.DDx_m),
-    DDy_m(rhs.DDy_m),
-    hr_m(rhs.hr_m),
-    nr_m(rhs.nr_m),
-    fs_m(nullptr),
-    couplingConstant_m(rhs.couplingConstant_m),
-    qi_m(rhs.qi_m),
-    distDump_m(rhs.distDump_m),
-    fieldDBGStep_m(rhs.fieldDBGStep_m),
-    dh_m(rhs.dh_m),
-    tEmission_m(rhs.tEmission_m),
-    bingamma_m(nullptr),
-    binemitted_m(nullptr),
-    lPath_m(rhs.lPath_m),
-    stepsPerTurn_m(rhs.stepsPerTurn_m),
-    localTrackStep_m(rhs.localTrackStep_m),
-    globalTrackStep_m(rhs.globalTrackStep_m),
-    numBunch_m(rhs.numBunch_m),
-    SteptoLastInj_m(rhs.SteptoLastInj_m),
-    globalPartPerNode_m(nullptr),
-    dist_m(nullptr),
-    dcBeam_m(rhs.dcBeam_m),
-    periodLength_m(rhs.periodLength_m),
-    pbase(rhs.pbase)
+    : PartBunchBase(pb)
 {
+    reference = ref;
 }
 
-
-// template <class T, unsigned Dim>
-// AbstractParticle<T, Dim>* PartBunchBase<T, Dim>::getParticleBase() {
-//     return pbase;
-// }
-//
-//
-// template <class T, unsigned Dim>
-// const AbstractParticle<T, Dim>* PartBunchBase<T, Dim>::getParticleBase() const {
-//     return pbase;
-// }
-
-
 /*
  * Bunch common member functions
  */
@@ -1900,12 +1620,12 @@ double PartBunchBase<T, Dim>::calcMeanPhi() {
     for(int ii = 0; ii < emittedBins; ii++) {
         phi[ii] = calculateAngle(px[ii], py[ii]);
         meanPhi += phi[ii];
-        INFOMSG("Bin " << ii  << "  mean phi = " << phi[ii] * 180.0 / pi - 90.0 << "[degree]" << endl);
+        INFOMSG("Bin " << ii  << "  mean phi = " << phi[ii] * 180.0 / Physics::pi - 90.0 << "[degree]" << endl);
     }
 
     meanPhi /= emittedBins;
 
-    INFOMSG("mean phi of all particles " <<  meanPhi * 180.0 / pi - 90.0 << "[degree]" << endl);
+    INFOMSG("mean phi of all particles " <<  meanPhi * 180.0 / Physics::pi - 90.0 << "[degree]" << endl);
 
     return meanPhi;
 }
@@ -2472,12 +2192,12 @@ const ParticleLayout<T, Dim>& PartBunchBase<T, Dim>::getLayout() const {
 }
 
 template <class T, unsigned Dim>
-bool PartBunchBase<T, Dim>::getUpdateFlag(UpdateFlags f) const {
+bool PartBunchBase<T, Dim>::getUpdateFlag(UpdateFlags_t f) const {
     return pbase->getUpdateFlag(f);
 }
 
 template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::setUpdateFlag(UpdateFlags f, bool val) {
+void PartBunchBase<T, Dim>::setUpdateFlag(UpdateFlags_t f, bool val) {
     pbase->setUpdateFlag(f, val);
 }
 

src/Classic/Structure/PeakFinder.cpp

@@ -7,6 +7,7 @@
 #include "AbstractObjects/OpalData.h"
 #include "Ippl.h"
 
+#include "Algorithms/PartBunchBase.h"
 
 PeakFinder::PeakFinder(std::string elem):
     element_m(elem), nBins_m(10), binWidth_m(1.0 /*mm*/),
@@ -28,11 +29,19 @@ void PeakFinder::save() {
     
     createHistogram_m();
 
-    bool found = findPeaks(smoothingNumber_m,
-                           minArea_m,
-                           minFractionalArea_m,
-                           minAreaAboveNoise_m,
-                           minSlope_m);
+    bool found = false;
+    
+    std::size_t nParticles = (OpalData::getInstance()->getPartBunch())->getTotalNum();
+    
+    if ( Ippl::getNodes() == 1 && nParticles == 1) {
+        found = findPeaks();
+    } else {
+        found = findPeaks(smoothingNumber_m,
+                          minArea_m,
+                          minFractionalArea_m,
+                          minAreaAboveNoise_m,
+                          minSlope_m);
+    }    
     
     if ( found ) {
         fn_m = element_m + std::string(".peaks");
@@ -59,6 +68,13 @@ void PeakFinder::save() {
 }
 
 
+bool PeakFinder::findPeaks() {
+    for (const auto &radius : radius_m)
+        peakRadii_m.push_back(radius);
+    return !peakRadii_m.empty();
+}
+
+
 bool PeakFinder::findPeaks(int smoothingNumber,
                            double minAreaFactor,
                            double minFractionalAreaFactor,

src/Classic/Structure/PeakFinder.h

@@ -54,6 +54,12 @@ public:
                    double minFractionalArea,
                    double minAreaAboveNoise,
                    double minSlope);
+    
+    
+    /**
+     * Single particle peak finder.
+     */
+    bool findPeaks();
 
 private:
     

src/Optimize/OpalSimulation.cpp

@@ -147,6 +147,24 @@ bool OpalSimulation::hasResultsAvailable() {
 }
 
 
+void OpalSimulation::createSymlink_m(const std::string& path) {
+    struct dirent **files;
+    int count = scandir(path.c_str(), &files, 0, alphasort);
+
+    for(int i=0; i<count; i++) {
+	if (files[i]->d_name == std::string(".") ||
+	    files[i]->d_name == std::string("..")) continue;
+	std::string source = path + "/" + files[i]->d_name;
+	std::string target = simulationDirName_ + '/' + files[i]->d_name;
+	int err = symlink(source.c_str(), target.c_str());
+	if (err != 0) {
+	    throw OptPilotException("OpalSimulation::createSymlink()",
+				    "Cannot create symbolic link '" + source + "' to " +
+				    target);
+	}
+    }
+}
+
 void OpalSimulation::setupSimulation() {
     namespace fs = boost::filesystem;
 
@@ -170,31 +188,18 @@ void OpalSimulation::setupSimulation() {
                              simulationName_ + ".in";
         gs_->writeInputFile(infile);
 
-        // linking fieldmaps
+        // linking fieldmaps + distributions
         if(getenv("FIELDMAPS") == NULL) {
             throw OptPilotException("OpalSimulation::OpalSimulation",
                 "Environment variable FIELDMAPS not defined!");
         }
+
         std::string fieldmapPath = getenv("FIELDMAPS");
+	this->createSymlink_m(fieldmapPath);
 
-        struct dirent **files;
-        int count = scandir(fieldmapPath.c_str(), &files, 0, alphasort);
-
-        for(int i=0; i<count; i++) {
-            if (files[i]->d_name == std::string(".") ||
-                files[i]->d_name == std::string("..")) continue;
-            std::string source = fieldmapPath + "/" + files[i]->d_name;
-            std::string target = simulationDirName_ + '/' + files[i]->d_name;
-	    int err = symlink(source.c_str(), target.c_str());
-	    if (err != 0) {
-	      // FIXME properly handle error
-	      std::cout << "Cannot symlink fieldmap "
-			<< source.c_str() << " to "
-			<< target.c_str() << " error no " << err << std::endl;
-	      std::cout << "fieldmapPath " << fieldmapPath << " i= " << i << std::endl;
-	      std::cout << "target       " << simulationDirName_ + '/' + files[i]->d_name << std::endl;
-
-	    }
+        if ( getenv("DISTRIBUTIONS") != NULL ) {
+	    std::string distPath = getenv("DISTRIBUTIONS");
+	    this->createSymlink_m(distPath);
 	}
     }
 
@@ -504,4 +509,4 @@ void OpalSimulation::cleanUp() {
     } catch(...) {
         std::cout << "Can't remove directory '" << simulationDirName_ << "'" << std::endl;
     }
-}
\ No newline at end of file
+}

src/Optimize/OpalSimulation.h

@@ -103,6 +103,10 @@ private:
 
     /// check if we already have simulated the current set of design vars
     bool hasResultsAvailable();
+    
+    /// create symbolic links
+    void createSymlink_m(const std::string& path);
+
     /// create directories, input files, fieldmaps...
     void setupSimulation();
 
@@ -117,4 +121,4 @@ private:
     void restoreOut();
 };
 
-#endif
\ No newline at end of file
+#endif

src/Optimize/OptimizeCmd.cpp

@@ -33,6 +33,7 @@
 #include "Expression/SDDSVariable.h"
 #include "Expression/RadialPeak.h"
 #include "Expression/MaxNormRadialPeak.h"
+#include "Expression/NumberOfPeaks.h"
 #include "Expression/SumErrSqRadialPeak.h"
 #include "Expression/ProbeVariable.h"
 
@@ -73,6 +74,7 @@ namespace {
         SIMTMPDIR,
         TEMPLATEDIR,
         FIELDMAPDIR,
+        DISTDIR,
         SIZE
     };
 }
@@ -132,7 +134,8 @@ OptimizeCmd::OptimizeCmd():
         ("TEMPLATEDIR", "Directory where templates are stored");
     itsAttr[FIELDMAPDIR] = Attributes::makeString
         ("FIELDMAPDIR", "Directory where field maps are stored");
-
+    itsAttr[DISTDIR] = Attributes::makeString
+        ("DISTDIR", "Directory where distributions are stored", "");
     registerOwnership(AttributeHandler::COMMAND);
 }
 
@@ -195,7 +198,11 @@ void OptimizeCmd::execute() {
     ff = MaxNormRadialPeak();
     funcs.insert(std::pair<std::string, client::function::type>
                  ("maxNormRadialPeak", ff));
-
+    
+    ff = NumberOfPeaks();
+    funcs.insert(std::pair<std::string, client::function::type>
+            ("numberOfPeaks", ff));
+    
     ff = SumErrSqRadialPeak();
     funcs.insert(std::pair<std::string, client::function::type>
                  ("sumErrSqRadialPeak", ff));
@@ -314,6 +321,18 @@ void OptimizeCmd::execute() {
         setenv("FIELDMAPS", dir.c_str(), 1);
     }
 
+    if (Attributes::getString(itsAttr[DISTDIR]) != "") {
+        fs::path dir(Attributes::getString(itsAttr[DISTDIR]));
+        if (dir.is_relative()) {
+            fs::path path = fs::path(std::string(getenv("PWD")));
+            path /= dir;
+            dir = path;
+        }
+
+        setenv("DISTRIBUTIONS", dir.c_str(), 1);
+    }
+
+
     *gmsg << endl;
     for (size_t i = 0; i < arguments.size(); ++ i) {
         argv.push_back(const_cast<char*>(arguments[i].c_str()));

src/Sample/SampleCmd.cpp

@@ -46,6 +46,7 @@ namespace {
         NUMCOWORKERS,
         TEMPLATEDIR,
         FIELDMAPDIR,
+        DISTDIR,
         RASTER,
         SEED,
         SIZE
@@ -73,6 +74,8 @@ SampleCmd::SampleCmd():
         ("TEMPLATEDIR", "Directory where templates are stored");
     itsAttr[FIELDMAPDIR] = Attributes::makeString
         ("FIELDMAPDIR", "Directory where field maps are stored");
+    itsAttr[DISTDIR] = Attributes::makeString
+        ("DISTDIR", "Directory where distributions are stored");
     itsAttr[RASTER] = Attributes::makeBool
         ("RASTER", "Scan full space given by design variables (default: true)", true);
     itsAttr[SEED] = Attributes::makeReal
@@ -275,6 +278,17 @@ void SampleCmd::execute() {
 
         setenv("FIELDMAPS", dir.c_str(), 1);
     }
+    
+    if (Attributes::getString(itsAttr[DISTDIR]) != "") {
+        fs::path dir(Attributes::getString(itsAttr[DISTDIR]));
+        if (dir.is_relative()) {
+            fs::path path = fs::path(std::string(getenv("PWD")));
+            path /= dir;
+            dir = path;
+        }
+
+        setenv("DISTRIBUTIONS", dir.c_str(), 1);
+    }
 
     *gmsg << endl;
     for (size_t i = 0; i < arguments.size(); ++ i) {